{ "cells": [ { "cell_type": "markdown", "id": "24e91247-ee5c-43c0-9116-c0cf9e87951b", "metadata": {}, "source": [ "# Introduction\n", "\n", "
\n", "\n", "Summary\n", "\n", "The model parameterisations refer to the set of parameters that are varied and estimated in a retrieval. Given that the computational time of retrievals largely depend on the number of fitted parameters, the input information in the atmospheric classes is often parameterised to reduce this number. ArchNEMESIS includes a large set of model parameterisations heritage from the original NEMESIS code and which might be useful for many applications. Nevertheless, new parameterisations can also be easily introduced to develop new methodologies for the analysis of planetary spectra.\n", "\n", "
\n", "\n", "The information about the model parameterisations is stored in the *Variables* class. Given the different forms parameterisations may take, this class needs to be versatile to store all the relevant information about the model parameters and the indications of what information must be changed in the forward model based on a particular parameterisation.\n", "\n", "The type of parameterisation is indicated by the VARIDENT parameter of the *Variables* class. For each parameterisation included in the model, VARIDENT is a three-element integer array. The logic behind this array is as follows:\n", "\n", "- Let's call the three elements of the array I$_1$, I$_2$, I$_3$.\n", "\n", "- I$_3$ indicates the ID for the particular parameterisation, following the identifiers listed in the next section.\n", "\n", "- If the selected model ID does not correspond to an atmospheric model, then I$_1$ and I$_2$ do not provide any useful information, and they can get arbitrary numbers. However, if the model ID corresponds to an atmospheric parameterisation, then I$_1$ and I$_2$ indicate the atmospheric parameter that must be modified. In particular:\n", "\n", " - If I$_1$ = 0, then the model parameterisation applies to the temperature profile.\n", " - If I$_1$ > 0, then the model parameterisation applies to the volume mixing ratio of a particular gas, with I$_1$ and I$_2$ representing the gas ID and isotope ID (e.g., I$_1$=1, I$_2$=0 for H$_2$O; I$_1$=1, I$_2$=4 for HDO).\n", " - If I$_1$ < 0, then we can have:\n", " - If 1 $\\leq$ |I$_1$| $\\leq$ NDUST, then the model parameterisation applies to the I$_1$ aerosol population.\n", " - If |I$_1$| = NDUST + 1, then the model parameterisation applies to the para-H$_2$ fraction profile.\n", " - If |I$_1$| = NDUST + 1, then the model parameterisation applies to the fractional cloud cover profile.\n", "\n", "\n", "
\n", "\n", "Examples of VARIDENT\n", "\n", "Let's assume we have an atmosphere with 3 aerosol populations. We include 5 model parameterisations, defined by:\n", "\n", "- I$_1$=2, I$_2$=0, I$_3$=0: Continuous vertical profile of CO$_2$ volume mixing ratio.\n", "- I$_1$=-2, I$_2$=0, I$_3$=2: Scaling factor for abundance of the second aerosol population in the atmosphere.\n", "- I$_1$=-4, I$_2$=0, I$_3$=0: Scaling factor for para-H2 fraction profile.\n", "- I$_1$=-1, I$_2$=0, I$_3$=32: Abundance of first aerosol population represented by base pressure, optical depth and fractional scale height.\n", "- I$_1$=0, I$_2$=0, I$_3$=231: Multiplication of spectrum by polynomial function\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "99af2a9a-5a7a-421b-93d8-2757322b1dd0", "metadata": {}, "source": [ "# Available models\n", "\n", "The model parameterisations have the purpose of modifying the information in the reference classes according to a smaller number of parameters. This is especially important for retrievals, when these parameters are modified in each iteration to find the best fit between the modelled and measured spectra. Each of the model parameterisations is defined by some Variable IDs and Variable parameters, which are defined in the *.apr* file. This file can be read using the *read_apr* function on the Variables class. \n", "\n", "You can run the command `python -c 'from archnemesis.Models import Models; print(Models.info())'` to display some summary information for each model.\n", "\n", "At this stage, only a number of models have been implemented in archNEMESIS with respect to NEMESIS. Detailed information about each of the model parameterisations is included in the Examples. In the following, a brief summary describing each of these models is included:\n" ] }, { "cell_type": "raw", "id": "78529677-a6c8-462d-9a3d-8730c71661f0", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 0: Continuous atmospheric vertical profile.\n", " \n", " \n", "
\n", "
" ] }, { "cell_type": "raw", "id": "6d0ff96d-c477-4bb5-b403-a80c7e01d528", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 2: Scaling of an atmospheric vertical profile.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "338d9d9f-0921-401d-ac99-9fe3c49b7e40", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 3: Scaling of an atmospheric vertical profile (carried in log-space).\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "a13bf768-245b-4445-94d9-436f53a7528e", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 9: Cloud represented by base height, optical depth and fractional scale height.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "fcc9a0f5-1548-4020-b44c-f155c9a70fe1", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 32: Cloud represented by base pressure, optical depth and fractional scale height.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "2b3b3945-ad78-4cf0-b5dd-82ffee076296", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 43: Temperature profile from double grey analytic formulation from Parmentier and Guillot (2014).\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "47410b28", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 45: Irwin Ice Giant CH4 model.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "d000d99c", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 47: Gaussian cloud profile.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "fe790fce", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 51: Scaling of another profile.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "4931f3a4-d11e-4810-852a-7e19f567b61b", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 110: Model for Venus cloud from Haus et al. (2016) with an altitude offset.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "f075942a-6cea-4dfa-a5bb-87d9b1b0d01b", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 111: Model for Venus cloud from Haus et al. (2016) with an altitude offset and coupled SO2 volume mixing ratio profile.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "460dbced-7af5-496b-820f-495d122ea704", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 229: Model for double-Gaussian instrument line shape for TGO/ACS (Alday et al. 2019).\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "c6320cb9-1ee7-4b45-b3b8-a614b8fa5b37", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 231: Multiplication of spectrum by polynomial function.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "14491917", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 444: Size distribution and imaginary refractive index spectrum of aerosol species.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "15639de1-411d-498d-9776-3fcd747c966a", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 447: Specification of Doppler velocity.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "003a16e4-a9de-49c5-89f3-4f315d3f7582", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 666: Retrieval of pressure at specified height.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "68a3b629-d7b8-4471-bfa6-8401b79977cc", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 667: Scaling of spectrum by a dilution factor.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "raw", "id": "a37a607d-1d42-4876-86f4-afbad9bac37b", "metadata": {}, "source": [ "\n", "\n", "
\n", "
\n", " Model 999: Retrieval of surface temperature.\n", "\n", " \n", "\n", "
\n", "
" ] }, { "cell_type": "markdown", "id": "2a8e1c2c-c292-46de-a401-1432b77509d4", "metadata": {}, "source": [ "| Model ID | No. of parameters | No. of extra parameters | Description |\n", "| --- | --- | --- | --- |\n", "| 49 | NP | 0 | Continuous atmospheric vertical profile in linear scale (i.e., allowing negative abundances) |\n", "| 50 | NP | 0 | Continuous profile of scaling factors in linear scale |\n", "| 2310 | NGEOM $\\times$ (NDEGREE+1) $\\times$ NWINDOWS | 3 | Same as model 231 but applying it to different spectral windows |" ] }, { "cell_type": "markdown", "id": "c3e37799-735a-4dae-ab11-b0a071de1627", "metadata": {}, "source": [ "# Custom Model Class\n", "\n", "Users may want to define their own custom parameterisations ande create new retrieval schemes for their specific purposes. Implementing new parameterisations in archNEMESIS is relatively straightforward, at the simplest level you just need to add a new class to one of the *Models/PreRTModels.py*, or *Models/PostRTModels.py* files. \n", "\n", "We are in the process of rationalising the model classes, so we recommend using the below model classes as a guide:\n", "* In the *Models/PreRTModels.py* file:\n", " - *Model2*\n", " - *Model3*\n", " - *Model9*\n", " - *Model32*\n", " - *Model45*\n", " - *Model47*\n", " - *Model444*\n", "* In the *Models/PostRTModels.py* file:\n", " - *Model231*\n", "\n", "NOTE: In each file there is a template model class, that is the result of the steps (2) and (3) (detailed below), that has the methods that need to be implemented. The body of each method body consists of `raise NotImplementedError('This is a template model and should never be used')` as the first expression and some example code with explanatory comments, just copy the template model class, rename it, give it a unique ID, and replace the method bodies with the required code for your model. \n", "\n", "The templates model classes are:\n", "* *TemplatePostRTModel* in the *Models/PreRTModels.py* file\n", "* *TemplatePreRTModel* in the *Models/PostRTModels.py* file\n", "\n", "\n", "## General Steps\n", "\n", "1. Work out what kind of model you are creating, this informs where you put your new model class.\n", " - Models applied before the radiative transfer calculation are in *Models/PreRTModels.py*, there are generally two types:\n", " - Atmospheric models affect one of the following atmosphere profiles: \n", " * gas volume mixing ratio\n", " * aerosol species (dust) density\n", " * temperature\n", " * para-h2 fraction\n", " * fractional cloud cover\n", " - Non-Atmospheric models affect other things, e.g.:\n", " * Optical properties of an aerosol species\n", " * Instrument lineshape\n", " * Doppler shift of the observer w.r.t the target\n", " - Models applied after the radiative transfer calculation are in *Models/PostRTModels.py*, currently there is only one type:\n", " - Spectral models affect the the modelled output spectrum directly.\n", "\n", "2. Look at the *ModelBase* class in the *Models/ModelBase.py* file. List any **abstract methods** on the class (they are marked with `@abstractmethod`)\n", "\n", "3. In the file that corresponds to the type of model you are creating, look at the \"model type base class\" (at the top of the file) : *PreRTModelBase* for *Models/PreRTModels.py*; *PostRTModelBase* for *Models/PostRTModels.py*. Do the following:\n", " * Any methods with the same name as an **abstract method** of the *ModelBase* class does **NOT** have to be implemented, as the \"model type base class\" already provides an implementation for you. You can cross these off your list.\n", " \n", " * Any **abstract methods** on the \"model type base class\" (they are marked with `@abstractmethod`), **DO** have to be implemented. Add these to your list.\n", "\n", "4. Write a new class that implements all of the **abstract methods** on your list, and give it an integer class attribute `id` that is a unique ID number that will be used to identify your custom model, there is a test to ensure all models have a unique ID number so run the tests after adding a model to ensure everything is set up properly. Use the other implementations in the file as a guide.\n", "\n", "\n", "## Example Custom Model Class\n", "\n", "As it is much easier to see the above process in action, we present a walkthrough example of implementing a custom model. The example we are using is not going to be *useful*, but it serves as a guide that shows how to do something non-trivial.\n", "\n", "Our custom model will parameterise the volume mixing ratio of a gas as a stepwise function. This is not really that *useful*, but it gives us the opportunity to explore how to do a few different things.\n", "\n", "Following the above procedure:\n", "\n", "* Step 1 is to work out what kind of model we are creating. In this case the model is an Atmospheric model (as it alters a gas volume mixing ratio), so we will put the model in the *Models/PreRTModels.py* file.\n", "\n", "* Step 2 is to look at *ModelBase* and *PreRTModelBase* and write down the **abstract methods** that require an implementation. At the time of writng (30/06/2025), these are the relevant methods on the two classes:\n", " \n", " - *ModelBase* **abstract methods**:\n", " + `is_varident_valid`\n", " + `from_apr_to_state_vector`\n", " + `calculate`\n", " + `calculate_from_subprofretg`\n", " + `patch_from_subprofretg`\n", " + `calculate_from_subspecret`\n", " \n", " - *PreRTModelBase* **concrete methods** (i.e. non-abstract):\n", " + `is_varident_valid`\n", " + `patch_from_subprofretg`\n", " + `calculate_from_subspecret`\n", " \n", " - *PreRTModelBase* **abstract methods**:\n", " + `calculate_from_subprofretg`\n", "\n", "Therefore, if we combine the **abstract methods** from both classes, and then remove the **concrete methods** of *PreRTModelBase*, we have the following list of methods that our new class needs to implement: `from_apr_to_state_vector`, `calculate`, `calculate_from_subprofretg`. This is step 3.\n", "\n", "We can use the template class in *Models/PreRTModels.py* to speed up steps (2) and (3), we just copy the template class, rename it, and re-write its methods to do what we need.\n", "\n", "Step 4 is to actually implement the methods, and is performed in the following sections.\n", "\n", "\n", "### first definitions of the model class\n", "\n", "First off, we need a name. The current model names are things like *Model45*, they are holdovers from the FORTRAN code and are not very descriptive. Ideally we want something that describes what the model represents, \"PiecewiseGasVMR\" is a decent choice so we will use that. \n", "\n", "Secondly, we need an ID number, we should choose something that is somewhat consistent with the current ID numbers if possible. I will use 10 for this example.\n", "\n", "Thirdly, we should add a docstring to the class that describes it.\n", "\n", "At this point our model class looks like (some code omitted for clarity):\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " \"\"\"\n", " Parameterises the gas volume mixing ratio as a piecewise function with `n_chunk` chunks.\n", " The chunk values are retrieved.\n", " \"\"\"\n", " id : int = 10\n", " \n", " def __init__(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def calculate(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " # OMITTED CODE #\n", " \n", " def calculate_from_subprofretg(\n", " # OMITTED CODE #\n", "```\n", "\n", "### defining the model scope\n", "\n", "The model class will need to know the following pieces of information\n", "* The number of chunks of the piecewise function\n", "* the gas ID we want to apply the profile to\n", "* the apriori values of the chunks\n", "\n", "I find it easiest to work in the methods in the order they are called by the program,\n", "the methods are called in the order:\n", "* `from_apr_to_state_vector` - is called first from `Variables_0::read_apr` when the .apr file is read.\n", "* `__init__` - is called inside `from_apr_to_state_vector` as that method needs to construct the model class instance.\n", "* `calculate_from_subprofretg` is called when forward modelling, inside `ForwardModel_0::subprofretg`.\n", "* `calculate` is called inside `calculate_from_subprofretg` as that method performs the actual calculation.\n", "\n", "\n", "### writing `from_apr_to_state_vector`\n", "\n", "Reading the docstring of the `from_apr_to_state_vector` method shows what arguments are passed to us from the `Variables_0::read_apr` method. \n", "\n", "#### unpacking values from VARIDENT\n", "\n", "Considering the `varident` for the model. Reading the documentation, `PreRTModelBase::is_varident_valid` method, and other model class's `from_apr_to_state_vector` implementations shows that, when we are dealing with an atmospheric model, the three integers in `varident` are defined as:\n", "\n", "* `varident[0]` \n", " - `== 0` denotes a temperature profile\n", " - `> 0` denotes a gas ID that is equal to `varident[0]`.\n", " - `< 0`\n", " + is an aerosol species profile when the magnitude is less than or equal to the number of aerosol species\n", " + a para-h2 fraction profile when the magnitude is equal to the number of aerosol species + 1\n", " + a fractional cloud cover profile when the magnitude is equal to the number of aerosol species + 2\n", "* `varident[1]` is the ID of the gas *isotope* **if** `varident[0]` is a gas ID.\n", "* `model_id` denotes the ID number of the model to be used.\n", "\n", "Therefore,\n", "1) we can assume that `varident[2]` is equal to the `cls.id` of our model. This is what `PreRTModelBase::is_varident_valid` ensures.\n", "2) our model requires `varident[0]` to be positive as it must be a gas ID\n", "3) `varident[1]` is an isotope ID.\n", "\n", "We can add the following code to the `from_apr_to_state_vector` method:\n", "```python\n", "# unpack values from VARIDENT\n", "gas_id = varident[0]\n", "assert gas_id > 0, f'model {cls.__name__} with id={cls.id} must have a gas ID as the first VARIDENT value'\n", "iso_id = varident[1]\n", "```\n", "\n", "#### reading apriori parameters\n", "\n", "We want to read the apriori coefficients from the .apr file. The simplest way to do this is to add some extra lines to the .apr file after the VARIDENT for our model. The first line will be the number of chunks, the next `n_chunk` lines be a space separated list of the chunk values and their errors.\n", "\n", "The `f` argument is the file descriptor of the open .apr file at the line *after* the VARIDENT for our model. Therefore we just need to read `n_chunk+1` lines from the file and unpack the values in the way we just said.\n", "\n", "The `npro` argument is the number of points in each profile, we want our piecewise function to have at least two points shared in each chunk, so we want to make sure we have enough points in the profile.\n", "\n", "We can add the following code to the `from_apr_to_state_vector` method:\n", "```python\n", "# read in apriori parameters\n", "# read number of chunks\n", "n_chunk = int(f.readline().strip())\n", "\n", "# ensure we have at least 2 profile points per chunk\n", "assert npro >= 2*n_chunk, f'{cls.__name__} id={cls.id} must have at most half as many chunks as there are points in a profile, have {n_chunk} vs {npro}'\n", "\n", "# create holder arrays\n", "chunk = np.zeros((n_chunk,),float)\n", "chunk_err = np.zeros((n_chunk,),float)\n", "\n", "# read value and error of each chunk\n", "for i in range(coeff.size):\n", " chunk[i], chunk_err[i] = tuple(map(float, f.readline().strip().split()))\n", "```\n", "\n", "#### packing the state vector and the covariance matrix\n", "\n", "We need to pack the apriori parameter values into the statevector and pack the parameter errors into the covariance vector. Note, as we have not set any correlation between our chunks, we can just set the diagonal covariance matrix elements.\n", "\n", "The `ix` argument is the next unused index of the state vector. We will need `ix` along with:\n", "* `lx` - the log flag for the state vector, denotes if a parameter is stored as a log of its \"real\" value.\n", "* `x0` - the actual state vector\n", "* `sx` - the covariance matrix\n", "* `inum` - flag for the state vector, denotes if the derivative of a parameter should be numerically calculated or not.\n", "\n", "Therefore, we want to:\n", "1) store the parameters in the state vector. NOTE: we should store the log of their value as the chunk values cannot be negative.\n", "2) store the errors in the covariance matrix.\n", "3) set the log flag for our parameters\n", "4) for simplicity, set to numerically calculate the derivatives of our parameters\n", "\n", "We can add the following code to the `from_apr_to_state_vector` method:\n", "```python\n", "# Packing the state vector and covariance matrix\n", "ix_0 = ix # store the initial value\n", "for i in range(n_chunk):\n", " # store the log(chunk) values\n", " lx[ix] = 1\n", " inum[ix] = 1\n", " x0[ix] = np.log(chunk[i])\n", " \n", " # store the chunk variance, must use relative errors as we are storing log(value)\n", " sx[ix,ix] = (chunk_err[i]/chunk[i])**2\n", " \n", " ix += 1\n", "```\n", "\n", "#### construct the model class instance and return\n", "\n", "Finally, we need to construct the model instance. To enable us to easily select the model when (for example) creating a plot, it is a good idea to classify the model this is generally only done for \"Atmospheric\" models (ones that influence the `Atmosphere` attribute of the forward model). We should also pass any constants we will need later to the constructor so we can store them on the instance.\n", "\n", "We always need to pass the index of the state vector where the parameters of the model start, and the number of parameters the model stores in the state vector. Then we can pass any constants. As they might be useful later, we will also pass `n_chunk`, `gas_id`, and `iso_id` to be stored on the model.\n", "\n", "We can add the following code to the `from_apr_to_state_vector` method:\n", "```python\n", "\n", "# Get the model classification\n", "model_classification = variables.classify_model_type_from_varident(varident, ngas, ndust)\n", "assert issubclass(cls, model_classification[0]), \"Model base class must agree with the classification from Variables_0::classify_model_type_from_varident\"\n", "\n", "# return the constructed model instance\n", "return cls(ix_0, ix-ix_0, model_classification[1], n_chunk, gas_id, iso_id)\n", "```\n", "\n", "\n", "#### code listing for `from_apr_to_state_vector`\n", "\n", "The altered code is below, I have omitted docstrings and comments etc. for clarity.\n", "\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " # CODE UNCHANGED FROM EARLIER #\n", " \n", " def __init__(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def calculate(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " cls,\n", " variables : \"Variables_0\",\n", " f : IO,\n", " varident : np.ndarray[[3],int],\n", " varparam : np.ndarray[[\"mparam\"],float],\n", " ix : int,\n", " lx : np.ndarray[[\"mx\"],int],\n", " x0 : np.ndarray[[\"mx\"],float],\n", " sx : np.ndarray[[\"mx\",\"mx\"],float],\n", " inum : np.ndarray[[\"mx\"],int],\n", " npro : int,\n", " ngas : int,\n", " ndust : int,\n", " nlocations : int,\n", " runname : str,\n", " sxminfac : float,\n", " ) -> Self:\n", " \n", " # unpack values from VARIDENT\n", " gas_id = varident[0]\n", " assert gas_id > 0, f'model {cls.__name__} with id={cls.id} must have a gas ID as the first VARIDENT value'\n", " iso_id = varident[1]\n", " \n", " # read in apriori parameters\n", " # read number of chunks\n", " n_chunk = int(f.readline().strip())\n", " \n", " # ensure we have at least 2 profile points per chunk\n", " assert npro >= 2*n_chunk, f'{cls.__name__} id={cls.id} must have at most half as many chunks as there are points in a profile, have {n_chunk} vs {npro}'\n", " \n", " # create holder arrays\n", " chunk = np.zeros((n_chunk,),float)\n", " chunk_err = np.zeros((n_chunk,),float)\n", "\n", " # read value and error of each chunk\n", " for i in range(coeff.size):\n", " chunk[i], chunk_err[i] = tuple(map(float, f.readline().strip().split()))\n", " \n", " # Packing the state vector and covariance matrix\n", " ix_0 = ix # store the initial value\n", " for i in range(n_chunk):\n", " # store the log(chunk) values\n", " lx[ix] = 1\n", " inum[ix] = 1\n", " x0[ix] = np.log(chunk[i])\n", " \n", " # store the chunk variance, must use relative errors as we are storing log(value)\n", " sx[ix,ix] = (chunk_err[i]/chunk[i])**2\n", " \n", " ix += 1\n", " \n", " # Classify the model\n", " model_classification = variables.classify_model_type_from_varident(varident, ngas, ndust)\n", " assert issubclass(cls, model_classification[0]), \"Model base class must agree with the classification from Variables_0::classify_model_type_from_varident\"\n", "\n", " \n", " # return the constructed model instance\n", " return cls(ix_0, ix-ix_0, model_classification[1], n_chunk, gas_id, iso_id)\n", " \n", " def calculate_from_subprofretg(\n", " # OMITTED CODE #\n", "```\n", "\n", "\n", "### writing `__init__(...)`\n", "\n", "This method calls the parent class constructor, stores any data we want on the model instance, and defines the parameters of the model to make it easy to pull the model's information from the state vector and display it.\n", "\n", "\n", "#### calling the parent class constructor\n", "\n", "The template gives us the code for this, we just need to a copy it. NOTE: If we were not making an \"Atmospheric\" model, we would probably not bother forwarding `atm_profile_type` to the parent class, however it should not be a problem to do so.\n", "\n", "Add the following code to the `__init__` method:\n", "```python\n", "# initialise the parent class\n", "super().__init__(i_state_vector_start, n_state_vector_entries, atm_profile_type)\n", "```\n", "\n", "#### store data on the model instance\n", "\n", "Earlier we sent `n_chunk`, `gas_id`, and `iso_id` to the class constructor. We should make sure we accept these arguments, then store their values on the class instance.\n", "\n", "Alter the definition of the `__init__` method to read:\n", "```python\n", "def __init__(\n", " self, \n", " i_state_vector_start : int, \n", " n_state_vector_entries : int,\n", " atm_profile_type : AtmosphericProfileType,\n", " n_chunk : int,\n", " gas_id : int,\n", " iso_id : int,\n", "):\n", "```\n", "\n", "Add the following code to the `__init__` method:\n", "```python\n", "# storing constants\n", "self.n_chunk = n_chunk\n", "self.gas_id = gas_id\n", "self.iso_id = iso_id\n", "```\n", "\n", "#### define model parameters\n", "\n", "Finally, we define the model parameters.\n", "\n", "The piecewise model is fairly simple, we only actually have a single parameter `chunk`, but is is a multi-value parameter (i.e. it is an array not a single number). The template model class gives us an example of how to define this kind of parameter. As the `chunk` parameter takes up all of the space in the sub-section of the state-vector assigned to the model, we can just use `slice(None)` to select the entire sub-section of the state vector. We should also give a description, and the units of the parameter.\n", "\n", "Add the following code to the `__init__` method:\n", "```python\n", "# define the model parameters, specifically how they are stored in the state vector\n", "self.parameters : tuple[ModelParameter,...] = (\n", " ModelParameter(\n", " 'chunk',\n", " slice(None),\n", " f'the {n_chunk} values of a piecewise function that define a volume mixing ratio for gas : {gas_id} isotope : {iso_id}',\n", " 'RATIO'\n", " ), # NOTE: the comma here is required.\n", ")\n", "```\n", "\n", "#### code listing for `__init__`\n", "\n", "The altered code is below, I have omitted docstrings and comments etc. for clarity.\n", "\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " # CODE UNCHANGED FROM EARLIER #\n", " \n", " def __init__(\n", " self, \n", " i_state_vector_start : int, \n", " n_state_vector_entries : int,\n", " atm_profile_type : AtmosphericProfileType,\n", " n_chunk : int,\n", " gas_id : int,\n", " iso_id : int,\n", " ):\n", " # initialise the parent class\n", " super().__init__(i_state_vector_start, n_state_vector_entries, atm_profile_type)\n", " \n", " # storing constants\n", " self.n_chunk = n_chunk\n", " self.gas_id = gas_id\n", " self.iso_id = iso_id\n", " \n", " # define the model parameters, specifically how they are stored in the state vector\n", " self.parameters : tuple[ModelParameter,...] = (\n", " ModelParameter(\n", " 'chunk',\n", " slice(None),\n", " f'the {n_chunk} values of a piecewise function that define a volume mixing ratio for gas : {gas_id} isotope : {iso_id}',\n", " 'RATIO'\n", " ), # NOTE: the comma here is required.\n", " )\n", " \n", " \n", " @classmethod\n", " def calculate(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " # OMITTED CODE #\n", " \n", " def calculate_from_subprofretg(\n", " # OMITTED CODE #\n", "```\n", "\n", "### writing `calculate_from_subprofretg`\n", "\n", "This method is called from `ForwardModel_0::subprofretg` and is how an atmospheric model influences the state of the retrieval. We have to perform a few tasks in this method.\n", "\n", "1) Unpack parameters from the state vector and get other values we will need\n", "2) Use the `self.calculate(...)` method (see later) to calculate the gas volume mixing ratio according to our model\n", "3) Put the results of the calculation in their correct palce.\n", "\n", "The `forward_model` argument is the current `ForwardModel_0` instance, we can access all of its components using this argument. It is very useful as in principle a model class may alter practically anything about a retrieval by interacting with this argument.\n", "\n", "\n", "#### Unpack parameters and values\n", "\n", "In our case we only have one parameter, and there is a helpful method `self.get_parameter_values_from_state_vector` (defined on the `ModelBase` class that we have indirectly inherited from) that we can use to get the values of all parameters simultaniously. We need to pass the state vector and the state vector log flags to this method. Those are stored on the `Variables` component of the `forward_model` as `forward_model.Variables.XN` and `forward_model.Variables.LX` respectively.\n", "\n", "We can also find the index of the gas volume mixing ratio we need to alter. There are actually two ways to do this. The first uses the `ipar` argument, which encodes the profile type and profile index that an atmospheric profile should alter. We can get the profile type and index from `ipar` using the `Atmosphere_0::ipar_to_atm_profile_type(...)` method on the `forward_model.AtmosphereX` object. We can also get it by using the `Atmosphere_0::locate_gas(...)` method on the same object as we have stored the `gas_id` and `iso_id` values on the model instance. In this case we should do both and check they agree. \n", "\n", "We can also make sure we have been given the correct profile type. Checks like this may seem superfulous, but they make sure that if something does go wrong the error is caught as early as possible, which makes debugging much easier. The profile types are defined in the `enums.py` file.\n", "\n", "Add the following code to the `calculate_from_subprofretg` method:\n", "```python\n", "# unpack parameters from state vector\n", "chunk = self.get_parameter_values_from_state_vector(forward_model.Variables.XN, forward_model.Variables.LX)\n", "\n", "# get the profile type and index, and check the values\n", "atm = forward_model.AtmosphereX\n", "atm_profile_type, atm_profile_idx = atm.ipar_to_atm_profile_type(ipar)\n", "\n", "assert atm_profile_type == AtmosphericProfileType.GAS_VOLUME_MIXING_RATIO, f'model {self.__class__.__name__} with id={self.id}. Only accepts gas volume mixing ratio profiles'\n", "\n", "assert atm_profile_idx == atm.locate_gas(self.gas_id, self.iso_id), f'model {self.__class__.__name__} with id={self.id}. The two different ways of getting the gas_vmr_idx must agree'\n", "```\n", "\n", "\n", "#### send to `self.calculate(...)`\n", "\n", "We have a lot of freedom to define how the class method `self.calculate` works. However, the conventions used by the other models are a good place to start. We can always alter the method later if we need to. For now we will send the atmosphere class, the profiel type and index, and the model parameters. We will assume that we will get back an updated atmosphere class and an array of functional derivatives.\n", "\n", "Add the following code to the `calculate_from_subprofretg` method:\n", "```python\n", "atm, xmap1 = self.calculate(\n", " atm,\n", " atm_profile_type,\n", " atm_profile_idx,\n", " chunk\n", ")\n", "```\n", "\n", "\n", "#### store the results\n", "\n", "We need to put the return `Atmosphere_0` class into the same bit of the `forward_model` argument we got it from. And we need to update the `xmap` argument with the calculated values in `xmap1`.\n", "\n", "`xmap` is a numpy array of shape `(nx, NVMR+2+NDUST, NP, NLOCATIONS)` and holds the functional derivatives of the `nx` values of the state vector, w.r.t the `NVMR+2+NDUST` profiles, `NP` profile points, and `NLOCATIONS` locations. As we are only considering one location right now, and we know we will only be altering one profile (in this case indexed by `ipar`) we can set the correct section of `xmap` from `xmap1`.\n", "\n", "Add the following code to the `calculate_from_subprofretg` method:\n", "```python\n", "forward_model.AtmosphereX = atm\n", "xmap[self.state_vector_slice, ipar, 0:atm.NP] = xmap1\n", "```\n", "\n", "\n", "#### code listing for `calculate_from_subprofretg`\n", "\n", "The altered code is below, I have omitted docstrings and comments etc. for clarity.\n", "\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " # CODE UNCHANGED FROM EARLIER #\n", " \n", " def __init__(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def calculate(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " # OMITTED CODE #\n", " \n", " def calculate_from_subprofretg(\n", " self,\n", " forward_model : \"ForwardModel_0\",\n", " ix : int,\n", " ipar : int,\n", " ivar : int,\n", " xmap : np.ndarray,\n", " ) -> None:\n", " \n", " # unpack parameters from state vector\n", " chunk = self.get_parameter_values_from_state_vector(forward_model.Variables.XN, forward_model.Variables.LX)\n", "\n", " # get the profile type and index, and check the values\n", " atm = forward_model.AtmosphereX\n", " atm_profile_type, atm_profile_idx = atm.ipar_to_atm_profile_type(ipar)\n", "\n", " assert atm_profile_type == AtmosphericProfileType.GAS_VOLUME_MIXING_RATIO, f'model {self.__class__.__name__} with id={self.id}. Only accepts gas volume mixing ratio profiles'\n", "\n", " assert atm_profile_idx == atm.locate_gas(self.gas_id, self.iso_id), f'model {self.__class__.__name__} with id={self.id}. The two different ways of getting the gas_vmr_idx must agree'\n", " \n", " atm, xmap1 = self.calculate(\n", " atm,\n", " atm_profile_type,\n", " atm_profile_idx,\n", " chunk\n", " )\n", " \n", " forward_model.AtmosphereX = atm\n", " xmap[self.state_vector_slice, ipar, 0:atm.NP] = xmap1\n", " \n", "```\n", "\n", "\n", "### writing `calculate`\n", "\n", "Finally we have to perform the actual calculation. All we have to do is work out the new profile values, and the functional derivatives. The interface for the `self.calcualte(...)` method has been effectively set for us in the previous section. Therefore we can go right ahead and define it.\n", "\n", "Alter the `calculate` class method definition to be:\n", "```python\n", "@classmethod\n", "def calculate(\n", " cls,\n", " atm : \"Atmosphere_0\",\n", " atm_profile_type : AtmosphericProfileType,\n", " atm_profile_idx : int | None,\n", " chunk : np.ndarray[[int],float],\n", " ):\n", "```\n", "\n", "#### find new profile values\n", "\n", "Our piecewise function sets the profile values in chunks. If we have 2 chunks for example, we set one half of the profile to the first chunk's value, and the second half to the second chunk's value.\n", "\n", "Add the following code to the `calculate` method:\n", "```python\n", "# update profile values\n", "temp = np.array(atm.VMR)\n", "fidx = np.linspace(0,chunk.size,temp.size)\n", "j=0\n", "for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " temp[j, atm_profile_idx] = chunk[i]\n", " j += 1\n", "atm.edit_VMR(temp)\n", "```\n", "\n", "#### find functional derivatives\n", "\n", "For functional derivatives `xmap` of shape `(nx, NVMR+2+NDUST, NP, NLOCATIONS)` assume the state vector values are labelled by 'K', and the profile points are labelled by 'J'. Then `xmap[K, :, J, :]` is the rate of change in the 'J^th' profile point w.r.t the 'K^th' state vector entry for each profile and location. We are restricting outselves to 1 location, and we are operating on 1 specific gas profile, so `xmap[K,ipar,J] = d(P[J])/d(X[K])` where `P[J]` is the J^th profile point and `X[K]` is the K^th state vector entry. \n", "\n", "Therefore as our piecewise function just sets profile values to the same value as the parameter, `d(P[J])/d(X[K]) = 1` when `J` influences `K` and zero otherwise.\n", "\n", "After we have found `xmap` we can return the results\n", "\n", "Add the following code to the `calculate` method:\n", "```python\n", "# find the functional derivatives\n", "xmap = np.zeros((chunk.size,atm.VMR.size[0]),float)\n", "j=0\n", "for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " xmap[i, j] = 1\n", " j += 1\n", "\n", "# return results\n", "return atm, xmap\n", "```\n", "\n", "\n", "#### code listing for `calculate`\n", "\n", "The altered code is below, I have omitted docstrings and comments etc. for clarity.\n", "\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " # CODE UNCHANGED FROM EARLIER #\n", " \n", " def __init__(\n", " # OMITTED CODE #\n", " \n", " @classmethod\n", " def calculate(\n", " cls,\n", " atm : \"Atmosphere_0\",\n", " atm_profile_type : AtmosphericProfileType,\n", " atm_profile_idx : int | None,\n", " chunk : np.ndarray[['mparam'],float],\n", " ):\n", " # update profile values\n", " temp = np.array(atm.VMR)\n", " fidx = np.linspace(0,chunk.size,temp.size)\n", " j=0\n", " for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " temp[j, atm_profile_idx] = chunk[i]\n", " j += 1\n", " atm.edit_VMR(temp)\n", " \n", " # find the functional derivatives\n", " xmap = np.zeros((chunk.size,atm.VMR.size[0]),float)\n", " j=0\n", " for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " xmap[i, j] = 1\n", " j += 1\n", "\n", " # return results\n", " return atm, xmap\n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " # OMITTED CODE #\n", " \n", " def calculate_from_subprofretg(\n", " # OMITTED CODE #\n", " \n", "```\n", "\n", "\n", "### Full final code listing\n", "\n", "Below is the full code listing for this example model class. You can add it to the *Models/PreRTModels.py* file and run `python -c 'from archnemesis.Models import Models; print(Models.info(10))'` to see a summary of this model. \n", "\n", "NOTE: This is an example model created as an learning aid. It has not been tested so the exact implementation may have some bugs, but the process of creation is correct. Whenever creating a new model, think about adding a test to check the new model performs as expected.\n", "\n", "```python\n", "class PiecewiseGasVMR(PreRTModelBase):\n", " \"\"\"\n", " Parameterises the gas volume mixing ratio as a piecewise function with `n_chunk` chunks.\n", " The chunk values are retrieved.\n", " \"\"\"\n", " id : int = 10\n", " \n", " def __init__(\n", " self, \n", " i_state_vector_start : int, \n", " n_state_vector_entries : int,\n", " n_chunk : int,\n", " gas_id : int,\n", " iso_id : int,\n", " ):\n", " # initialise the parent class\n", " super().__init__(i_state_vector_start, n_state_vector_entries)\n", " \n", " # storing constants\n", " self.n_chunk = n_chunk\n", " self.gas_id = gas_id\n", " self.iso_id = iso_id\n", " \n", " # define the model parameters, specifically how they are stored in the state vector\n", " self.parameters : tuple[ModelParameter,...] = (\n", " ModelParameter(\n", " 'chunk',\n", " slice(None),\n", " f'the {n_chunk} values of a piecewise function that define a volume mixing ratio for gas : {gas_id} isotope : {iso_id}',\n", " 'RATIO'\n", " ), # NOTE: the comma here is required.\n", " )\n", " \n", " \n", " @classmethod\n", " def calculate(\n", " cls,\n", " atm : \"Atmosphere_0\",\n", " atm_profile_type : AtmosphericProfileType,\n", " atm_profile_idx : int | None,\n", " chunk : np.ndarray[['mparam'],float],\n", " ):\n", " # update profile values\n", " temp = np.array(atm.VMR)\n", " fidx = np.linspace(0,chunk.size,temp.size)\n", " j=0\n", " for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " temp[j, atm_profile_idx] = chunk[i]\n", " j += 1\n", " atm.edit_VMR(temp)\n", " \n", " # find the functional derivatives\n", " xmap = np.zeros((chunk.size,atm.VMR.size[0]),float)\n", " j=0\n", " for i in range(chunk.size):\n", " while fidx[j] < i+1:\n", " xmap[i, j] = 1\n", " j += 1\n", "\n", " # return results\n", " return atm, xmap\n", " \n", " \n", " @classmethod\n", " def from_apr_to_state_vector(\n", " cls,\n", " variables : \"Variables_0\",\n", " f : IO,\n", " varident : np.ndarray[[3],int],\n", " varparam : np.ndarray[[\"mparam\"],float],\n", " ix : int,\n", " lx : np.ndarray[[\"mx\"],int],\n", " x0 : np.ndarray[[\"mx\"],float],\n", " sx : np.ndarray[[\"mx\",\"mx\"],float],\n", " inum : np.ndarray[[\"mx\"],int],\n", " npro : int,\n", " nlocations : int,\n", " runname : str,\n", " sxminfac : float,\n", " ) -> Self:\n", " \n", " # unpack values from VARIDENT\n", " gas_id = varident[0]\n", " assert gas_id > 0, f'model {cls.__name__} with id={cls.id} must have a gas ID as the first VARIDENT value'\n", " iso_id = varident[1]\n", " \n", " # read in apriori parameters\n", " # read number of chunks\n", " n_chunk = int(f.readline().strip())\n", " \n", " # ensure we have at least 2 profile points per chunk\n", " assert npro >= 2*n_chunk, f'{cls.__name__} id={cls.id} must have at most half as many chunks as there are points in a profile, have {n_chunk} vs {npro}'\n", " \n", " # create holder arrays\n", " chunk = np.zeros((n_chunk,),float)\n", " chunk_err = np.zeros((n_chunk,),float)\n", "\n", " # read value and error of each chunk\n", " for i in range(coeff.size):\n", " chunk[i], chunk_err[i] = tuple(map(float, f.readline().strip().split()))\n", " \n", " # Packing the state vector and covariance matrix\n", " ix_0 = ix # store the initial value\n", " for i in range(n_chunk):\n", " # store the log(chunk) values\n", " lx[ix] = 1\n", " inum[ix] = 1\n", " x0[ix] = np.log(chunk[i])\n", " \n", " # store the chunk variance, must use relative errors as we are storing log(value)\n", " sx[ix,ix] = (chunk_err[i]/chunk[i])**2\n", " \n", " ix += 1\n", " \n", " # return the constructed model instance\n", " return cls(ix_0, ix-ix_0, n_chunk, gas_id, iso_id)\n", " \n", " \n", " def calculate_from_subprofretg(\n", " self,\n", " forward_model : \"ForwardModel_0\",\n", " ix : int,\n", " ipar : int,\n", " ivar : int,\n", " xmap : np.ndarray,\n", " ) -> None:\n", " \n", " # unpack parameters from state vector\n", " chunk = self.get_parameter_values_from_state_vector(forward_model.Variables.XN, forward_model.Variables.LX)\n", "\n", " # get the profile type and index, and check the values\n", " atm = forward_model.AtmosphereX\n", " atm_profile_type, atm_profile_idx = atm.ipar_to_atm_profile_type(ipar)\n", "\n", " assert atm_profile_type == AtmosphericProfileType.GAS_VOLUME_MIXING_RATIO, f'model {self.__class__.__name__} with id={self.id}. Only accepts gas volume mixing ratio profiles'\n", "\n", " assert atm_profile_idx == atm.locate_gas(self.gas_id, self.iso_id), f'model {self.__class__.__name__} with id={self.id}. The two different ways of getting the gas_vmr_idx must agree'\n", " \n", " atm, xmap1 = self.calculate(\n", " atm,\n", " atm_profile_type,\n", " atm_profile_idx,\n", " chunk\n", " )\n", " \n", " forward_model.AtmosphereX = atm\n", " xmap[self.state_vector_slice, ipar, 0:atm.NP] = xmap1\n", " \n", "```\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }