In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt
from pycvodes import fpes
from pyodesys.tests._robertson import run_integration
from pyodesys.plotting import info_vlines
from pyodesys.native import native_sys
%matplotlib inline


    

In [ ]:

    

def integrate_and_plot(ax, title=None, info_vlines_kw=None, **kwargs):
    xout, yout, info, systems = run_integration(
        integrator='cvode', record_rhs_xvals=True, record_jac_xvals=True,
        record_order=True, record_fpe=True, **kwargs)
    info_vlines(ax, xout, info, **(info_vlines_kw or {
                'fpes': fpes, 'post_proc': (lambda x: np.exp(x)) if kwargs.get('logt', False) else None}))
    for idx, name in enumerate('ABC'):
        ax.plot(xout, yout[:, idx], label=name)
    ax.set_title((title or '') + ('%d steps, %d rhs evals, %.4f s CPU t' % (info['n_steps'], info['nfev'], info['time_cpu'])) +
                 ('' if info['success'] else ' - failed!'))
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_ylim([1e-21, 1e14])
    ax.set_xlim([1e-12, 1e19])
    ax.legend()
    return systems


    

In [ ]:

    

kw = dict(nsteps=18000, atol={'A': 1e-8, 'B': 1e-10, 'C': 1e-4}, rtol=1e-6,
          dep_scaling=1e4, zero_conc=1e-24, t0=1e-12, return_on_error=True)
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
integrate_and_plot(axes[0], **kw)
_ = integrate_and_plot(axes[1], first_step=1e-12, **kw)


    

In [ ]:

    

_ = integrate_and_plot(plt.subplot(1, 1, 1), wrapping_class=native_sys['cvode'], **kw)


    

In [ ]:

    

def integrate_and_plot_reduced(**kwargs):
    kw2 = kw.copy()
    kw2.update(kwargs)
    kw2['wrapping_class'] = native_sys['cvode']
    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    systems = []
    for i in range(4):
        title = 'A B C'.split()
        if i > 0:
            title.pop(i - 1)
        systems.append(integrate_and_plot(
                axes.flat[i], reduced=i, title=''.join(title) + ', ', **kw2))
    return systems


    

In [ ]:

    

rlin = integrate_and_plot_reduced()


    

In [ ]:

    

kw_log = dict(atol=1e-7, logc=True, logt=True, dep_scaling=1)
rlog1 = integrate_and_plot_reduced(**kw_log)


    

In [ ]:

    

rloglin = integrate_and_plot_reduced(atol=1e-7, logc=True, logt=False, dep_scaling=1)


    

In [ ]:

    

rlinlog = integrate_and_plot_reduced(atol=1e-7, logc=False, logt=True, dep_scaling=1)


    

In [ ]:

    

rlog2 = integrate_and_plot_reduced(powsimp=True, **kw_log)


    

In [ ]:

    

import sympy as sp
sp.init_printing()


    

In [ ]:

    

rlog1[0][-1].exprs


    

In [ ]:

    

rlog1[0][-1].get_jac()


    

In [ ]:

    

rlog2[0][-1].exprs


    

In [ ]:

    

rlog2[0][-1].get_jac()


    

In [ ]:

    

rlog2[0][-1]._native._written_files


    

In [ ]:

    

def rhs_code_in_odesys(odesys):
    seen_rhs = False
    for line in open(next(filter(lambda n: n.endswith('.cpp'), odesys._native._written_files))):
        if not seen_rhs:
            if 'rhs(' in line:
                seen_rhs = True
            else:
                continue
        if 'nfev++' in line:
            yield '...'
            break        
        yield line


    

In [ ]:

    

print(''.join(rhs_code_in_odesys(rlog1[0][-1])))


    

In [ ]:

    

print(''.join(rhs_code_in_odesys(rlog2[0][-1])))


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: