<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step2:  Overriding root-finding in the generated C++ code -- Chemical kinetics
This is an advanced notebook exploring how the user can override the root-finding routine in the native C++ code which is exported.
If you only want to do some simple modelling of chemical kinetics you may want look into ChemPy or Cantera (Cantera is more advanced and is the preferred choice when extensive thermodynamic data is available for the reacting species).
Step3:  There are some linking issues with boost's program options in the below (commented) cells.
Step4:  Time to reach steady state
If we define steady state to occur when the change in chemical concentrations is below a certain threshold, then the obtained time will depend on that threshold. Here we investigate how that choice affects the answer (with respect to numerical precision etc.)
<ASSISTANT_TASK:>
Python Code:
from operator import mul
from functools import reduce
import subprocess
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from pyodesys.results import Result
from pyodesys.symbolic import SymbolicSys
from pyodesys.native import native_sys
from pyodesys.native.util import parse_standalone_output
from _chem_kinet import get_odesys
sp.init_printing()
%matplotlib inline
# Fe+3 + SCN- <-> FeSCN+2
stoich_reac, stoich_prod = [(1, 1, 0), (0, 0, 1)], [(0, 0, 1), (1, 1, 0)]
kineticsys = get_odesys(stoich_reac, stoich_prod, 'Fe+3 SCN- FeSCN+2'.split(), SymbolicSys,
                        steady_state_root=True, latex_names=['%s' % s for s in 'Fe^{3+} SCN^- FeSCN^{2+}'.split()])
def integrate_and_plot(odesys, plot=True, **kwargs):
    tend = 2
    result = odesys.integrate(tend, [1e-2, 2e-3, 0], [800, 8], integrator='cvode', **kwargs)
    if plot:
        fig, axes = plt.subplots(1, 2, figsize=(14, 4))
        if result.xout[-1] != tend:
            axes[0].axvline(result.xout[-1], linestyle='--', label='t = %.4f' % result.xout[-1])
        result.plot(ax=axes[0])
        result.plot(ax=axes[1], deriv=True)
        axes[1].set_yscale('symlog', linthreshy=1e-9)
        axes[1].axhline(1e-9, linestyle='--')
        axes[1].axhline(-1e-9, linestyle='--')
        for ax in axes:
            ax.set_xlim([0, tend])
    return result
integrate_and_plot(kineticsys)
integrate_and_plot(kineticsys, atol=1e-14, rtol=1e-14)
integrate_and_plot(kineticsys, atol=1e-14, rtol=1e-14, return_on_root=True)
kineticsys.roots
native_override = {
    'p_nroots':  return 1; ,
    'p_roots': 
    const indextype ny = get_ny();
    std::vector<realtype> f(ny);
    realtype tot=0.0;
    rhs(x, y, &f[0]);
    for (indextype i=0; i<ny; ++i){
        tot += std::min(std::abs(f[i]/m_atol[i]), std::abs(f[i]/y[i]/m_rtol));
    }
    out[0] = tot/ny - m_special_settings[0];
    this->nrev++;
    return AnyODE::Status::success;
}
native_extend={
    'p_constructor': [
        'if (special_settings.size() != 1) std::cerr << "len(special_settings) != 1" << std::endl;'
    ]
}
nativesys = native_sys['cvode'].from_other(
    kineticsys, namespace_override=native_override, namespace_extend=native_extend)
for path in nativesys._native._written_files:
    if path.endswith('.cpp'):
        print(path)
        print('...\n' + ''.join(open(path).readlines()[-20:]))
        print("")
Explanation: Overriding root-finding in the generated C++ code -- Chemical kinetics
This is an advanced notebook exploring how the user can override the root-finding routine in the native C++ code which is exported.
If you only want to do some simple modelling of chemical kinetics you may want look into ChemPy or Cantera (Cantera is more advanced and is the preferred choice when extensive thermodynamic data is available for the reacting species).
End of explanation
#standalone_prog = nativesys.as_standalone('chem_kinet', compile_kwargs=dict(options=['warn', 'pic', 'openmp', 'debug']))
#standalone_prog
#p = subprocess.Popen([standalone_prog, '--return-on-root', '1', '--special-settings', '1000'],
#                     stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
#out, err = p.communicate(input='2 1e-2 2e-3 0 800 8 0 0 0 0 1'.encode('utf-8'))
#retc = p.wait()
#assert retc == 0
#print(err.decode('utf-8'))
#res_sa, = [Result(*args, kineticsys) for args in parse_standalone_output(out.decode('utf-8').split('\n'))]
#res_sa.plot()
Explanation: There are some linking issues with boost's program options in the below (commented) cells.
End of explanation
native = native_sys['cvode'].from_other(kineticsys, namespace_override=native_override)
def plot_tss_conv(factor, tols, ax):
    tol_kw = dict(plot=False, return_on_root=True, nsteps=2000, special_settings=[factor])
    tss = [integrate_and_plot(native, atol=tol, rtol=tol, **tol_kw).xout[-1] for tol in tols]
    ax.semilogx(tols, tss, label=factor)
fig, ax = plt.subplots(figsize=(14, 6))
tols = np.logspace(-15, -10, 200)
for factor in [1e2, 1e3, 1e4, 1.1e4, 1e5, 1e6, 1e7]:
    plot_tss_conv(factor, tols, ax)
ax.legend()
Explanation: Time to reach steady state
If we define steady state to occur when the change in chemical concentrations is below a certain threshold, then the obtained time will depend on that threshold. Here we investigate how that choice affects the answer (with respect to numerical precision etc.)
End of explanation
<END_TASK>