#!/usr/bin/env python # coding: utf-8 # In[ ]: import numpy as np import sympy as sp import matplotlib.pyplot as plt from pyodesys.symbolic import PartiallySolvedSystem, ScaledSys, symmetricsys, get_logexp from pyodesys.plotting import plot_result as _plot_res from pycvodes import fpes from chempy import ReactionSystem, Substance from chempy.kinetics.ode import get_odesys sp.init_printing() get_ipython().run_line_magic('matplotlib', 'inline') # In[ ]: rsys = ReactionSystem.from_string( """ A -> B; 0.04 B + C -> A + C; 1e4 2 B -> B + C; 3e7 """, substance_factory=lambda formula: Substance(formula, composition={1: 1}), substances='ABC') # In[ ]: rsys.composition_balance_vectors() # In[ ]: dep_scaling = 1e8 orisys, extra = get_odesys(rsys, description='original') scaledsys = ScaledSys.from_other(orisys, dep_scaling=dep_scaling, description='scaled') # In[ ]: orisys.linear_invariants, orisys.names, orisys.latex_names # In[ ]: indep_end = 1e18 init_conc = {'A': 1, 'B': 0, 'C': 0} integrate_kw = dict(integrator='cvode', atol=1e-6, rtol=1e-6, nsteps=5000, record_rhs_xvals=True, record_jac_xvals=True) def integrate_systems(systems, **kwargs): for k, v in integrate_kw.items(): if k not in kwargs: kwargs[k] = v return [odesys.integrate(indep_end, init_conc, record_order=True, record_fpe=True, return_on_error=True, **kwargs) for odesys in systems] def descr_str(info): keys = 'n_steps nfev time_cpu success'.split() if isinstance(info, dict): nfo = info else: nfo = {} for d in info: for k in keys: if k in nfo: if k == 'success': nfo[k] = nfo[k] and d[k] else: nfo[k] += d[k] else: nfo[k] = d[k] return ' (%d steps, %d rhs evals., %.5g s CPUt)' % tuple([nfo[k] for k in keys[:3]]) + ( '' if nfo['success'] else ' - failed!') def plot_result(description, res, ax=None, vline_post_proc=None, colors=('k', 'r', 'g'), xlim=None): if ax is None: ax = plt.subplot(1, 1, 1) vk = 'steps rhs_xvals jac_xvals'.split() # 'fe_underflow fe_overflow fe_invalid fe_divbyzero' res.plot(xscale='log', yscale='log', ax=ax, c=colors, info_vlines_kw=dict(vline_keys=vk, post_proc=vline_post_proc)) lines = ax.get_lines() for idx, val in enumerate([1 - 0.04*1e-7, 0.04*1e-7]): ax.plot(1e-7, val, marker='o', c=colors[idx]) for idx, val in enumerate([0.2083340149701255e-7, 0.8333360770334713e-13, 0.9999999791665050]): ax.plot(1e11, val, marker='o', c=colors[idx]) ax.legend(loc='best') ax.set_title(description + descr_str(res.info)) if xlim is not None: ax.set_xlim(xlim) def plot_results(systems, results, axes=None, **kwargs): if axes is None: _fig, axes = plt.subplots(2, 2, figsize=(14, 14)) for idx, (odesys, res) in enumerate(zip(systems, results)): plot_result(odesys.description, res, ax=axes.flat[idx], **kwargs) # In[ ]: plot_results([orisys, scaledsys], integrate_systems([orisys, scaledsys], first_step=1e-23), axes=plt.subplots(1, 2, figsize=(14, 7))[1], xlim=(1e-23, indep_end)) # In[ ]: psys = [PartiallySolvedSystem.from_linear_invariants(scaledsys, [k], description=k) for k in 'ABC'] scaled_linear = [scaledsys] + psys [cs.dep for cs in psys] # Switch formulation using roots: # In[ ]: orisys['A'] # In[ ]: rssys_A = PartiallySolvedSystem.from_linear_invariants( scaledsys, ['A'], description='A root finding', roots=[10*scaledsys['A'] - scaledsys['C']]) # In[ ]: res_A = rssys_A.integrate(indep_end, init_conc, return_on_root=True, **integrate_kw) res_C = psys[2].integrate(indep_end - res_A.xout[-1], res_A.yout[-1, :], **integrate_kw) # In[ ]: fig, axes = plt.subplots(1, 3, figsize=(16, 4)) plt_kw = dict(xscale='log', yscale='log') res_A.plot(**plt_kw, ax=axes[0]) res_C.plot(**plt_kw, ax=axes[1]) _plot_res(np.concatenate((res_A.xout, res_C.xout[1:] + res_A.xout[-1])), np.concatenate((res_A.yout, res_C.yout[1:])), ax=axes[2], names='ABC', **plt_kw) for descr, ax, info in zip(['A', 'C', 'A/C'], axes, [res_A.info, res_C.info, (res_A.info, res_C.info)]): ax.legend(loc='best', prop={'size': 11}) ax.set_title(descr + descr_str(info)) ax.set_xlim([1e-9, indep_end]) # Not switching: # In[ ]: linresults = integrate_systems(scaled_linear) # In[ ]: plot_results(scaled_linear, linresults, xlim=(1e-9, indep_end)) # In[ ]: plot_result('scaled A/C ', rssys_A.integrate( indep_end, init_conc, return_on_root=True, **integrate_kw).extend_by_integration( indep_end, odesys=scaled_linear[3], **integrate_kw), xlim=(1e-9, indep_end)) # ## Logarithmic # In[ ]: indep0 = 1e-26 def symmetric_factory(exprs_process_cb=lambda exprs: [sp.powsimp(expr.expand(), force=True) for expr in exprs]): return symmetricsys(get_logexp(1, 1e-20, b2=0), get_logexp(1, indep0, b2=0), exprs_process_cb=exprs_process_cb, check_transforms=False) LogLogSys_together = symmetric_factory() stlogsystems = [LogLogSys_together.from_other(ls, description='stl ' + ls.description) for ls in scaled_linear] # In[ ]: stlogresults = integrate_systems(stlogsystems, nsteps=500) # In[ ]: plot_results(stlogsystems, stlogresults, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0), xlim=[indep0*10, indep_end]) # In[ ]: stlogsystems[0].exprs # In[ ]: unscaled_linear = [orisys] + [PartiallySolvedSystem.from_linear_invariants( orisys, [k], description=k) for k in 'ABC'] utlogsystems = [LogLogSys_together.from_other(ls, description='utl ' + ls.description) for ls in unscaled_linear] utlogresults = integrate_systems(utlogsystems, nsteps=500) plot_results(utlogsystems, utlogresults, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0), xlim=[indep0*10, indep_end]) # In[ ]: utlogsystems[0].exprs # In[ ]: LogLogSys_apart = symmetric_factory(None) salogsystems = [LogLogSys_apart.from_other(ls, description='sal ' + ls.description) for ls in scaled_linear] salogresults = integrate_systems(salogsystems, nsteps=500) plot_results(salogsystems, salogresults, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0), xlim=[indep0*10, indep_end]) # In[ ]: ualogsystems = [LogLogSys_apart.from_other(ls, description='ual ' + ls.description) for ls in unscaled_linear] ualogresults = integrate_systems(ualogsystems, nsteps=500) plot_results(ualogsystems, ualogresults, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0), xlim=[indep0*10, indep_end]) # In[ ]: rsalogsys_A = LogLogSys_apart.from_other(rssys_A, autonomous_interface=True) assert rsalogsys_A.autonomous_interface rsalogsys_A.roots, rsalogsys_A.exprs # In[ ]: res_log_switch = rsalogsys_A.integrate(indep_end, init_conc, return_on_root=True, **integrate_kw) res_log_switch.extend_by_integration(indep_end, odesys=salogsystems[3], return_on_error=True, autonomous=True, **integrate_kw) fig, ax = plt.subplots(1, 1, figsize=(14, 4)) ax.axvline(res_log_switch.xout[res_log_switch.info['root_indices'][0]], ls='--') res_log_switch.plot(ax=ax, xscale='log', yscale='log') ax.legend(loc='best', prop={'size': 11}) ax.set_title('log A/C ' + descr_str(res_log_switch.info)) _ = ax.set_xlim([1e-9, indep_end]) # In[ ]: fig, ax = plt.subplots(1, 1, figsize=(16, 4)) plot_result('log A/C', res_log_switch, ax, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0)) _ = ax.set_xlim([1e-20, 1e12]) # In[ ]: