Skip to content

Commit

Permalink
updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@bjodah
bjodah committed Feb 7, 2025
1 parent 34ed552 commit 635df1c
Show file tree
Hide file tree
Showing 4 changed files with 177 additions and 31 deletions.
16 changes: 10 additions & 6 deletions chempy/kinetics/_native.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,13 +62,18 @@
const int ny = get_ny();
std::vector<double> f(ny);
double tot=0.0;
rhs(x, y, &f[0]);
auto flag_rhs = rhs(x, y, &f[0]);
if (flag_rhs != AnyODE::Status::success) { return AnyODE::Status::unrecoverable_error; }
for (int i=0; i<ny; ++i){
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!
//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));
//fflush(stderr);
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!
}
//tot = cbrt(tot);
//fprintf(stderr, "tot = %12.5g\\n", tot);
//fflush(stderr);
//return AnyODE::Status::unrecoverable_error;
out[0] = tot/ny - m_special_settings[0];
this->nrev++;
return AnyODE::Status::success;
"""

_constr_special_settings = r"""
Expand Down Expand Up @@ -153,8 +158,7 @@ def get_native(rsys, odesys, integrator, skip_keys=(0,), steady_state_root=False
native_code_kw['namespace_override']['p_nroots'] = ' return %d; ' % len(conc_roots)
native_code_kw['namespace_override']['p_roots'] = (
''.join([' out[%(i)d] = y[%(j)d] - m_special_settings[%(i)d];\n' %
dict(i=i, j=odesys.names.index(k)) for i, k in enumerate(conc_roots)]) +
' return AnyODE::Status::success;\n'
dict(i=i, j=odesys.names.index(k)) for i, k in enumerate(conc_roots)])
)
if 'p_constructor' not in ns_extend:
ns_extend['p_constructor'] = []
Expand Down
25 changes: 17 additions & 8 deletions examples/Ammonia.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,10 @@
"metadata": {},
"outputs": [],
"source": [
"substance_names = ['H+', 'OH-', 'NH4+', 'NH3', 'H2O']\n",
"substance_names = ['H+', 'OH-', 'H2O2', 'HO2-', 'H2O']\n",
"subst = {n: Species.from_formula(n) for n in substance_names}\n",
"assert [subst[n].charge for n in substance_names] == [1, -1, 1, 0, 0], \"Charges of substances\"\n",
"print(u'Composition of %s: %s' % (subst['NH3'].unicode_name, subst['NH3'].composition))"
"assert [subst[n].charge for n in substance_names] == [1, -1, 0, -1, 0], \"Charges of substances\"\n",
"print(u'Composition of %s: %s' % (subst['H2O2'].unicode_name, subst['H2O2'].composition))"
]
},
{
Expand All @@ -61,11 +61,11 @@
"metadata": {},
"outputs": [],
"source": [
"init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}\n",
"init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'H2O2': 1.0, 'HO2-': 1e-20, 'H2O': 55.5}\n",
"x0 = [init_conc[k] for k in substance_names]\n",
"H2O_c = init_conc['H2O']\n",
"w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)\n",
"NH4p_pr = Equilibrium({'NH4+': 1}, {'H+': 1, 'NH3': 1}, 10**-9.26)\n",
"NH4p_pr = Equilibrium({'H2O2': 1}, {'H+': 1, 'HO2-': 1}, 10**-11.8)\n",
"equilibria = w_autop, NH4p_pr\n",
"[(k, init_conc[k]) for k in substance_names]"
]
Expand Down Expand Up @@ -97,6 +97,15 @@
"x, sol['success'], sane"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dict(zip(eqsys.substance_names(), x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -300,7 +309,7 @@
"metadata": {},
"outputs": [],
"source": [
"nc=60\n",
"nc=500\n",
"Hp_0 = np.logspace(-3, 0, nc)\n",
"def plot_rref(**kwargs):\n",
" fig, axes = plt.subplots(2, 2, figsize=(16, 6), subplot_kw=dict(xscale='log', yscale='log'))\n",
Expand All @@ -327,7 +336,7 @@
"for col_id in range(len(substance_names)):\n",
" for i in range(1, 4):\n",
" plt.subplot(1, 3, i, xscale='log')\n",
" plt.gca().set_yscale('symlog', linthreshy=1e-14)\n",
" plt.gca().set_yscale('symlog', linthresh=1e-14)\n",
" plt.plot(Hp_0, res_lin[0][0][:, col_id] - res_lin[i][0][:, col_id])\n",
"plt.tight_layout()"
]
Expand Down Expand Up @@ -436,7 +445,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
"version": "3.12.8"
}
},
"nbformat": 4,
Expand Down
11 changes: 10 additions & 1 deletion examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -636,6 +636,15 @@
"## Reduced non-linearity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys; sys.version, sys.executable, np.__version__"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -2640,7 +2649,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
"version": "3.12.8"
}
},
"nbformat": 4,
Expand Down
156 changes: 140 additions & 16 deletions examples/protein_binding_unfolding_4state_model.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,15 @@
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"MassAction.inactive_reactants_modulation = lambda self, variables, backend=math, reaction=None: 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -149,7 +158,8 @@
"rsys = ReactionSystem(\n",
" eq_dis.as_reactions(kb=kinetics_as, new_name='ligand-protein association') +\n",
" eq_u.as_reactions(kb=kinetics_f, new_name='protein folding') +\n",
" (r_agg,), substances, name='4-state CETSA system')"
" (r_agg,), substances, name='4-state CETSA system')\n",
"print(rsys.string(with_param=False))"
]
},
{
Expand Down Expand Up @@ -257,32 +267,120 @@
"def pretty_replace(s, subs=None):\n",
" if subs is None:\n",
" subs = {\n",
" 'Ha_(\\w+)': r'\\\\Delta_{\\1}H^{\\\\neq}',\n",
" 'Sa_(\\w+)': r'\\\\Delta_{\\1}S^{\\\\neq}',\n",
" 'He_(\\w+)': r'\\\\Delta_{\\1}H^\\\\circ',\n",
" 'Se_(\\w+)': r'\\\\Delta_{\\1}S^\\\\circ',\n",
" 'Cp_(\\w+)': r'\\\\Delta_{\\1}\\,C_p',\n",
" 'Tref_(\\w+)': r'T^{\\\\circ}_{\\1}',\n",
" r'Ha_(\\w+)': r'\\\\Delta_{\\\\rm \\1}H^{\\\\neq}',\n",
" r'Sa_(\\w+)': r'\\\\Delta_{\\\\rm \\1}S^{\\\\neq}',\n",
" r'He_(\\w+)': r'\\\\Delta_{\\\\rm \\1}H^\\\\circ',\n",
" r'Se_(\\w+)': r'\\\\Delta_{\\\\rm \\1}S^\\\\circ',\n",
" r'Cp_(\\w+)': r'\\\\Delta_{\\\\rm \\1}\\,C_p',\n",
" r'Tref_(\\w+)': r'T^{\\\\circ}_{\\\\rm \\1}',\n",
" 'k_B': r'k_{\\\\rm B}',\n",
" 'Tmelt': r'T_{\\\\rm m}',\n",
" }\n",
" for pattern, repl in subs.items():\n",
" s = re.sub(pattern, repl, s)\n",
" return s\n",
"\n",
"_created = {}\n",
"def mk_Symbol(key):\n",
" if key in substances:\n",
" arg = substances[key].latex_name\n",
" else:\n",
" arg = pretty_replace(key.replace('temperature', 'T'))\n",
" if key == 'T' or 'temperature' in key or key.startswith(\"Tref\") or key.startswith(\"Tmelt\"):\n",
" positive = True\n",
" else:\n",
" positive = None\n",
" _created[key] = sympy.Symbol(arg, real=True, positive=positive)\n",
" return _created[key]\n",
"\n",
" return sympy.Symbol(arg)\n",
"def smart_factor(expr):\n",
" numer, denom = expr.as_numer_denom()\n",
" if denom != 1:\n",
" return (smart_factor(numer)/denom).expand(deep=False)\n",
" if expr.is_Mul and len(expr.args) > 1:\n",
" return reduce(mul, map(smart_factor, expr.args))\n",
" const, nonconst = expr.as_coeff_Add()\n",
" cand0 = sympy.factor(nonconst) + const\n",
" candidates = {cand0: cand0.count_ops()}\n",
" for fs in expr.free_symbols:\n",
" i, d = expr.as_independent(fs, as_Add=True)\n",
" cand = sympy.factor(i) + sympy.factor(d)\n",
" candidates[cand] = cand.count_ops()\n",
" #print(f\"{candidates=}\") # debug print, remove in production\n",
" result = min(candidates, key=candidates.__getitem__)\n",
" const, terms = result.as_coeff_add()\n",
" if len(terms) > 1:\n",
" return const + reduce(add, map(smart_factor, terms))\n",
" return result\n",
" \n",
"\n",
"def log_simp(e):\n",
" e = sympy.logcombine(e)\n",
" def _log(*args):\n",
" arg, = args\n",
" pdict = arg.as_powers_dict()\n",
" pv = list(pdict.values())\n",
" if len(pv) == 1:\n",
" ret = pv[0]*sympy.log(*pdict.keys())\n",
" elif all(_ == pv[0] for _ in pv[1:]):\n",
" ret = pv[0]*sympy.log(reduce(mul, pdict.keys()))\n",
" else:\n",
" ret = sympy.log(arg)\n",
" #print(f\"{arg=} {pdict=} {ret=}\")\n",
" return ret\n",
" e = e.replace(sympy.log, _log)\n",
" return e\n",
"x,y,z=sympy.symbols('x,y,z',real=True, positive=True)\n",
"_ = log_simp(z*(sympy.log(x/y)-sympy.log(z/y)))\n",
"assert _ == z*sympy.log(x/z)\n",
"assert _.is_Mul\n",
" \n",
"from operator import mul, add\n",
"from functools import reduce\n",
"autosymbols = defaultkeydict(mk_Symbol)\n",
"rnames = {}\n",
"lbl = {\n",
" 'ligand-protein dissociation': r'\\ref{eq:e_NL}, \\mathrm{fwd}',\n",
" 'ligand-protein association': r'\\ref{eq:e_NL}, \\mathrm{rev}',\n",
" 'protein unfolding': r'\\ref{eq:e_N__U}, \\mathrm{fwd}',\n",
" 'protein folding' : r'\\ref{eq:e_N__U}, \\mathrm{rev}',\n",
" 'protein aggregation': r'\\ref{eq:r_U}'\n",
"}\n",
"lstrs = []\n",
"for rxn in rsys.rxns:\n",
" rnames[rxn.name] = rxn.name.replace(' ', '~').replace('-','-')\n",
" rate_expr_str = sympy.latex(rxn.rate_expr()(autosymbols, backend=sympy, reaction=rxn))\n",
" lstr = r'$r(\\mathrm{%s}) = %s$' % (rnames[rxn.name], rate_expr_str)\n",
" display(Latex(lstr))"
" e1 = rxn.rate_expr()(autosymbols, backend=sympy, reaction=rxn)\n",
" if 'Se_u' in _created:\n",
" Tm = autosymbols['Tmelt']\n",
" dCp = autosymbols['Cp_u']\n",
" T0 = autosymbols['Tref_u']\n",
" _Se_u = autosymbols['He_u']/Tm + (Tm-T0)*dCp/Tm - dCp*sympy.log(Tm/T0)\n",
" e1=e1.subs(_created['Se_u'], _Se_u)\n",
" e2 = sympy.powsimp(e1.expand())\n",
" e3 = e2.replace(sympy.exp, lambda arg: sympy.exp(smart_factor(arg)).expand(deep=False))\n",
" other = [[],[]]\n",
" exps = defaultdict(list)\n",
" for i, numden in enumerate(e3.as_numer_denom()):\n",
" for fact in numden.as_ordered_factors():\n",
" if fact.func is sympy.exp:\n",
" arg, = fact.args\n",
" exnum, exdenom = arg.as_numer_denom()\n",
" exps[exdenom].append((-1)**i * exnum)\n",
" else:\n",
" other[i].append(fact)\n",
" #print(f\"{denom=}\")\n",
" #print(f\"{other=}\")\n",
" #print(f\"{exps=}\")\n",
" e4 = reduce(mul, other[0]+[sympy.exp(sum(v)/k) for k, v in exps.items()])/reduce(mul, other[1])\n",
" 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",
" kB_T__h = autosymbols['k_B']*autosymbols['T']/autosymbols['h']\n",
" e6 = log_simp(e5/kB_T__h)\n",
" assert (e6*kB_T__h - e1).expand().factor().simplify() == 0\n",
" display(Latex(r'$r(\\mathrm{%s}) = %s$' % (rnames[rxn.name], sympy.latex(e6*sympy.UnevaluatedExpr(kB_T__h)))))\n",
" lstrs.append(r\"r_{%s} &= %s\" % (lbl[rxn.name], sympy.latex(e6).replace(r'\\neq', r'\\ddag').replace(\n",
" '] e^{', r'] \\frac{k_{\\rm B} T}{%s} e^{' % {2: r'h c^\\circ', 1: 'h'}[len(rxn.reac)])))\n",
"with open('protein_binding_unfolding_4state_model.txt', 'wt') as ofh:\n",
" ofh.write('\\\\\\\\\\n'.join(lstrs))"
]
},
{
Expand All @@ -291,7 +389,16 @@
"metadata": {},
"outputs": [],
"source": [
"ratexs = [autosymbols['r(\\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns]\n",
"!cat protein_binding_unfolding_4state_model.txt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ratexs = [autosymbols[r'r(\\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns]\n",
"rates = rsys.rates(autosymbols, backend=sympy, ratexs=ratexs)\n",
"for k, v in rates.items():\n",
" display(Latex(r'$\\frac{[%s]}{dt} = %s$' % (k, sympy.latex(v))))"
Expand Down Expand Up @@ -824,7 +931,7 @@
"outputs": [],
"source": [
"native_tLN = NativeCvodeSys.from_other(tsysLN)\n",
"_, _, info_tLN = integrate_and_plot(native_tLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-9)\n",
"_, _, info_tLN = integrate_and_plot(native_tLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-11)\n",
"{k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')}"
]
},
Expand All @@ -834,7 +941,7 @@
"metadata": {},
"outputs": [],
"source": [
"_, _, info_tLN = integrate_and_plot(tsysLN, first_step=1e-9, nsteps=18000, atol=1e-8, rtol=1e-8)\n",
"_, _, info_tLN = integrate_and_plot(tsysLN, first_step=1e-9, nsteps=18000, atol=1e-9, rtol=1e-11)\n",
"{k: info_tLN[k] for k in ('nfev', 'njev', 'n_steps')}"
]
},
Expand All @@ -861,8 +968,25 @@
"metadata": {},
"outputs": [],
"source": [
"print(open(next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))).read())"
"cpp_file = next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))\n",
"print(open(cpp_file).read())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"!clang-format <$cpp_file"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -885,7 +1009,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
"version": "3.12.8"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 635df1c

Please sign in to comment.