Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CLUE tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does CLUE do?\n",
    "\n",
    "CLUE (Constrained LUmping for differential Eqiations) implements in Python an algorithm that takes as **input**\n",
    "* a system of ODEs with polynomial or rational right-hand side\n",
    "* a list of linear combinations of the unknown functions to be preserved (*observables*)\n",
    "\n",
    "and **returns** the maximal exact reduction of the system by a linear transformation that preserves the given combinations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic usage\n",
    "\n",
    "We will demonstrate CLUE on the following system\n",
    "\n",
    "$\\begin{cases} \\dot{x}_1  = x_2^2 + 4x_2x_3 + 4x_3^2,\\\\ \\dot{x}_2  =  4x_3 - 2x_1,\\\\ \\dot{x}_3  = x_1 + x_2 \\end{cases}$\n",
    "\n",
    "with the combination to preserve being just $x_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. import relevant functions from sympy and the function that does lumping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import vring, QQ\n",
    "from clue import FODESystem, SparsePolynomial, RationalFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Introduce the variables $x_1, x_2, x_3$ by defining the ring of polynomials in these variables (QQ refers to the coefficients being rational numbers; for other options, see below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = vring([\"x1\", \"x2\", \"x3\"], QQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Construct a list of right-hand sides of the ODE. The right-hand sides must be in the same order as the variables in the definition of the ring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode = FODESystem([\n",
    "    x2**2 + 4 * x2 * x3 + 4 * x3**2, # derivative of x1\n",
    "    4 * x3 - 2 * x1,                 # derivative of x2\n",
    "    x1 + x2                          # derivative of x3\n",
    "    ], variables=['x1','x2','x3'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Call method `lumping` from the system poviding the combinations to preserve, that is, `[x1]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1\n",
      "y1 = x2 + 2*x3\n",
      "Lumped system:\n",
      "y0' = y1**2\n",
      "y1' = 2*y1\n"
     ]
    }
   ],
   "source": [
    "ode.lumping([x1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computation shows that, in the new variables,\n",
    "\n",
    "$y_1 = x_1 \\quad \\text{ and }\\quad y_2 = x_2 + 2x_3,$\n",
    "\n",
    "the system will be reduced to \n",
    "\n",
    "$\\begin{cases} \\dot{y}_1 = y_2^2, \\\\ \\dot{y_2} = 2y_2. \\end{cases}$\n",
    "\n",
    "*Remark:* to supress the output, use `print_reduction=False` keyword argument. If you would like to print the original system, use `print_system=True`. Check the documentation of the method `lumping` for further information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This software also allows systems with rational functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "varnames = ['x','y']\n",
    "rational_ode = FODESystem([\n",
    "    RationalFunction.from_string(\"y/(x-y)\", varnames),\n",
    "    RationalFunction.from_string(\"x/(x-y)\", varnames)\n",
    "], variables = varnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-01-28 13:27:26 WARNING  Method [polynomial] selected but input is not SparsePolynomial. Trying random...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x + -y\n",
      "Lumped system:\n",
      "y0' = (-1)/(1)\n"
     ]
    }
   ],
   "source": [
    "rational_ode.lumping([SparsePolynomial.from_string(\"x-y\", varnames)]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subtleties and extra options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading models from \\*.ode files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CLUE accepts models in \\*.ode format which is especially convenient for large models (can be used for both ODE systems and chemical reaction networks). For more details in the format and how to convert SBML-encoded models to it, we refer to [this paper](https://doi.org/10.1007/978-3-030-31304-3_13).\n",
    "\n",
    "The model can be read from a file as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Ap', 'ApB', 'Au', 'AuB', 'B', 'r1', 'r2']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = FODESystem(file=\"./models/polynomial/model_example.ode\")\n",
    "system.variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observables can be specified as string expressions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = Ap + Au\n",
      "y1 = ApB + AuB\n",
      "y2 = B\n",
      "y3 = r1\n",
      "Lumped system:\n",
      "y0' = -3*y0*y2 + 4*y1\n",
      "y1' = 3*y0*y2 + -4*y1\n",
      "y2' = -3*y0*y2 + 4*y1\n",
      "y3' = 0\n"
     ]
    }
   ],
   "source": [
    "from clue.rational_function import SparsePolynomial\n",
    "\n",
    "obs = [\n",
    "    SparsePolynomial.from_string(\"Au + Ap\", system.variables),\n",
    "    SparsePolynomial.from_string(\"r1\", system.variables),\n",
    "]\n",
    "\n",
    "system.lumping(obs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rational numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the model involves rational numbers, e.g. $\\frac{2}{3}$, writing `2/3` will yield a floating point number and make the subsequent computation inexact.\n",
    "\n",
    "We recommend to define a rational number $\\frac{p}{q}$ as `QQ(p, q)`. For example, consider the system\n",
    "\n",
    "$\\begin{cases}\\dot{x}_1 = x_1 + \\frac{x_2}{3},\\\\ \\dot{x}_2 = -\\frac{2x_1}{3}, \\end{cases}$\n",
    "\n",
    "where the observable to preserve is $x_1 + x_2$.\n",
    "The maximal lumping can be found as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1 + x2\n",
      "Lumped system:\n",
      "y0' = 1/3*y0\n"
     ]
    }
   ],
   "source": [
    "R = vring([\"x1\", \"x2\"], QQ)\n",
    "ode = FODESystem([\n",
    "    x1 + QQ(1, 3) * x2,\n",
    "    QQ(-2, 3) * x1\n",
    "    ], variables=['x1','x2'])\n",
    "ode.lumping([x1 + x2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unknown scalar parameters\n",
    "\n",
    "It often happens that the system of interest involves unknown scalar parameters, for example:\n",
    "\n",
    "$\\begin{cases}\n",
    "  \\dot{x}_1 = ax_1 + bx_2,\\\\\n",
    "  \\dot{x}_2 = bx_1 + ax_2\n",
    "\\end{cases}$\n",
    "\n",
    "with observable $x_1 + x_2$.\n",
    "\n",
    "We will describe two different ways of applying CLUE to such models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convert the parameters into states with zero derivative**\n",
    "\n",
    "One can rewrite the system above as\n",
    "\n",
    "$\\begin{cases}\n",
    "  \\dot{x}_1 = a x_1 + b x_2,\\\\\n",
    "  \\dot{x}_2 = b x_1 + a x_2,\\\\\n",
    "  \\dot{a} = 0,\\\\\n",
    "  \\dot{b} = 0.\n",
    "\\end{cases}$\n",
    "\n",
    "And then apply CLUE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1 + x2\n",
      "y1 = a + b\n",
      "Lumped system:\n",
      "y0' = y0*y1\n",
      "y1' = 0\n"
     ]
    }
   ],
   "source": [
    "R = vring([\"x1\", \"x2\", \"a\", \"b\"], QQ)\n",
    "ode = FODESystem([\n",
    "    a * x1 + b * x2,\n",
    "    b * x1 + a * x2,\n",
    "    R(0), # this means that zero is interpreted as zero polynomial, not the zero number\n",
    "    R(0)\n",
    "], variables=['x1', 'x2', 'a', 'b'])\n",
    "ode.lumping([x1 + x2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that:\n",
    "* CLUE has also discovered a *reduction for the parameters* suggesting a single macroparameter $y_1 = a + b$\n",
    "* the new variables will be always linear combinations of the original variables and parameters. If you would like to search for a combination with coefficients involving parameters, use the next method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Include symbols $a$ and $b$ into the coefficient field**.\n",
    "\n",
    "Another option would be to set the coefficient field to be rational functions in the parameters. This method will generally be slower than the previous one because we have an improved version of algorithm for the rational coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1 + x2\n",
      "Lumped system:\n",
      "y0' = (a + b)*y0\n"
     ]
    }
   ],
   "source": [
    "# defining the coefficient field\n",
    "from sympy import FractionField\n",
    "coef_field = FractionField(QQ, [\"a\", \"b\"])\n",
    "a, b = coef_field.gens\n",
    "\n",
    "# defining a system over this field\n",
    "R = vring([\"x1\", \"x2\"], coef_field)\n",
    "ode = FODESystem([\n",
    "    a * x1 + b * x2,\n",
    "    b * x1 + a * x2\n",
    "], variables=['x1', 'x2'])\n",
    "ode.lumping([x1 + x2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As has been mentioned above, this method allows to find lumpings in which the new variables are expressed as combinations of the original ones with the *coefficients involving parameters*. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1\n",
      "y1 = x2 + (b/a)*x3\n",
      "Lumped system:\n",
      "y0' = (a)*y1\n",
      "y1' = y1\n"
     ]
    }
   ],
   "source": [
    "# defining the coefficient field\n",
    "from sympy import FractionField\n",
    "coef_field = FractionField(QQ, [\"a\", \"b\"])\n",
    "a, b = coef_field.gens\n",
    "\n",
    "# defining a system over this field\n",
    "R = vring([\"x1\", \"x2\", \"x3\"], coef_field)\n",
    "ode = FODESystem([\n",
    "    a * x2 + b * x3,\n",
    "    x2, \n",
    "    x3\n",
    "], variables = ['x1', 'x2', 'x3'])\n",
    "ode.lumping([x1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Irrational coefficients\n",
    "\n",
    "A system may involve irrational coefficients (e.g., $\\sqrt{2}$). For example, consider a reaction system from the paper [\"A general analysis of exact lumping in chemical kinetics\" (Example 2)](<https://doi.org/10.1016/0009-2509(89)85014-6>) by Li and Rabitz:\n",
    "\n",
    "$\\begin{cases}\n",
    "  \\dot{x}_1 = -2 x_1 - 2 x_1 x_2 + 4 x_3 x_4,\\\\\n",
    "  \\dot{x}_2 = -2 x_2 - 2 x_1 x_2 + 4 x_3 x_4,\\\\\n",
    "  \\dot{x}_3 = -2 x_3 - 4 x_3 x_4 + 2 x_1 x_2,\\\\\n",
    "  \\dot{x}_4 = -2 x_4 - 4 x_3 x_4 + 2 x_1 x_2,\\\\\n",
    "  \\dot{x}_5 = -x_5 + x_1 + 2 x_2 + \\sqrt{2} x_6,\\\\\n",
    "  \\dot{x}_6 = -\\sqrt{2} x_6 + 2 x_3 + x_5,\\\\\n",
    "  \\dot{x}_7 = -\\sqrt{2} x_7 + x_1 + x_8,\\\\\n",
    "  \\dot{x}_8 = -x_8 + 2 x_4 + \\sqrt{2} x_7\n",
    "\\end{cases}$\n",
    "\n",
    "with the observable $x_5 - \\sqrt{2} x_6$.\n",
    "In order to apply CLUE to the system, one can include $\\sqrt{2}$ into the field of coefficients as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New variables:\n",
      "y0 = x1\n",
      "y1 = x2\n",
      "y2 = x3\n",
      "y3 = x4\n",
      "y4 = x7 + (-sqrt(2)/2)*x8\n",
      "Lumped system:\n",
      "y0' = (-2)*y0 + (-2)*y0*y1 + (4)*y2*y3\n",
      "y1' = (-2)*y1 + (-2)*y0*y1 + (4)*y2*y3\n",
      "y2' = (-2)*y2 + (-4)*y2*y3 + (2)*y0*y1\n",
      "y3' = (-2)*y3 + (-4)*y2*y3 + (2)*y0*y1\n",
      "y4' = (-sqrt(2) - 1)*y4 + y0 + (-sqrt(2))*y3\n"
     ]
    }
   ],
   "source": [
    "from sympy import sqrt\n",
    "R = vring([f\"x{i + 1}\" for i in range(8)], QQ.algebraic_field(sqrt(2)))\n",
    "\n",
    "ode = FODESystem([\n",
    "    -2 * x1 - 2 * x1 * x2 + 4 * x3 * x4,\n",
    "    -2 * x2 - 2 * x1 * x2 + 4 * x3 * x4,\n",
    "    -2 * x3 - 4 * x3 * x4 + 2 * x1 * x2,\n",
    "    -2 * x4 - 4 * x3 * x4 + 2 * x1 * x2,\n",
    "    -x5 + x1 + 2 * x2 + sqrt(2) * x6,\n",
    "    -sqrt(2) * x6 + 2 * x3 + x5,\n",
    "    -sqrt(2) * x7 + x1 + x8,\n",
    "    -x8 + 2 * x4 + sqrt(2) * x7\n",
    "], variables = ['x1','x2','x3','x4','x5','x6','x7','x8'])\n",
    "res = ode.lumping([x8 - sqrt(2) * x7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}