{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# ASM1: true and false positive ramps\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Copyright (C) 2024 Juan Pablo Carbajal\n# Copyright (C) 2024 Mariane Yvonne Schneider\n#\n# This program is free software; you can redistribute it and/or modify\n# it under the terms of the GNU General Public License as published by\n# the Free Software Foundation; either version 3 of the License, or\n# (at your option) any later version.\n#\n# This program is distributed in the hope that it will be useful,\n# but WITHOUT ANY WARRANTY; without even the implied warranty of\n# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n# GNU General Public License for more details.\n#\n# You should have received a copy of the GNU General Public License\n# along with this program. If not, see <http://www.gnu.org/licenses/>.\n\n# Author: Juan Pablo Carbajal <ajuanpi@gmail.com>\n# Author: Mariane Yvonne Schneider <myschneider@meiru.ch>\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom numpy import ma\nfrom sympy import *\n\ntry:\n    import dsafeatures\nexcept ModuleNotFoundError:\n    import sys\n    import os\n\n    sys.path.insert(0, os.path.abspath(\"../..\"))\n\nimport dsafeatures.odemodel as m\nimport dsafeatures.symtools as st\nimport dsafeatures.signal_processing as sp\n\nfrom dsafeatures.printing import *\n\ninit_printing()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Load the model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "asm = m.ODEModel.from_csv(name=\"ASM1\")\n\n# States\nX = asm.states\n# Process rates\nr_sym, r = asm.rates, asm.rates_expr\nJr = r.jacobian(X)\n# Matrix\nM = asm.coef_matrix\n# State derivative\ndotX = st.dXdt(M, r)\n# add areation\nO2, idxO2 = asm.state_by_name(\"S_O2\")\noxyflow = 500.0  # oxygen intake flow\naer = m.oxygen_PC(O2, max_inflow=oxyflow)\ndotX[idxO2] += aer  #\n# Inflection point curve\nipc = st.ddXdt(M[idxO2, :], Jr, dotX)[0] + aer.diff(O2) * dotX[idxO2]\n\n# Parameter values\nparam_values = dict(asm.parameters(and_values=True).set_index(\"sympy\").value)\nparam_values = {k: float(v.subs(param_values)) for k, v in param_values.items()}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "replace parameter symbols with values\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ipc_pv = ipc.subs(param_values)\nipc_states = sorted(tuple(x for x in X if ipc_pv.has(x)), key=str)\n\n# same for ODE system\ndotX_pv = dotX.subs(param_values)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Numerical simulation\nState values initial values\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "states_iv = asm.component[\"states\"][[\"sympy\", \"initial_value\"]].set_index(\"sympy\").loc[list(X), \"initial_value\"]\n# ignore all states not affecting the ODE,\nstates_iv[asm.state_by_symbol_name(\"S_b\")[0]] = 128.94744\nstates_iv[asm.state_by_symbol_name(\"X_cb\")[0]] = 11.48889\nstates_iv[asm.state_by_symbol_name(\"X_h\")[0]] = 1585.25958\nstates_iv[asm.state_by_symbol_name(\"X_a\")[0]] = 0.96299\nstates_iv[O2] = 8.76297\nstates_iv[asm.state_by_symbol_name(\"S_NOx\")[0]] = 0.00235\nNHx, idxNHx = asm.state_by_name(\"S_NHx\")\nstates_iv[NHx] = 0.0065\nstates_iv[asm.state_by_symbol_name(\"S_bN\")[0]] = 0.62963\ncbN = asm.state_by_symbol_name(\"X_cbN\")[0]\nstates_iv[cbN] = 1.01402"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Integrate numerically\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Tend = 7\ndotO2_min = 4e-2\n\nipc_np = lambdify([tuple(X)], ipc_pv, cse=True)\n\ndef ev(t, y):\n    return ipc_np(y)\nev.terminal = 0\nev.direction = -1\n\nsol, dX_np, J_np = sp.integrate_numer(dotX_pv, X, X_iv=states_iv, t_span=(0, Tend), max_step=1e-3,\n                                      events=[ev])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = sol.t\ny = sol.sol(t)\nipc_val = ipc_np(y)\ndx_val = np.zeros_like(t)\nfor k, x_ in enumerate(y.T):\n    dx_val[k] = dX_np(x_).flatten()[idxO2]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "other_states = [x for x in ipc_states if x not in (O2, NHx)] + [cbN]\nfig = plt.figure(layout=\"constrained\")\naxd = fig.subplot_mosaic(\n    [\n        [O2] * 4,\n        other_states[:4],\n        other_states[4:] + ['.'],\n        ['ipc'] * 4,\n    ])\n\nax = axd[O2]\n# add shared axis for NHx\nax.plot(t, y[idxO2])\nclc1 = ax.get_lines()[-1].get_color()\nax.set_ylabel(latex(O2, mode=\"inline\")+' [g$\\mathrm{m}^{-3}$]', color=clc1, fontsize=12)\nax.set_xlabel(\"time [d]\", fontsize=12)\n\naxNHx = ax.twinx()\naxNHx.plot(t, y[idxNHx], color='C1')\nclc = axNHx.get_lines()[-1].get_color()\naxNHx.set_ylabel(latex(NHx, mode=\"inline\")+' [g$\\mathrm{m}^{-3}$]', color=clc, fontsize=12)\naxNHx.spines['right'].set_color(clc)\naxNHx.tick_params(axis='y', colors=clc)\n\nXl_ = list(X)\nfor x in other_states:\n    ax = axd[x]\n    ax.plot(t, y[Xl_.index(x)], label=latex(x, mode=\"inline\"))\n    ax.legend()\n    ax.set_xticklabels([])\n\ntrans = lambda x, n=0.3: np.sign(x) * np.minimum(np.abs(x), 10)\nax = axd['ipc']\nax.axhline(0.0, color='k', alpha=0.3)\nax.plot(t, trans(ipc_val))\nipc_ = ma.masked_where(dx_val <= 4e-2, ipc_val)\nax.plot(t, trans(ipc_), c='r', lw=3, alpha=0.4)\n\nax.set_xlabel(\"time [d]\", fontsize=12)\nax.set_ylabel(\"ramp condition\", fontsize=12)\n\n# Add lines at ramps\nfor t_, y_ in zip(sol.t_events[0], sol.y_events[0]):\n    if dX_np(y_)[idxO2] > dotO2_min:\n        for ax in axd.values():\n            ax.axvline(t_, color='k', alpha=0.3)\n\nif False:  # DEBUG\n    # Spline representation\n    dx_s = sp.make_interp_spline(t, dx_val, k=3)\n    ipc_s = dx_s.derivative()\n\n    ax.plot(t, ipc_s(t), ':')\n    ax_ = ax.twinx()\n    ax_.plot(t, dx_val, label=\"symbolic\")\n    ax_.plot(t, dx_s(t), ':', label=\"numerical\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}