{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# ASM1: samples in the SO ramp manifold in phase space\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 .\n\n# Author: Juan Pablo Carbajal \n# Author: Mariane Yvonne Schneider \n\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport seaborn as sns\nfrom sympy import *\n\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.levelsample as ls\nfrom dsafeatures.printing import *\nfrom dsafeatures.odemodel import ODEModel\n\ninit_printing()\n\nplt.style.use('figures.mplstyle')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Select model\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"model = ODEModel.from_csv(name='ASM1')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sample the ramp manifold\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"max_O2 = 8.0 # WARNING: if changed delete pickle file\nmax_O2_inflow = 1e3 # WARNING: if changed delete pickle file\nfilepath = Path('example_asm1_SO_ramp_hypersurface.pkl')\nif not filepath.exists():\n points, rng = ls.sample_ramp(model=model, file=filepath, npoints=1024,\n max_O2_inflow=max_O2_inflow,\n S_O2=(1e-9, max_O2),\n append=False, show=False)\nelse:\n points = pd.read_pickle(filepath)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot samples on the ramp manifold\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nhx_T = 'NHx level'\nnT = points[nhx_T].sum()\nlbls = {True: fr\"$S_\\text{{NHx}}$ < T ({nT})\", False: fr\"$S_\\text{{NHx}}$ \u2265 T ({points.shape[0] - nT})\"}\ntmp_ = points.replace({nhx_T: lbls})\n\nhue_order = [lbls[True], lbls[False]]\nfd, fig = ls.plot_samples_ramp(points=tmp_, hue=nhx_T, hue_order=hue_order, show=False, plot_kws={'alpha': 0.1})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for ax_ in fd.axes[-1, :]:\n for ax in (ax_.xaxis, ax_.yaxis):\n lbl = ax.get_label()\n if lbl.get_text() == 'dotO2':\n s = r'$\\dot{S}_\\text{O2}$'\n lbl.set_text(s)\nfig.tight_layout(pad=0.2)\nfig.subplots_adjust(left=0.05, bottom=0.075, right=0.99, top=1.0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for l in fd.legend.get_lines():\n l.set_markersize(16)\n l.set_alpha(1.0)\n\nfd.legend.set_title('')\n\n# Inset with regression and other coloring\nn = len(points.columns) - 2\ninset = fig.add_subplot(n, n, 3)\ninset.set_position([0.55, 0.55, 0.29, 0.39])\n\nXa_lvl = r'$X_\\text{a}$ level'\nbins = [0, 150, 250, 1e6]\ntmp_[Xa_lvl] = pd.cut(tmp_[r'$X_\\text{a}$'], bins,\n labels=[f\"\u2264 {bins[1]}\", f\" [{bins[1]}, {bins[-2]}]\", f\"> {bins[-2]}\"])\n\nsns.kdeplot(tmp_, x='dotO2', y=r'$S_\\text{NHx}$', hue=Xa_lvl, palette='tab10',\n fill=True, cut=0, levels=2, bw_adjust=1.5, ax=inset)\nsns.scatterplot(tmp_, x='dotO2', y=r'$S_\\text{NHx}$', marker='.', color='k',\n alpha=0.25, ax=inset, legend=False)\ninset.axhline(1, c='r')\n\n# # Regression plot\n# args_ = dict(x='dotO2', y=r'$S_\\text{NHx}$', ax=inset, marker='.')\n# sns.regplot(points[points[nhx_T]], **args_)\n# sns.regplot(points[~points[nhx_T]], **args_)\n\ninset.set_xlabel(r'$\\dot{S}_\\text{O2} \\left[g m^{-3} d^{-1}\\right]$')\ninset.set_ylabel(r'$S_\\text{NHx} \\left[g m^{-3}\\right]$')\ninset.set_ylim(-1e-2, inset.get_ylim()[1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot the quality of the ramp conditions\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axs = ls.plot_quality_ramp(points=tmp_, hue=nhx_T, hue_order=hue_order, stat=\"proportion\", bins=25,\n common_norm=False, show=False)\naxs[0].set_xlabel(r'$\\ddot{S}_\\text{O2} \\left[g m^{-3} d^{-2}\\right]$')\naxs[1].set_xlabel(r'$\\dot{S}_\\text{O2} \\left[g m^{-3} d^{-1}\\right]$')"
]
},
{
"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
}