Skip to content

Commit 635df1c

Browse files
bjodahbjodah
committed
updates
1 parent 34ed552 commit 635df1c

File tree

4 files changed

+177
-31
lines changed

4 files changed

+177
-31
lines changed

‎chempy/kinetics/_native.py

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -62,13 +62,18 @@
6262
const int ny = get_ny();
6363
std::vector<double> f(ny);
6464
double tot=0.0;
65-
rhs(x, y, &f[0]);
65+
auto flag_rhs = rhs(x, y, &f[0]);
66+
if (flag_rhs != AnyODE::Status::success) { return AnyODE::Status::unrecoverable_error; }
6667
for (int i=0; i<ny; ++i){
67-
tot += std::min(std::abs(f[i]/m_atol[i]), std::abs(f[i]/y[i]/m_rtol)); // m_atol needs to be of size ny!
68+
//fprintf(stderr, "f[i] = %12.5g y[i] = %12.5g f[i]/(|y[i]|+1e-100) = %12.5g\\n", f[i], y[i], f[i]/(fabs(y[i])+1e-100));
69+
//fflush(stderr);
70+
tot += fmin(fabs(f[i]/(m_atol[i]+1e-100)), fabs(f[i]/(fabs(y[i])+1e-100)/(m_rtol+1e-100))); // m_atol needs to be of size ny!
6871
}
72+
//tot = cbrt(tot);
73+
//fprintf(stderr, "tot = %12.5g\\n", tot);
74+
//fflush(stderr);
75+
//return AnyODE::Status::unrecoverable_error;
6976
out[0] = tot/ny - m_special_settings[0];
70-
this->nrev++;
71-
return AnyODE::Status::success;
7277
"""
7378

7479
_constr_special_settings = r"""
@@ -153,8 +158,7 @@ def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False
153158
native_code_kw['namespace_override']['p_nroots'] = ' return %d; ' % len(conc_roots)
154159
native_code_kw['namespace_override']['p_roots'] = (
155160
''.join([' out[%(i)d] = y[%(j)d] - m_special_settings[%(i)d];\n' %
156-
dict(i=i, j=odesys.names.index(k)) for i, k in enumerate(conc_roots)]) +
157-
' return AnyODE::Status::success;\n'
161+
dict(i=i, j=odesys.names.index(k)) for i, k in enumerate(conc_roots)])
158162
)
159163
if 'p_constructor' not in ns_extend:
160164
ns_extend['p_constructor'] = []

‎examples/Ammonia.ipynb

