{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A script to facilitate analysis of the LSD model in Section 6.5.1.\n",
    "\n",
    "#The Data: Sampling times, in hours\n",
    "data = [(1/12, 11.1), (1/4, 7.4), (1/2, 6.3), (1, 6.9), (2, 5), (4, 3.1),\n",
    "(8, 0.80), (1/12, 10.6), (1/4, 7.6), (1/2, 7), (1, 4.8), (2, 2.8),\n",
    "(4, 2.5), (8, 2), (1/12, 8.7), (1/4, 6.7), (1/2, 5.9), (1, 4.3),\n",
    "(2, 4.4), (8, 0.3), (1/12, 10.9), (1/4, 8.2), (1/2, 7.9), (1, 6.6),\n",
    "(2, 5.3), (4, 3.8), (8, 1.2), (1/12, 6.4), (1/4, 6.3), (1/2, 5.1),\n",
    "(1, 4.3), (2, 3.4), (4, 1.9), (8, 0.7)];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot of the data:\n",
    "plt1 = line([data[i] for i in range(7)],color='red',legend_label='Subject 1')\n",
    "plt2 = line([data[i] for i in range(7,14)],color='blue',legend_label='Subject 2')\n",
    "plt3 = line([data[i] for i in range(14,20)],color='black',legend_label='Subject 3')\n",
    "plt4 = line([data[i] for i in range(20,27)],color='green',legend_label='Subject 4')\n",
    "plt5 = line([data[i] for i in range(27,34)],color='purple',legend_label='Subject 5')\n",
    "p = plt1+plt2+plt3+plt4+plt5;\n",
    "p.set_legend_options(loc='upper center');\n",
    "p.axes_labels(['Time (hours)','Concentration (mg/ml)'])\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The ODE Model: We use cP(t) for the plasma concentration and cT(t) as the tissue concentration.  We have cP(0) = 12.27 \n",
    "#for each subject and cT(0) = 0. The model for cP(t) and cT(t) of Section 6.5.1 (equation (6.68)) is\n",
    "var('t ka kb ke');\n",
    "cP = function('cP')(t);\n",
    "cT = function('cT')(t);\n",
    "deP = diff(cP,t) ==  (-kb - ke)*cP + 0.706*ka*cT;\n",
    "deT = diff(cT,t) == 1.407*kb*cP - ka*cT;\n",
    "#where ka, kb, and ke are the constants to be determined from the concentration data.\n",
    "\n",
    "#Unfortunately, Sage's desolve_system() will not solve this system, nor will the Laplace transform commands work \n",
    "#(the symbolic constants ka, kb, and ke in the ODEs cause problems.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Instead, let's use eigenvalue/vector techniques to solve this ODE system. The system can be expressed as a homogeneous \n",
    "#linear system of the form x' = A*x where\n",
    "A = matrix(SR,[[-(kb+ke),0.706*ka],[1.407*kb,-ka]]); #Matrix over SR, the reals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Eigenvalues/vectors are\n",
    "eigs = A.eigenvectors_right() #Compute right eigenvectors, A*v = lambda*v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvalues are\n",
    "lambda1 = eigs[0][0];\n",
    "lambda2 = eigs[1][0];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be defined as vectors as\n",
    "v1 = vector(eigs[0][1][0]);\n",
    "v2 = vector(eigs[1][1][0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a matrix P with the eigenvectors as the columns:\n",
    "P = transpose(matrix([v1, v2]))\n",
    "#as was done in equation (6.18) in the text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let\n",
    "\n",
    "x0 = vector(SR,[12.27,0]); #Symbolic vector with real entries\n",
    "\n",
    "#denote the initial condition vector. As per the procedure in Section 6.2.2, we can construct the solution to this ODE\n",
    "#system by first solving Pc = x0 for c, which in Sage takes the form\n",
    "c = P\\x0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From (6.17) in the text the solution is then given by \"xsol\" where\n",
    "var('t');\n",
    "xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assign the solution components to cPsol(t) and cTsol(t):\n",
    "cPsol(t,ka,kb,ke) = xsol[0];\n",
    "cTsol(t,ka,kb,ke) = xsol[1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We could form a sum of squares and minimize using Calculus 3 techniques, but in Sage it is easier\n",
    "#to use Sage's built-in command \"find_fit\" to find the least-squares model\n",
    "#fit to the data. See the worksheet \"paramest_demo.sage\"."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}