TUTORIAL_parameter_evolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter evolution\n",
    "\n",
    "Experimental conditions might change during an experiment. For example, the concentration of an oxidant could increase or decrease after a certain time. There could be a correction factor which needs to be applied to the gas phase concentration of an oxidant in order to define the near-surface concentration. The temperature could change - the list goes on. \n",
    "\n",
    "Some parameters may depend on others via some empirical relationship. \n",
    "\n",
    "MultilayerPy facilitates the inclusion model parameter evolution. An additional set of varying parameters associated with these can also be included. To do this, the user is required to supply a function which takes the parameter dictionary as an input (along with any other additional varying parameters) and returns the modified parameter dictionary as an output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the model used in the crash course notebook\n",
    "\n",
    "# importing the necessaries\n",
    "import numpy as np\n",
    "import multilayerpy\n",
    "import multilayerpy.build as build \n",
    "import multilayerpy.simulate as simulate\n",
    "import multilayerpy.optimize as optimize\n",
    "\n",
    "# get the version of multilayerpy\n",
    "print(multilayerpy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the ModelType class\n",
    "from multilayerpy.build import ModelType\n",
    "\n",
    "# import the ReactionScheme class\n",
    "from multilayerpy.build import ReactionScheme\n",
    "\n",
    "# define the model type (KM-SUB in this case) and geometry (spherical or film)\n",
    "mod_type = ModelType('km-sub','spherical')\n",
    "\n",
    "# build the reaction tuple list, in this case only 1 tuple in the list (for 1 reaction)\n",
    "# component 1 (oleic acid) reacts with component 2 (ozone)\n",
    "reaction_tuple_list = [(1,2)]\n",
    "\n",
    "# build the product tuple list, only component 3 (products) is a product\n",
    "# a tuple with a single value inside is defined (value,)\n",
    "product_tuple_list = [(3,)]\n",
    "\n",
    "# now construct the reaction scheme\n",
    "# we can give it a name and define the nuber of components as below\n",
    "reaction_scheme = ReactionScheme(mod_type,name='Oleic acid ozonolysis',\n",
    "                                                   reactants=reaction_tuple_list,\n",
    "                                                products=product_tuple_list)\n",
    "\n",
    "# let's print out a representation of the reaction scheme\n",
    "reaction_scheme.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ModelComponent class\n",
    "from multilayerpy.build import ModelComponent\n",
    "\n",
    "# making model components\n",
    "\n",
    "# oleic acid\n",
    "OA = ModelComponent(1,reaction_scheme,name='Oleic acid')\n",
    "\n",
    "# ozone, declare that it is a gas\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone') \n",
    "\n",
    "# products\n",
    "prod = ModelComponent(3,reaction_scheme, name='Reaction products')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':OA,\n",
    "                        '2':O3,\n",
    "                        '3':prod}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import DiffusionRegime class\n",
    "from multilayerpy.build import DiffusionRegime\n",
    "\n",
    "# making the diffusion dictionary\n",
    "diff_dict = None\n",
    "\n",
    "# make diffusion regime\n",
    "diff_regime = DiffusionRegime(mod_type,model_components_dict,diff_dict=diff_dict)\n",
    "\n",
    "# call it to build diffusion code ready for the builder\n",
    "diff_regime()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ModelBuilder class\n",
    "from multilayerpy.build import ModelBuilder\n",
    "\n",
    "# create the model object\n",
    "model = ModelBuilder(reaction_scheme,model_components_dict,diff_regime)\n",
    "\n",
    "# build the model. Will save a file, don't include the date in the model filename\n",
    "model.build(date=False)\n",
    "\n",
    "# print out the parameters required for the model to run\n",
    "print(model.req_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up parameter evolution\n",
    "\n",
    "Below is where the model is finally constructed and run. I have set `alpha_s_0_2` to half the \"true\" value. We know that this value of `alpha_s_0_2`should be twice what it is now. Let's pretend we don't know that. We do know that `alpha_s_0_2` needs to be multiplied by some factor. This can be achieved by writing a parameter evolution function as outlined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the Simulate class\n",
    "from multilayerpy.simulate import Simulate\n",
    "\n",
    "# import the Parameter class\n",
    "from multilayerpy.build import Parameter\n",
    "\n",
    "# make the parameter dictionary\n",
    "# these are parameters from the KM-SUB model description paper (see reference at the start of this notebook)\n",
    "param_dict = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2':Parameter(4.2e-4 / 2.0),\n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(1e-5),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(1e-10),  # cm2 s-1\n",
    "              'Db_3':Parameter(1e-10),  # cm2 s-1\n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'H_2':Parameter(4.8e-4),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(7.0e13),  # cm-3\n",
    "              'Td_2': Parameter(1e-2,vary=True,bounds=(1e-3,1e-1)),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'k_1_2_surf':Parameter(6.0e-12)}  # cm2 s-1\n",
    "\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict)\n",
    "\n",
    "# let's load in the fake noisy data - columns are (time, y, y_error)\n",
    "fake_data = np.genfromtxt('noisy_data.txt')\n",
    "\n",
    "# set the Simulate.data attribute to be the fake data, using the Data object and normalised\n",
    "from multilayerpy.simulate import Data\n",
    "sim.data = Data(fake_data)\n",
    "\n",
    "#=======PARAMETER EVOLUTION FUNCTION=========\n",
    "\n",
    "# import copy module for deepcopying (see the function below)\n",
    "import copy\n",
    "\n",
    "# give it a sensible name\n",
    "# 1st argument MUST be time \n",
    "# 2nd argument MUST be the current state of the ODE integration at that time \n",
    "# (i.e. the number concentration (KM-SUB) or total number (KM-GAP) of each component in each layer - same format as y0 below)\n",
    "# 3rd argument MUST be the parameter dictionary\n",
    "# 4th *keyword* argument (optional) is a list of any additional varying parameters (useful to set the default as None)\n",
    "# These are the parameters that will vary\n",
    "\n",
    "def param_evo_func(t,y,param_dict,additional_params=None):\n",
    "    # we MUST make a deepcopy of the parameter dictionary\n",
    "    param_dict = copy.deepcopy(param_dict)\n",
    "    \n",
    "    # extract the coefficient we want to vary from the param list\n",
    "    # the list has only one member in this case, so we just extract that member\n",
    "    alpha_coeff = additional_params[0] \n",
    "    \n",
    "    # now let's multiply alpha_s_0_2 by this coefficient and set that value in the parameter dictionary\n",
    "    alpha = param_dict['alpha_s_0_2'].value\n",
    "    param_dict['alpha_s_0_2'].value = alpha * alpha_coeff\n",
    "    \n",
    "    # return the parameter dictionary\n",
    "    return param_dict\n",
    "\n",
    "# let's create the additional parameter, vary the value between 1 and 2.2. Initial guess of 1.\n",
    "alpha_coeff = Parameter(1.0,vary=True,bounds=(1.0,2.2),name='alpha_coeff')\n",
    "\n",
    "# now we need a list of additional varying Parameter objects\n",
    "add_params = [alpha_coeff]\n",
    "\n",
    "# now supply the parameter evolution function and the additional parameters to the Simulate object\n",
    "sim.param_evo_func = param_evo_func\n",
    "sim.param_evo_additional_params = add_params\n",
    "#============================================\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 10\n",
    "rp = 0.2e-4 # radius in cm\n",
    "time_span = [0,40] # in s\n",
    "n_time = 999 # number of timepoints to save to output\n",
    "\n",
    "#spherical V and A\n",
    "# use simulate.make_layers function\n",
    "V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "# initial conc. of everything\n",
    "\n",
    "bulk_conc_dict = {'1':1.21e21,'2':0,'3':0} # key=model component number, value=bulk conc, units: cm-3\n",
    "surf_conc_dict = {'1':9.68e13,'2':0,'3':0} # key=model component number, value=surf conc, units: cm-2\n",
    "\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "    \n",
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "sim.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot_bulk_concs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimize the model\n",
    "from multilayerpy.optimize import Optimizer\n",
    "\n",
    "fitter = Optimizer(sim)\n",
    "\n",
    "res = fitter.fit()\n",
    "\n",
    "# check sim.plot shows optimised result\n",
    "sim.plot()\n",
    "\n",
    "print('Optimised value of alpha_coeff: ', sim.param_evo_additional_params[0].value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excellent! The optimised value of `alpha_coeff` is ~2.0, which is what we know to be the \"true\" answer (we divided `alpha_s_0_2` by 2 at the beginning). \n",
    "\n",
    "### Summary\n",
    "This method of indtroducing additional varying parameters can prove useful when a user wants to include special parameterisations of model input parameters. \n",
    "\n",
    "A few points:\n",
    "* You can introduce time-dependent changes in model parameters (time is one of the arguments supplied to the function).\n",
    "* You can link a pair of model parameters together through one additional varying parameter. \n",
    "* This function can also be used as a sort of callback function. You can monitor how model outputs change over time. \n",
    "* Beware that adding more varying parameters to the model will increase the time needed to optimise it. "
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}