{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "74541811",
   "metadata": {},
   "source": [
    "The adaptive integrator\n",
    "===================\n",
    "\n",
    "The ``taylor_adaptive`` class provides an easy-to-use interface to heyoka.py's\n",
    "main capabilities. Objects of this class can be constructed from a system\n",
    "of ODEs and a set of initial conditions (plus a number of optional configuration parameters\n",
    "with - hopefully - sensible defaults). Methods are provided to\n",
    "propagate in time the state of the system, either step-by-step or by specifying\n",
    "time limits.\n",
    "\n",
    "Let's see how we can use ``taylor_adaptive`` to integrate the ODE\n",
    "system of the [simple pendulum](https://en.wikipedia.org/wiki/Pendulum_(mathematics)),\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "x^\\prime = v \\\\\n",
    "v^\\prime = -9.8 \\sin x\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "with initial conditions\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "x\\left( 0 \\right) = 0.05 \\\\\n",
    "v\\left( 0 \\right) = 0.025\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Construction\n",
    "-------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7b0ebb82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Tolerance               : 2.2204460492503131e-16\n",
       "High accuracy           : false\n",
       "Compact mode            : false\n",
       "Taylor order            : 20\n",
       "Dimension               : 2\n",
       "Time                    : 0.0000000000000000\n",
       "State                   : [0.050000000000000003, 0.025000000000000001]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import heyoka as hy\n",
    "\n",
    "# Create the symbolic variables x and v.\n",
    "x, v = hy.make_vars(\"x\", \"v\")\n",
    "\n",
    "# Create the integrator object.\n",
    "ta = hy.taylor_adaptive(\n",
    "                        # Definition of the ODE system:\n",
    "                        # x' = v\n",
    "                        # v' = -9.8 * sin(x)\n",
    "                        sys = [(x, v),\n",
    "                         (v, -9.8 * hy.sin(x))],\n",
    "                        # Initial conditions for x and v.\n",
    "                        state = [0.05, 0.025])\n",
    "\n",
    "ta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f1f15b3",
   "metadata": {},
   "source": [
    "After creating the symbolic variables ``x`` and ``v``, we\n",
    "construct an instance of ``taylor_adaptive`` called ``ta``.\n",
    "By default, ``taylor_adaptive`` operates using double-precision arithmetic. As (mandatory) construction arguments, we pass in the system of differential equations and a set\n",
    "of initial conditions for ``x`` and ``v`` respectively.\n",
    "\n",
    "By default, the error tolerance of an adaptive integrator is set to the\n",
    "machine epsilon, which, for double precision, is $\\sim 2.2\\times10^{-16}$.\n",
    "From this value, heyoka.py deduces an optimal Taylor order of 20, as indicated\n",
    "by the screen output. ``taylor_adaptive`` manages its own copy of the state vector and the time variable.\n",
    "Both the state vector and the time variable are updated automatically by the timestepping\n",
    "methods. Note also how, by default, the time variable is initially set to zero.\n",
    "\n",
    "Single timestep\n",
    "---------------\n",
    "\n",
    "Let's now try to perform a single integration timestep:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bd2a9ac0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome : taylor_outcome.success\n",
      "Timestep: 0.21605277478009474\n"
     ]
    }
   ],
   "source": [
    "# Perform a single step.\n",
    "oc, h = ta.step()\n",
    "\n",
    "# Print the outcome flag and the timestep used.\n",
    "print(\"Outcome : {}\".format(oc))\n",
    "print(\"Timestep: {}\".format(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b24e2dd0",
   "metadata": {},
   "source": [
    "First, we invoke the ``step()`` method, which returns a pair of values.\n",
    "The first value is a status flag indicating the outcome of the integration timestep,\n",
    "while the second value is the step size that was selected by heyoka.py in order\n",
    "to respect the desired error tolerance.\n",
    "\n",
    "Let's also print again the integrator object to screen in order to inspect how state and time have changed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f75fed0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Tolerance               : 2.2204460492503131e-16\n",
       "High accuracy           : false\n",
       "Compact mode            : false\n",
       "Taylor order            : 20\n",
       "Dimension               : 2\n",
       "Time                    : 0.21605277478009474\n",
       "State                   : [0.043996448369926382, -0.078442455470687983]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fda0d15",
   "metadata": {},
   "source": [
    "It is also possible to perform a single timestep backward in time\n",
    "via the ``step_backward()`` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dc54bdd2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome : taylor_outcome.success\n",
      "Timestep: -0.21312300047513288\n"
     ]
    }
   ],
   "source": [
    "# Perform a step backward.\n",
    "oc, h = ta.step_backward()\n",
    "\n",
    "# Print the outcome flag and the timestep used.\n",
    "print(\"Outcome : {}\".format(oc))\n",
    "print(\"Timestep: {}\".format(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78144bc3",
   "metadata": {},
   "source": [
    "The ``step()`` method can also be called with an argument representing\n",
    "the maximum step size ``max_delta_t``: if the adaptive timestep\n",
    "selected by heyoka.py is larger (in absolute value) than ``max_delta_t``,\n",
    "then the timestep will be clamped to ``max_delta_t``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0186f288",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome : taylor_outcome.time_limit\n",
      "Timestep: 0.01\n",
      "\n",
      "Outcome : taylor_outcome.time_limit\n",
      "Timestep: -0.02\n"
     ]
    }
   ],
   "source": [
    "# Perform a step forward in time, clamping\n",
    "# the timestep size to 0.01.\n",
    "oc, h = ta.step(max_delta_t = 0.01)\n",
    "\n",
    "# Print the outcome flag and the timestep used.\n",
    "print(\"Outcome : {}\".format(oc))\n",
    "print(\"Timestep: {}\\n\".format(h))\n",
    "\n",
    "# Perform a step backward in time, clamping\n",
    "# the timestep size to 0.02.\n",
    "oc, h = ta.step(max_delta_t = -0.02)\n",
    "\n",
    "# Print the outcome flag and the timestep used.\n",
    "print(\"Outcome : {}\".format(oc))\n",
    "print(\"Timestep: {}\".format(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2411cb35",
   "metadata": {},
   "source": [
    "Note that the integration outcome is now ``time_limit``, instead of ``success``.\n",
    "\n",
    "Before moving on, we need to point out an important caveat when using the single\n",
    "step functions:\n",
    "\n",
    "> **WARNING**: if the exact solution of the ODE system can be expressed as a polynomial function\n",
    "> of time, the automatic timestep deduction algorithm may return a timestep of infinity.\n",
    "> This is the case, for instance, when integrating the rectilinear motion of a free\n",
    "> particle or the constant-gravity free-fall equation. In such cases, the step functions\n",
    "> should be called with a finite ``max_delta_t`` argument, in order to clamp the timestep\n",
    "> to a finite value.\n",
    ">\n",
    "> Note that the ``propagate_*()`` functions (described below)\n",
    "> are not affected by this issue.\n",
    "\n",
    "Accessing state and time\n",
    "------------------------\n",
    "\n",
    "It is possible to read from and write to both the time variable and the state\n",
    "vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ebf10bac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current time        : -0.007070225695038143\n",
      "Current state vector: [0.04981102 0.02845657]\n",
      "\n",
      "Current time        : 0.0\n",
      "Current state vector: [0.05  0.025]\n"
     ]
    }
   ],
   "source": [
    "# Print the current time.\n",
    "print(\"Current time        : {}\".format(ta.time))\n",
    "\n",
    "# Print out the current state vector.\n",
    "print(\"Current state vector: {}\\n\".format(ta.state))\n",
    "\n",
    "# Reset the time and state to the initial values.\n",
    "ta.time = 0.\n",
    "ta.state[:] = [0.05, 0.025]\n",
    "\n",
    "# Print them again.\n",
    "print(\"Current time        : {}\".format(ta.time))\n",
    "print(\"Current state vector: {}\".format(ta.state))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351f0c4a",
   "metadata": {},
   "source": [
    "Note that the time is stored as a scalar value, while the state is stored as a NumPy array.\n",
    "\n",
    "Because of technical reasons related to the management of the lifetime of arrays when interacting with the underlying heyoka C++ library, it is **not** possible to directly set the state via the syntax\n",
    "\n",
    "```python\n",
    ">>> ta.state = [0.05, 0.025] # Won't work!\n",
    "```\n",
    "\n",
    "Thus, the array assignment syntax ``ta.state[:] = ...`` must be used instead. Similarly, it is also possible to set directly the values of the components of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b62eb907",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta.state[0] = 0.05\n",
    "ta.state[1] = 0.025"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b3e66e",
   "metadata": {},
   "source": [
    "Time-limited propagation\n",
    "------------------------\n",
    "\n",
    "In addition to the step-by-step integration methods,\n",
    "``taylor_adaptive`` also provides methods to propagate\n",
    "the state of the system for a specified amount of time.\n",
    "These methods are called ``propagate_for()`` and\n",
    "``propagate_until()``: the former integrates\n",
    "the system for a specified amount of time, the latter\n",
    "propagates the state up to a specified epoch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a84060d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome      : taylor_outcome.time_limit\n",
      "Min. timestep: 0.20213323505293765\n",
      "Max. timestep: 0.21813566576411725\n",
      "Num. of steps: 24\n",
      "Current time : 5.0\n",
      "\n",
      "Outcome      : taylor_outcome.time_limit\n",
      "Min. timestep: 0.20212172864807665\n",
      "Max. timestep: 0.2181392923080563\n",
      "Num. of steps: 72\n",
      "Current time : 20.0\n"
     ]
    }
   ],
   "source": [
    "# Propagate for 5 time units.\n",
    "status, min_h, max_h, nsteps, _ = ta.propagate_for(delta_t = 5.)\n",
    "\n",
    "print(\"Outcome      : {}\".format(status))\n",
    "print(\"Min. timestep: {}\".format(min_h))\n",
    "print(\"Max. timestep: {}\".format(max_h))\n",
    "print(\"Num. of steps: {}\".format(nsteps))\n",
    "print(\"Current time : {}\\n\".format(ta.time))\n",
    "\n",
    "# Propagate until t = 20.\n",
    "status, min_h, max_h, nsteps, _ = ta.propagate_until(t = 20.)\n",
    "\n",
    "print(\"Outcome      : {}\".format(status))\n",
    "print(\"Min. timestep: {}\".format(min_h))\n",
    "print(\"Max. timestep: {}\".format(max_h))\n",
    "print(\"Num. of steps: {}\".format(nsteps))\n",
    "print(\"Current time : {}\".format(ta.time))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adfbf18a",
   "metadata": {},
   "source": [
    "The time-limited propagation methods return\n",
    "a tuple of 5 values, which represent, respectively:\n",
    "\n",
    "* the outcome of the integration (which will always be\n",
    "  ``time_limit``, unless error conditions arise),\n",
    "* the minimum and maximum integration timesteps\n",
    "  that were used in the propagation,\n",
    "* the total number of steps that were taken,\n",
    "* the [continuous output](<./Dense output.ipynb>) function object,\n",
    "  if requested (off by default).\n",
    "\n",
    "The time-limited propagation methods can be used\n",
    "to propagate both forward and backward in time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "aa877a4e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Outcome      : taylor_outcome.time_limit\n",
      "Min. timestep: 0.20207792808238695\n",
      "Max. timestep: 0.21818982934810394\n",
      "Num. of steps: 97\n",
      "Current time : 0.0\n",
      "\n",
      "Tolerance               : 2.2204460492503131e-16\n",
      "High accuracy           : false\n",
      "Compact mode            : false\n",
      "Taylor order            : 20\n",
      "Dimension               : 2\n",
      "Time                    : 0.0000000000000000\n",
      "State                   : [0.050000000000000044, 0.024999999999999991]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Propagate back to t = 0.\n",
    "status, min_h, max_h, nsteps, _ = ta.propagate_until(t = 0.)\n",
    "\n",
    "print(\"Outcome      : {}\".format(status))\n",
    "print(\"Min. timestep: {}\".format(min_h))\n",
    "print(\"Max. timestep: {}\".format(max_h))\n",
    "print(\"Num. of steps: {}\".format(nsteps))\n",
    "print(\"Current time : {}\\n\".format(ta.time))\n",
    "\n",
    "print(ta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d012675f-273d-4387-b780-c841b0891104",
   "metadata": {},
   "source": [
    "Note also that the time-limited propagation methods will stop\n",
    "integrating if a non-finite value is detected in the state vector\n",
    "at the end of the timestep. In such case, the outcome of the\n",
    "integration will be ``err_nf_state``.\n",
    "\n",
    "The ``propagate_for()`` and ``propagate_until()`` methods\n",
    "can be invoked with two additional optional keyword arguments:\n",
    "\n",
    "- ``max_delta_t``: similarly to the ``step()`` function, this value\n",
    "  represents the maximum timestep size in absolute value;\n",
    "- ``callback``: this is a callable which will be invoked at the end of\n",
    "  each timestep, with the integrator object as only argument. If the callback returns ``True`` then the integration\n",
    "  will continue after the invocation of the callback, otherwise the integration\n",
    "  will be interrupted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ed030c82-59c9-4b45-8abe-75c90818cc13",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current time: 0.1\n",
      "Current time: 0.2\n",
      "Current time: 0.30000000000000004\n",
      "Current time: 0.4\n",
      "Current time: 0.5\n"
     ]
    }
   ],
   "source": [
    "# Propagate to t = .5 using a max_delta_t and\n",
    "# providing a callback that prints the current time.\n",
    "\n",
    "# The callback.\n",
    "def cb(ta):\n",
    "    print(\"Current time: {}\".format(ta.time))\n",
    "    \n",
    "    return True\n",
    "\n",
    "ta.propagate_until(t = .5, max_delta_t = .1, callback = cb);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45f92420",
   "metadata": {},
   "source": [
    "Propagation over a time grid\n",
    "----------------------------\n",
    "\n",
    "Another way of propagating the state of a system in a ``taylor_adaptive``\n",
    "integrator is over a time grid. In this mode, the integrator\n",
    "uses [dense output](<./Dense output.ipynb>) to compute the state of the system\n",
    "over a grid of time coordinates provided by the user. If the grid is denser\n",
    "than the typical timestep size, this can be noticeably more efficient than\n",
    "repeatedly calling ``propagate_until()`` on the grid points, because\n",
    "propagating the system state via dense output is much faster than taking\n",
    "a full integration step.\n",
    "\n",
    "Let's see a simple usage example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b07369e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<taylor_outcome.time_limit: -4294967299>,\n",
       " 0.20291845444801257,\n",
       " 0.216140019757225,\n",
       " 5,\n",
       " array([[ 0.05003035, -0.024398  ],\n",
       "        [ 0.04519961, -0.07142727],\n",
       "        [ 0.03597685, -0.11152037],\n",
       "        [ 0.02325783, -0.14078016],\n",
       "        [ 0.00827833, -0.15635952],\n",
       "        [-0.00750582, -0.15674117],\n",
       "        [-0.02256041, -0.14188793],\n",
       "        [-0.03542229, -0.11324639],\n",
       "        [-0.04484178, -0.07360369],\n",
       "        [-0.04990399, -0.02681336]]))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Reset the time and state to the initial values.\n",
    "ta.time = 0.\n",
    "ta.state[:] = [0.05, 0.025]\n",
    "\n",
    "# Propagate over a time grid from 0 to 1\n",
    "# at regular intervals.\n",
    "out = ta.propagate_grid(grid = [.1, .2, .3, .4, .5, .6, .7, .8, .9, 1.])\n",
    "out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fc39da8",
   "metadata": {},
   "source": [
    "``propagate_grid()`` takes in input a grid of time points,\n",
    "and returns a tuple of 5 values. The first 4 values are the same\n",
    "as in the other ``propagate_*()`` functions:\n",
    "\n",
    "* the outcome of the integration,\n",
    "* the minimum and maximum integration timesteps\n",
    "  that were used in the propagation,\n",
    "* the total number of steps that were taken.\n",
    "\n",
    "The fifth value returned by ``propagate_grid()`` is a 2D array containing\n",
    "the state of the system at the time points in the grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "70aac3ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "State at t = 0.4: [ 0.02325783 -0.14078016]\n"
     ]
    }
   ],
   "source": [
    "# Print the state at t = 0.4 (index 3 in the time grid).\n",
    "print(\"State at t = 0.4: {}\".format(out[4][3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8824388",
   "metadata": {},
   "source": [
    "As you can see from the screen output above (where we printed the ``out`` variable), the use of ``propagate_grid()`` resulted in 5 integration timesteps being taken. Had we used ``propagate_until()``, we would have needed 10 integration timesteps to obtain the same result.\n",
    "\n",
    "There are no special requirements on the time values in the grid (apart from the fact that they must be finite and ordered monotonically).\n",
    "\n",
    "The ``propagate_grid()`` method\n",
    "can be invoked with two additional optional keyword arguments:\n",
    "\n",
    "- ``max_delta_t``: similarly to the ``step()`` function, this value\n",
    "  represents the maximum timestep size in absolute value;\n",
    "- ``callback``: this is a callable which will be invoked at the end of\n",
    "  each timestep, with the integrator object as only argument. If the callback returns ``True`` then the integration\n",
    "  will continue after the invocation of the callback, otherwise the integration\n",
    "  will be interrupted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4eb16baf-5196-4e8d-9a57-de2ed21ba935",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current time: 1.2000000000000002\n",
      "Current time: 1.3\n",
      "Current time: 1.4000000000000001\n",
      "Current time: 1.5\n",
      "Current time: 1.6\n",
      "Current time: 1.7000000000000002\n",
      "Current time: 1.8\n",
      "Current time: 1.9000000000000001\n",
      "Current time: 2.0\n"
     ]
    }
   ],
   "source": [
    "# Propagate over a grid using a max_delta_t and\n",
    "# providing a callback that prints the current time.\n",
    "\n",
    "# The callback.\n",
    "def cb(ta):\n",
    "    print(\"Current time: {}\".format(ta.time))\n",
    "    \n",
    "    return True\n",
    "\n",
    "ta.propagate_grid(grid = [1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.],\n",
    "                  max_delta_t = .1, callback = cb);"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}