Lines changed: 17 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -42,10 +42,10 @@
4242
"metadata": {},
4343
"outputs": [],
4444
"source": [
45-
"substance_names = ['H+', 'OH-', 'NH4+', 'NH3', 'H2O']\n",
45+
"substance_names = ['H+', 'OH-', 'H2O2', 'HO2-', 'H2O']\n",
4646
"subst = {n: Species.from_formula(n) for n in substance_names}\n",
47-
"assert [subst[n].charge for n in substance_names] == [1, -1, 1, 0, 0], \"Charges of substances\"\n",
48-
"print(u'Composition of %s: %s' % (subst['NH3'].unicode_name, subst['NH3'].composition))"
47+
"assert [subst[n].charge for n in substance_names] == [1, -1, 0, -1, 0], \"Charges of substances\"\n",
48+
"print(u'Composition of %s: %s' % (subst['H2O2'].unicode_name, subst['H2O2'].composition))"
4949
]
5050
},
5151
{
@@ -61,11 +61,11 @@
6161
"metadata": {},
6262
"outputs": [],
6363
"source": [
64-
"init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}\n",
64+
"init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'H2O2': 1.0, 'HO2-': 1e-20, 'H2O': 55.5}\n",
6565
"x0 = [init_conc[k] for k in substance_names]\n",
6666
"H2O_c = init_conc['H2O']\n",
6767
"w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)\n",
68-
"NH4p_pr = Equilibrium({'NH4+': 1}, {'H+': 1, 'NH3': 1}, 10**-9.26)\n",
68+
"NH4p_pr = Equilibrium({'H2O2': 1}, {'H+': 1, 'HO2-': 1}, 10**-11.8)\n",
6969
"equilibria = w_autop, NH4p_pr\n",
7070
"[(k, init_conc[k]) for k in substance_names]"
7171
]
@@ -97,6 +97,15 @@
9797
"x, sol['success'], sane"
9898
]
9999
},
100+
{
101+
"cell_type": "code",
102+
"execution_count": null,
103+
"metadata": {},
104+
"outputs": [],
105+
"source": [
106+
"dict(zip(eqsys.substance_names(), x))"
107+
]
108+
},
100109
{
101110
"cell_type": "markdown",
102111
"metadata": {},
@@ -300,7 +309,7 @@
300309
"metadata": {},
301310
"outputs": [],
302311
"source": [
303-
"nc=60\n",
312+
"nc=500\n",
304313
"Hp_0 = np.logspace(-3, 0, nc)\n",
305314
"def plot_rref(**kwargs):\n",
306315
" fig, axes = plt.subplots(2, 2, figsize=(16, 6), subplot_kw=dict(xscale='log', yscale='log'))\n",
@@ -327,7 +336,7 @@
327336
"for col_id in range(len(substance_names)):\n",
328337
" for i in range(1, 4):\n",
329338
" plt.subplot(1, 3, i, xscale='log')\n",
330-
" plt.gca().set_yscale('symlog', linthreshy=1e-14)\n",
339+
" plt.gca().set_yscale('symlog', linthresh=1e-14)\n",
331340
" plt.plot(Hp_0, res_lin[0][0][:, col_id] - res_lin[i][0][:, col_id])\n",
332341
"plt.tight_layout()"
333342
]
@@ -436,7 +445,7 @@
436445
"name": "python",
437446
"nbconvert_exporter": "python",
438447
"pygments_lexer": "ipython3",
439-
"version": "3.8.10"
448+
"version": "3.12.8"
440449
}
441450
},
442451
"nbformat": 4,

‎examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -636,6 +636,15 @@
636636
"## Reduced non-linearity"
637637
]
638638
},
639+
{
640+
"cell_type": "code",
641+
"execution_count": null,
642+
"metadata": {},
643+
"outputs": [],
644+
"source": [
645+
"import sys; sys.version, sys.executable, np.__version__"
646+
]
647+
},
639648
{
640649
"cell_type": "code",
641650
"execution_count": null,
@@ -2640,7 +2649,7 @@
26402649
"name": "python",
26412650
"nbconvert_exporter": "python",
26422651
"pygments_lexer": "ipython3",
2643-
"version": "3.10.6"
2652+
"version": "3.12.8"
26442653
}
26452654
},
26462655
"nbformat": 4,

‎examples/protein_binding_unfolding_4state_model.ipynb

