From 635df1cfc0894d4a0abb7c0ebad986e13826d22c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bj=C3=B6rn=20Ingvar=20Dahlgren?= Date: Fri, 7 Feb 2025 18:15:39 +0100 Subject: [PATCH] updates --- chempy/kinetics/_native.py | 16 +- examples/Ammonia.ipynb | 25 ++- ...olution__modified_do_merge_to_master.ipynb | 11 +- ...otein_binding_unfolding_4state_model.ipynb | 156 ++++++++++++++++-- 4 files changed, 177 insertions(+), 31 deletions(-) diff --git a/chempy/kinetics/_native.py b/chempy/kinetics/_native.py index 3516e6b3..1fa20d89 100644 --- a/chempy/kinetics/_native.py +++ b/chempy/kinetics/_native.py @@ -62,13 +62,18 @@ const int ny = get_ny(); std::vector 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; inrev++; - return AnyODE::Status::success; """ _constr_special_settings = r""" @@ -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'] = [] diff --git a/examples/Ammonia.ipynb b/examples/Ammonia.ipynb index 7967a954..09890008 100644 --- a/examples/Ammonia.ipynb +++ b/examples/Ammonia.ipynb @@ -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))" ] }, { @@ -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]" ] @@ -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": {}, @@ -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", @@ -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()" ] @@ -436,7 +445,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.12.8" } }, "nbformat": 4, diff --git a/examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb b/examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb index faa6854f..9c880b64 100644 --- a/examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb +++ b/examples/ammonical_cupric_solution__modified_do_merge_to_master.ipynb @@ -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, @@ -2640,7 +2649,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.12.8" } }, "nbformat": 4, diff --git a/examples/protein_binding_unfolding_4state_model.ipynb b/examples/protein_binding_unfolding_4state_model.ipynb index 80cddc3e..8ed8e2d6 100644 --- a/examples/protein_binding_unfolding_4state_model.ipynb +++ b/examples/protein_binding_unfolding_4state_model.ipynb @@ -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": {}, @@ -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))" ] }, { @@ -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))" ] }, { @@ -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))))" @@ -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')}" ] }, @@ -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')}" ] }, @@ -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": { @@ -885,7 +1009,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.12.8" } }, "nbformat": 4,