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()
%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))
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 [ ]: