.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_analysis/dde.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_analysis_dde.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_analysis_dde.py:


Analysis of Delayed Differential Equations
==========================================

In this tutorial, you will learn how to implement and solve delayed differential equation (DDE) systems via PyRates.
DDEs take the general form

.. math::
    \dot{\mathbf{y}} = \mathbf{f}(\mathbf{y}, t, \mathbf{y_{    au}})

with state vector :math:`\mathbf{y}`, vector field :math:`f`, time point :math:`t` and state history
:math:`\mathbf{y_{      au}} = {\mathbf{y}(\tau) : \tau \leq t}`.
Here, we will address the case of implementing and solving a DDE system with discrete delays of the form

.. math::
    \dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y}(t), \mathbf{y}(t-\tau)),

with discrete delay :math:`\tau`. As a specific example, we will use a modified version of the Van der Pol oscillator
with additional delayed feedback:

.. math::
        \dot x(t) &= z(t), 

        \dot z(t) &= \mu (1-x(t)^2) z - x(t) + k x(t-\tau),

with delayed feedback strength k. For a detailed introduction of the Van der Pol oscillator model, see the gallery
example in the model introduction section.

Below, we will show how such a DDE model can be implemented and how it can be solved via the third-party
DDE solver `ddeint`.

.. GENERATED FROM PYTHON SOURCE LINES 31-36

.. code-block:: Python

    import numpy as np

    from pyrates import CircuitTemplate, OperatorTemplate, clear
    import matplotlib.pyplot as plt


.. GENERATED FROM PYTHON SOURCE LINES 37-43

Step 1: DDE definition
^^^^^^^^^^^^^^^^^^^^^^

DDEs can be defined in two different ways in PyRates. The first way would be the usage of the `past(y, tau)` function
call that indicates to PyRates that :math:`y(t-tau)` should be evaluated instead of :math:`y(t)`. An example of how to
implement an `OperatorTemplate` representing the delayed Van der Pol equations is provided below.

.. GENERATED FROM PYTHON SOURCE LINES 44-63

.. code-block:: Python


    # define parameters
    k = -2.0
    tau = 5.0

    # define operator
    eqs = [
        "x' = z",
        "z' = mu*(1-x**2)*z - x + k*past(x,tau)"
    ]
    variables = {
        'x': 'output(0.0)',
        'z': 'variable(1.0)',
        'mu': 1.0,
        'k': k,
        'tau': tau
    }
    op = OperatorTemplate(name='vdp_delayed', equations=eqs, variables=variables)


.. GENERATED FROM PYTHON SOURCE LINES 64-72

This `OperatorTemplate` could then be used to define a PyRates model (see the examples in the model introduction
section of the gallery for a tutorial on how to do that).

Alternatively, it is also possible to add an edge to a `CircuitTemplate`, including a discrete delay. This edge will
then automatically translated into a corresponding `past` call by PyRates. Note that this only works if the edge
source variable is a variable defined by a differential equation. So in the case of the Van der Pol model, both
:math:`x` and :math:`z` would be valid variables for such an edge definition. Below, we show how to do that using the
Van der Pol model that is already implemented in PyRates:

.. GENERATED FROM PYTHON SOURCE LINES 72-76

.. code-block:: Python


    vdp = CircuitTemplate.from_yaml("model_templates.oscillators.vanderpol.vdp")
    vdp.update_template(edges=[('p/vdp_op/x', 'p/vdp_op/inp', None, {'weight': k, 'delay': tau})], in_place=True)


.. GENERATED FROM PYTHON SOURCE LINES 77-81

Step 2: DDE simulation
^^^^^^^^^^^^^^^^^^^^^^

We can perform a simple simulation using a forward Euler algorithm to solve the DDE as follows:

.. GENERATED FROM PYTHON SOURCE LINES 82-95

.. code-block:: Python


    # define simulation time
    T = 100.0
    step_size = 1e-2

    # solve DDE via forward Euler
    results = vdp.run(step_size=step_size, simulation_time=T, outputs={'x': 'p/vdp_op/x'}, solver='euler',
                      in_place=False, clear=True)

    # plot resulting signal
    plt.plot(results)
    plt.show()


.. GENERATED FROM PYTHON SOURCE LINES 96-101

To use a different solver, for example a Runge-Kutta algorithm with automatic step-size adaptation, you would either
have to use :code:`backend='julia'` and :code:`solver='julia_dde'` as options to the :code:`run` method, or you could
generate the run function instead and use any DDE solver available in the backend of your choice. Below, we
demonstrate how to solve the DDE via the `ddeint` Python package. The latter is not part of the PyRates prerequisites
and thus has to be manually installed via :code:`pip install ddeint`.

.. GENERATED FROM PYTHON SOURCE LINES 101-126

.. code-block:: Python


    from ddeint import ddeint

    # generate DDE function
    func, args, _, _ = vdp.get_run_func(func_name='vdp_dde_run', file_name='vdp_dde', solver='scipy', step_size=step_size,
                                        vectorize=False, backend='default', in_place=False)

    # define ddeint wrapper function
    t, y, hist, *fargs = args
    def dde_run(y, t):
        return func(t, y(t), y, *fargs)

    # solve DDE via ddeint
    eval_time = np.linspace(0, T, num=int(T/step_size))
    res2 = ddeint(func=dde_run, g=hist, tt=eval_time)

    # compare results
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(eval_time, res2[:, 0])
    ax.plot(results)
    ax.set_xlabel('time')
    ax.set_ylabel('x')
    plt.legend(['euler', 'ddeint'])
    plt.show()


.. GENERATED FROM PYTHON SOURCE LINES 127-132

As can be seen, the model dynamics agree qualitatively, with quantitative differences being caused by the different
solving algorithms. The most important feature of the DDE function generated by PyRates is that it contains a call to
a history function in order to calculate the delayed state variable versions. These history function calls vary in
their syntax between backends in order to allow to interface the most common DDE solvers in each programming language.
In the case of the default backend, the generated DDE function looks as follows:

.. GENERATED FROM PYTHON SOURCE LINES 132-138

.. code-block:: Python


    f = open('vdp_dde.py', 'r')
    print('')
    print(f.read())
    f.close()


.. GENERATED FROM PYTHON SOURCE LINES 139-141

In order to interface any given DDE solver, a wrapper function as defined above may have to be provided to the solver
to makes sure that the correct variable is passed to the :code:`hist` argument.

.. GENERATED FROM PYTHON SOURCE LINES 141-143

.. code-block:: Python


    clear(vdp)


.. _sphx_glr_download_auto_analysis_dde.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: dde.ipynb <dde.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: dde.py <dde.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_