Lines changed: 140 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,15 @@
3939
"%matplotlib inline"
4040
]
4141
},
42+
{
43+
"cell_type": "code",
44+
"execution_count": null,
45+
"metadata": {},
46+
"outputs": [],
47+
"source": [
48+
"MassAction.inactive_reactants_modulation = lambda self, variables, backend=math, reaction=None: 1"
49+
]
50+
},
4251
{
4352
"cell_type": "markdown",
4453
"metadata": {},
@@ -149,7 +158,8 @@
149158
"rsys = ReactionSystem(\n",
150159
" eq_dis.as_reactions(kb=kinetics_as, new_name='ligand-protein association') +\n",
151160
" eq_u.as_reactions(kb=kinetics_f, new_name='protein folding') +\n",
152-
" (r_agg,), substances, name='4-state CETSA system')"
161+
" (r_agg,), substances, name='4-state CETSA system')\n",
162+
"print(rsys.string(with_param=False))"
153163
]
154164
},
155165
{
@@ -257,32 +267,120 @@
257267
"def pretty_replace(s, subs=None):\n",
258268
" if subs is None:\n",
259269
" subs = {\n",
260-
" 'Ha_(\\w+)': r'\\\\Delta_{\\1}H^{\\\\neq}',\n",
261-
" 'Sa_(\\w+)': r'\\\\Delta_{\\1}S^{\\\\neq}',\n",
262-
" 'He_(\\w+)': r'\\\\Delta_{\\1}H^\\\\circ',\n",
263-
" 'Se_(\\w+)': r'\\\\Delta_{\\1}S^\\\\circ',\n",
264-
" 'Cp_(\\w+)': r'\\\\Delta_{\\1}\\,C_p',\n",
265-
" 'Tref_(\\w+)': r'T^{\\\\circ}_{\\1}',\n",
270+
" r'Ha_(\\w+)': r'\\\\Delta_{\\\\rm \\1}H^{\\\\neq}',\n",
271+
" r'Sa_(\\w+)': r'\\\\Delta_{\\\\rm \\1}S^{\\\\neq}',\n",
272+
" r'He_(\\w+)': r'\\\\Delta_{\\\\rm \\1}H^\\\\circ',\n",
273+
" r'Se_(\\w+)': r'\\\\Delta_{\\\\rm \\1}S^\\\\circ',\n",
274+
" r'Cp_(\\w+)': r'\\\\Delta_{\\\\rm \\1}\\,C_p',\n",
275+
" r'Tref_(\\w+)': r'T^{\\\\circ}_{\\\\rm \\1}',\n",
276+
" 'k_B': r'k_{\\\\rm B}',\n",
277+
" 'Tmelt': r'T_{\\\\rm m}',\n",
266278
" }\n",
267279
" for pattern, repl in subs.items():\n",
268280
" s = re.sub(pattern, repl, s)\n",
269281
" return s\n",
270282
"\n",
283+
"_created = {}\n",
271284
"def mk_Symbol(key):\n",
272285
" if key in substances:\n",
273286
" arg = substances[key].latex_name\n",
274287
" else:\n",
275288
" arg = pretty_replace(key.replace('temperature', 'T'))\n",
289+
" if key == 'T' or 'temperature' in key or key.startswith(\"Tref\") or key.startswith(\"Tmelt\"):\n",
290+
" positive = True\n",
291+
" else:\n",
292+
" positive = None\n",
293+
" _created[key] = sympy.Symbol(arg, real=True, positive=positive)\n",
294+
" return _created[key]\n",
276295
"\n",
277-
" return sympy.Symbol(arg)\n",
296+
"def smart_factor(expr):\n",
297+
" numer, denom = expr.as_numer_denom()\n",
298+
" if denom != 1:\n",
299+
" return (smart_factor(numer)/denom).expand(deep=False)\n",
300+
" if expr.is_Mul and len(expr.args) > 1:\n",
301+
" return reduce(mul, map(smart_factor, expr.args))\n",
302+
" const, nonconst = expr.as_coeff_Add()\n",
303+
" cand0 = sympy.factor(nonconst) + const\n",
304+
" candidates = {cand0: cand0.count_ops()}\n",
305+
" for fs in expr.free_symbols:\n",
306+
" i, d = expr.as_independent(fs, as_Add=True)\n",
307+
" cand = sympy.factor(i) + sympy.factor(d)\n",
308+
" candidates[cand] = cand.count_ops()\n",
309+
" #print(f\"{candidates=}\") # debug print, remove in production\n",
310+
" result = min(candidates, key=candidates.__getitem__)\n",
311+
" const, terms = result.as_coeff_add()\n",
312+
" if len(terms) > 1:\n",
313+
" return const + reduce(add, map(smart_factor, terms))\n",
314+
" return result\n",
315+
" \n",
278316
"\n",
317+
"def log_simp(e):\n",
318+
" e = sympy.logcombine(e)\n",
319+
" def _log(*args):\n",
320+
" arg, = args\n",
321+
" pdict = arg.as_powers_dict()\n",
322+
" pv = list(pdict.values())\n",
323+
" if len(pv) == 1:\n",
324+
" ret = pv[0]*sympy.log(*pdict.keys())\n",
325+
" elif all(_ == pv[0] for _ in pv[1:]):\n",
326+
" ret = pv[0]*sympy.log(reduce(mul, pdict.keys()))\n",
327+
" else:\n",
328+
" ret = sympy.log(arg)\n",
329+
" #print(f\"{arg=} {pdict=} {ret=}\")\n",
330+
" return ret\n",
331+
" e = e.replace(sympy.log, _log)\n",
332+
" return e\n",
333+
"x,y,z=sympy.symbols('x,y,z',real=True, positive=True)\n",
334+
"_ = log_simp(z*(sympy.log(x/y)-sympy.log(z/y)))\n",
335+
"assert _ == z*sympy.log(x/z)\n",
336+
"assert _.is_Mul\n",
337+
" \n",
338+
"from operator import mul, add\n",
339+
"from functools import reduce\n",
279340
"autosymbols = defaultkeydict(mk_Symbol)\n",
280341
"rnames = {}\n",
342+
"lbl = {\n",
343+
" 'ligand-protein dissociation': r'\\ref{eq:e_NL}, \\mathrm{fwd}',\n",
344+
" 'ligand-protein association': r'\\ref{eq:e_NL}, \\mathrm{rev}',\n",
345+
" 'protein unfolding': r'\\ref{eq:e_N__U}, \\mathrm{fwd}',\n",
346+
" 'protein folding' : r'\\ref{eq:e_N__U}, \\mathrm{rev}',\n",
347+
" 'protein aggregation': r'\\ref{eq:r_U}'\n",
348+
"}\n",
349+
"lstrs = []\n",
281350
"for rxn in rsys.rxns:\n",
282351
" rnames[rxn.name] = rxn.name.replace(' ', '~').replace('-','-')\n",
283-
" rate_expr_str = sympy.latex(rxn.rate_expr()(autosymbols, backend=sympy, reaction=rxn))\n",
284-
" lstr = r'$r(\\mathrm{%s}) = %s$' % (rnames[rxn.name], rate_expr_str)\n",
285-
" display(Latex(lstr))"
352+
" e1 = rxn.rate_expr()(autosymbols, backend=sympy, reaction=rxn)\n",
353+
" if 'Se_u' in _created:\n",
354+
" Tm = autosymbols['Tmelt']\n",
355+
" dCp = autosymbols['Cp_u']\n",
356+
" T0 = autosymbols['Tref_u']\n",
357+
" _Se_u = autosymbols['He_u']/Tm + (Tm-T0)*dCp/Tm - dCp*sympy.log(Tm/T0)\n",
358+
" e1=e1.subs(_created['Se_u'], _Se_u)\n",
359+
" e2 = sympy.powsimp(e1.expand())\n",
360+
" e3 = e2.replace(sympy.exp, lambda arg: sympy.exp(smart_factor(arg)).expand(deep=False))\n",
361+
" other = [[],[]]\n",
362+
" exps = defaultdict(list)\n",
363+
" for i, numden in enumerate(e3.as_numer_denom()):\n",
364+
" for fact in numden.as_ordered_factors():\n",
365+
" if fact.func is sympy.exp:\n",
366+
" arg, = fact.args\n",
367+
" exnum, exdenom = arg.as_numer_denom()\n",
368+
" exps[exdenom].append((-1)**i * exnum)\n",
369+
" else:\n",
370+
" other[i].append(fact)\n",
371+
" #print(f\"{denom=}\")\n",
372+
" #print(f\"{other=}\")\n",
373+
" #print(f\"{exps=}\")\n",
374+
" e4 = reduce(mul, other[0]+[sympy.exp(sum(v)/k) for k, v in exps.items()])/reduce(mul, other[1])\n",
375+
" e5 = e4.replace(lambda arg: arg.func is sympy.exp, lambda arg: sympy.exp(smart_factor(arg.args[0].as_numer_denom()[0])/arg.args[0].as_numer_denom()[1]))\n",
376+
" kB_T__h = autosymbols['k_B']*autosymbols['T']/autosymbols['h']\n",
377+
" e6 = log_simp(e5/kB_T__h)\n",
378+
" assert (e6*kB_T__h - e1).expand().factor().simplify() == 0\n",
379+
" display(Latex(r'$r(\\mathrm{%s}) = %s$' % (rnames[rxn.name], sympy.latex(e6*sympy.UnevaluatedExpr(kB_T__h)))))\n",
380+
" lstrs.append(r\"r_{%s} &= %s\" % (lbl[rxn.name], sympy.latex(e6).replace(r'\\neq', r'\\ddag').replace(\n",
381+
" '] e^{', r'] \\frac{k_{\\rm B} T}{%s} e^{' % {2: r'h c^\\circ', 1: 'h'}[len(rxn.reac)])))\n",
382+
"with open('protein_binding_unfolding_4state_model.txt', 'wt') as ofh:\n",
383+
" ofh.write('\\\\\\\\\\n'.join(lstrs))"
286384
]
287385
},
288386
{
@@ -291,7 +389,16 @@
291389
"metadata": {},
292390
"outputs": [],
293391
"source": [
294-
"ratexs = [autosymbols['r(\\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns]\n",
392+
"!cat protein_binding_unfolding_4state_model.txt"
393+
]
394+
},
395+
{
396+
"cell_type": "code",
397+
"execution_count": null,
398+
"metadata": {},
399+
"outputs": [],
400+
"source": [
401+
"ratexs = [autosymbols[r'r(\\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns]\n",
295402
"rates = rsys.rates(autosymbols, backend=sympy, ratexs=ratexs)\n",
296403
"for k, v in rates.items():\n",
297404
" display(Latex(r'$\\frac{[%s]}{dt} = %s$' % (k, sympy.latex(v))))"
@@ -824,7 +931,7 @@
824931
"outputs": [],
825932
"source": [
826933
"native_tLN = NativeCvodeSys.from_other(tsysLN)\n",
827-
"_, _, info_tLN = integrate_and_plot(native_tLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-9)\n",
934+
"_, _, info_tLN = integrate_and_plot(native_tLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-11)\n",
828935
"{k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')}"
829936
]
830937
},
@@ -834,7 +941,7 @@
834941
"metadata": {},
835942
"outputs": [],
836943
"source": [
837-
"_, _, info_tLN = integrate_and_plot(tsysLN, first_step=1e-9, nsteps=18000, atol=1e-8, rtol=1e-8)\n",
944+
"_, _, info_tLN = integrate_and_plot(tsysLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-11)\n",
838945
"{k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')}"
839946
]
840947
},
@@ -861,8 +968,25 @@
861968
"metadata": {},
862969
"outputs": [],
863970
"source": [
864-
"print(open(next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))).read())"
971+
"cpp_file = next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))\n",
972+
"print(open(cpp_file).read())"
973+
]
974+
},
975+
{
976+
"cell_type": "code",
977+
"execution_count": null,
978+
"metadata": {},
979+
"outputs": [],
980+
"source": [
981+
"!clang-format <$cpp_file"
865982
]
983+
},
984+
{
985+
"cell_type": "code",
986+
"execution_count": null,
987+
"metadata": {},
988+
"outputs": [],
989+
"source": []
866990
}
867991
],
868992
"metadata": {
@@ -885,7 +1009,7 @@
8851009
"name": "python",
8861010
"nbconvert_exporter": "python",
8871011
"pygments_lexer": "ipython3",
888-
"version": "3.8.10"
1012+
"version": "3.12.8"
8891013
}
8901014
},
8911015
"nbformat": 4,

0 commit comments

Comments
 (0)