added new use example,

updated some of the plotting functions
pyrates-neuroscience · Oct 17, 2022 · b74e3a0 · b74e3a0
1 parent 4bedfeb
commit b74e3a0
Show file tree

Hide file tree

Showing 7 changed files with 163 additions and 63 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -4,6 +4,14 @@ Changelog
 0.7
 ---
 
+0.7.2
+~~~~~
+
+- added use example for the QIF-SFA model to the documentation
+- improved support for `pandas.DataFrames` as the main results storage data type
+- added the pyrates model template as an attribute to the `ODESystem`
+- added the option of clearing all generated fortran/auto files via the `ODESystem.close_session()` method
+
 0.7.1
 ~~~~~
 

diff --git a/documentation/doc/source/index.rst b/documentation/doc/source/index.rst
@@ -32,6 +32,13 @@ Content
    simple_example
    changelog
 
+.. toctree::
+    :maxdepth: 1
+    :numbered:
+    :caption: Use Examples
+    :hidden:
+
+    use_examples
 
 Reference
 ---------

diff --git a/documentation/doc/source/simple_example.rst b/documentation/doc/source/simple_example.rst
@@ -1,6 +1,6 @@
-************************************************
-Parameter continuation and bifurcation detection
-************************************************
+*********************************************************
+Example: Parameter continuation and bifurcation detection
+*********************************************************
 
 In this tutorial, you will learn how to perform a 1D
 `numerical parameter continuation <http://www.scholarpedia.org/article/Numerical_analysis#Numerical_solution_of_
@@ -11,6 +11,7 @@ the quadratic integrate-and-fire population model [1]_, a detailed introduction
 introductions example gallery. The dynamic equations of this model read the following:
 
 .. math::
+
     \tau \dot r &= \frac{\Delta}{\pi\tau} + 2 r v, \\
     \tau \dot v &= v^2 +\bar \eta + I(t) + J r \tau - (\pi r \tau)^2,
 
@@ -40,7 +41,7 @@ Part 1: Creating an ODESystem Instance
 
 In this first part, we will be concerned with how to create a model representation that is compatible with auto-07p,
 which is the software that is used for parameter continuations and bifurcation analysis in PyCoBi [2]_.
-To this end, we make use of the code generation software `PyRates <>`_, which allows us to generate the Fortran files
+To this end, we make use of the code generation software `PyRates <https://github.com/pyrates-neuroscience/PyRates>`_, which allows us to generate the Fortran files
 required to run auto-07p from `YAML` model definition files.
 The QIF mean-field model comes pre-implemented with `PyRates`, so we can simply point to the ready-to-use `YAML` definition file.
 

diff --git a/documentation/doc/source/use_examples.rst b/documentation/doc/source/use_examples.rst
@@ -0,0 +1,6 @@
+***************************************************
+Use Examples for Parameter Continuations via PyCoBi
+***************************************************
+
+.. include::
+   auto_examples/index.rst
diff --git a/documentation/use_examples/qif_sfa.py b/documentation/use_examples/qif_sfa.py
@@ -42,10 +42,7 @@
 References
 ^^^^^^^^^^
 
-.. [1] E. Montbrió, D. Pazó, A. Roxin (2015) *Macroscopic description for networks of spiking neurons.* Physical
-       Review X, 5:021028, https://doi.org/10.1103/PhysRevX.5.021028.
-
-.. [2] R. Gast, H. Schmidt, T.R. Knösche (2020) *A Mean-Field Description of Bursting Dynamics in Spiking Neural
+.. [1] R. Gast, H. Schmidt, T.R. Knösche (2020) *A Mean-Field Description of Bursting Dynamics in Spiking Neural
        Networks with Short-Term Adaptation.* Neural Computation 32 (9): 1615-1634, https://doi.org/10.1162/neco_a_01300.
 """
 
@@ -58,10 +55,12 @@
 from pycobi import ODESystem
 import numpy as np
 
-ode, _ = ODESystem.from_yaml("model_templates.neural_mass_models.qif.qif_sfa", auto_dir="~/PycharmProjects/auto-07p",
-                             node_vars={'p/qif_sfa_op/Delta': 2.0, 'p/qif_sfa_op/alpha': 1.0, 'p/qif_sfa_op/eta': 3.0},
-                             edge_vars=[('p/qif_sfa_op/r', 'p/qif_sfa_op/r_in', {'weight': 15.0*np.sqrt(2.0)})],
-                             NPR=100, NMX=30000)
+ode = ODESystem.from_yaml(
+    "model_templates.neural_mass_models.qif.qif_sfa", auto_dir="~/PycharmProjects/auto-07p",
+    node_vars={'p/qif_sfa_op/Delta': 2.0, 'p/qif_sfa_op/alpha': 1.0, 'p/qif_sfa_op/eta': 3.0},
+    edge_vars=[('p/qif_sfa_op/r', 'p/qif_sfa_op/r_in', {'weight': 15.0*np.sqrt(2.0)})],
+    NPR=100, NMX=30000
+)
 
 # %%
 # In the code above, we adjusted some default parameters of the model in accordance with the parameters reported
@@ -74,12 +73,12 @@
 plt.show()
 
 # %%
-# If the system didn't converge yet, try increasing the paramter :code:`NMX` provided to the :code:`ODESystem.from_yaml`
+# If the system didn't converge yet, try increasing the parameter :code:`NMX` provided to the :code:`ODESystem.from_yaml`
 # method, which controls the number of time integration steps, until the system is assumed to have converged.
 
 # %%
-# Step 2: Parameter continuation of the excitability
-# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+# Step 2: 1D parameter continuation of a steady-state solution
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
 # A 1D parameter continuation in :math:`\bar \eta` can be carried out as follows:
 
@@ -94,7 +93,12 @@
 # %%
 # In the above code, we provided all the `auto-07p` meta parameters as keyword arguments. As an alternative, you can
 # create constants file with name :code:`c.name` and provide it to :code:`ODESystem.run` via keyword argument
-# :code:`c=name`.
+# :code:`c=name`. See the `auto-07p documentation <https://github.com/auto-07p/auto-07p/doc>`_ for a detailed explanation
+# of each of these parameters.
+# The following keyword arguments are specific to `PyCoBi` and deviate from standard `auto-07p` syntax:
+# :code:`bidirectional = True` leads to automatic continuation of the solution branch into both directions of the
+# continuation parameter :math:`\bar \eta`. :code:`origin` indicates the key of the solution branch, where a solution with
+# the name :code:`starting_point='EP'` exists, which is used as the starting point of the parameter continuation.
 # We can plot the results of the 1D parameter continuation, to examine the results.
 
 ode.plot_continuation("PAR(4)", "U(1)", cont="eta")
@@ -107,6 +111,11 @@
 # marked by the green circles (`Hopf bifurcations <http://www.scholarpedia.org/article/Andronov-Hopf_bifurcation>`_)
 # and grey triangles (`fold bifurcations <http://www.scholarpedia.org/article/Saddle-node_bifurcation>`_),
 # at which the stability of the solution branch changes.
+
+# %%
+# Step 3: Branch switching at a Hopf bifurcation to the limit cycle solutions
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+#
 # We know from bifurcation theory that the Hopf bifurcation at :math:`\bar \eta \approx 4` is a subcritical Hopf
 # bifurcation, and that an unstable limit cycle branch must emerge from that bifurcation. We can switch to this branch
 # and continue the limit cycle solution curve as follows.
@@ -119,6 +128,28 @@
 # %%
 # Let's visualize the 1D bifurcation diagram again, this time with the limit cycle branch included:
 
-ode.plot_continuation("PAR(4)", "U(1)", cont="eta")
-ode.plot_continuation("PAR(4)", "U(1)", cont="eta_hopf")
+fig, ax = plt.subplots()
+ode.plot_continuation("PAR(4)", "U(1)", cont="eta", ax=ax)
+ode.plot_continuation("PAR(4)", "U(1)", cont="eta_hopf", ax=ax)
+plt.show()
+
+# %%
+# As we can see, the limit cycle solution curve that branches off the Hopf bifurcation is indeed unstable.
+# However, it undergoes a fold-of-limit-cycle bifurcation, giving rise to a stable limit cycle solution branch.
+# During continuation of the limit cycle branch, we set a user point at :math:`\bar \eta = 2.0` via
+# :code:`UZR={4: -2.0}`. In the code below, we integrate the system equations with respect to time to show that the
+# system dynamics are indeed governed by the stable limit cycle solution in that regime.
+
+ode.run(origin=hopf_cont, starting_point="UZ1", c="ivp", name="lc")
+ode.plot_continuation("PAR(14)", "U(1)", cont="lc")
 plt.show()
+
+# %%
+# This concludes our use example. As a final step, we clear all the temporary files we created and make sure all
+# working directory changes that have happened in the background during file creation and parameter continuation are
+# undone. This can be done via a single call:
+
+ode.close_session(clear_files=True)
+
+# %%
+# If you want to keep the Fortran files for future parameter continuations etc., just set :code:`clear_files=False`.
diff --git a/pycobi/pycobi.py b/pycobi/pycobi.py
@@ -1,9 +1,9 @@
 import os
 import matplotlib.pyplot as plt
 import numpy as np
-from pandas import DataFrame
+from pandas import DataFrame, MultiIndex
 from mpl_toolkits.mplot3d import Axes3D
-from pyrates import CircuitTemplate
+from pyrates import CircuitTemplate, clear
 from matplotlib.collections import LineCollection
 from mpl_toolkits.mplot3d.art3d import Line3DCollection
 from typing import Union, Any, Optional
@@ -52,6 +52,7 @@ def __init__(self, working_dir: str = None, auto_dir: str = None, init_cont: boo
                                     'BT': {'marker': 's', 'color': 'k'},
                                     'GH': {'marker': 'o', 'color': '#148F77'}
                                     }
+        self._temp = kwargs.pop("template", None)
 
         # perform initial continuation in time to ensure convergence to steady-state solution
         if init_cont:
@@ -63,12 +64,18 @@ def __getitem__(self, item):
         except KeyError:
             return self.results[self._results_map[item]]
 
-    def close_session(self):
+    @property
+    def pyrates_template(self):
+        return self._temp
+
+    def close_session(self, clear_files: bool = False, **kwargs):
+        if clear_files:
+            clear(self._temp, **kwargs)
         os.chdir(self._orig_dir)
 
     @classmethod
     def from_yaml(cls, path: str, working_dir: str = None, auto_dir: str = None, init_cont: bool = True,
-                  init_kwargs: dict = None, **kwargs) -> tuple:
+                  init_kwargs: dict = None, **kwargs):
         """
 
         Parameters
@@ -82,8 +89,8 @@ def from_yaml(cls, path: str, working_dir: str = None, auto_dir: str = None, ini
 
         Returns
         -------
-        tuple
-            `ODESystem` instance and `CircuitTemplate` instance
+        ODESystem
+            `ODESystem` instance.
         """
 
         # preparations
@@ -110,7 +117,7 @@ def from_yaml(cls, path: str, working_dir: str = None, auto_dir: str = None, ini
 
         # initialize ODESystem
         return cls(working_dir=working_dir, auto_dir=auto_dir, init_cont=init_cont, e=file_name, c="ivp",
-                   **init_kwargs), template
+                   template=template, **init_kwargs)
 
     def run(self, variables: list = None, params: list = None, get_stability: bool = True,
             get_period: bool = False, get_timeseries: bool = False, get_eigenvals: bool = False,
@@ -281,8 +288,8 @@ def merge(self, key: int, cont, summary: DataFrame, icp: tuple):
         start_old, start_new, end_old, end_new = points_old[0], points_new[0], points_old[-1], points_new[-1]
 
         # extract ICP values
-        icp_old = summary_old.at[end_old, f"PAR({icp[0]})"]
-        icp_new = summary.at[end_new, f"PAR({icp[0]})"]
+        icp_old = summary_old.loc[end_old, f"PAR({icp[0]})"][0]
+        icp_new = summary.loc[end_new, f"PAR({icp[0]})"][0]
 
         # connect starting points of both continuations and re-label points accordingly
         if icp_new < icp_old:
@@ -662,7 +669,7 @@ def plot_bifurcation_points(self, solution_types, x_vals, y_vals, ax, default_co
         plt.sca(ax)
 
         # draw bifurcation points
-        for bf, x, y in zip(solution_types, x_vals, y_vals):
+        for bf, x, y in zip(solution_types.values, x_vals.values, y_vals.values):
             if bf not in "EPMXRG" and bf not in ignore:
                 if bf in bf_styles:
                     m = bf_styles[bf]['marker']
@@ -714,43 +721,87 @@ def _create_summary(self, solution: Union[Any, dict], points: list, variables: l
         -------
         DataFrame
         """
-        summary = {}
+
+        columns_2d, columns_1d, indices, data_2d, data_1d = [], [], [], [], []
+        add_columns = True
         for point in points:
 
+            data_2d_tmp = []
+            data_1d_tmp = []
+
             # get solution
             solution_type, s = self.get_solution(cont=solution, point=point)
 
             if solution_type != 'No Label' and solution_type != 'MX':
 
-                summary[point] = {}
+                indices.append(point)
 
                 # extract variables and params from solution
                 var_vals = self.get_vars(s, variables, timeseries)
                 param_vals = self.get_params(s, params)
 
-                # store solution information in summary
-                summary[point]['bifurcation'] = solution_type
+                # store solution type
+                data_1d_tmp.append(solution_type)
+                if add_columns:
+                    columns_1d.append('bifurcation')
+
+                # store parameter and variable information
                 branch, _ = self.get_branch_info(s)
                 for var, val in zip(variables, var_vals):
-                    summary[point][var] = val
-                if len(var_vals) > len(variables) and timeseries:
-                    summary[point]['time'] = var_vals[-1]
+                    for i, v in enumerate(val):
+                        data_2d_tmp.append(v)
+                        if add_columns:
+                            columns_2d.append((var, i))
                 for param, val in zip(params, param_vals):
-                    summary[point][param] = val
+                    data_1d_tmp.append(val)
+                    if add_columns:
+                        columns_1d.append(param)
+
+                # store time information, if requested
+                if len(var_vals) > len(variables) and timeseries:
+                    data_1d_tmp.append(var_vals[-1])
+                    if add_columns:
+                        columns_1d.append('time')
+
+                # store stability information
                 if stability:
-                    summary[point]['stability'] = self.get_stability(solution, s, point)
+                    data_1d_tmp.append(self.get_stability(solution, s, point))
+                    if add_columns:
+                        columns_1d.append('stability')
+
                 if period or lyapunov_exp or eigenvals:
-                    p = summary[point]['PAR(11)'] if 'PAR(11)' in params else self.get_params(s, ['PAR(11)'])[0]
-                    if period:
-                        summary[point]['period'] = p
-                    if eigenvals or lyapunov_exp:
-                        evs = self.get_eigenvalues(solution, branch, point)
-                        if eigenvals:
-                            summary[point]['eigenvalues'] = evs
-                        if lyapunov_exp:
-                            summary[point]['lyapunov_exponents'] = self.get_lyapunov_exponent(evs, p)
-
-        return DataFrame(summary).T
+                    p = self.get_params(s, ['PAR(11)'])[0]
+
+                # store information about oscillation periods
+                if period:
+                    data_1d_tmp.append(p)
+                    if add_columns:
+                        columns_1d.append('period')
+
+                # store information about local eigenvalues/lyapunov exponents
+                if eigenvals or lyapunov_exp:
+                    evs = self.get_eigenvalues(solution, branch, point)
+                    if eigenvals:
+                        for i, v in enumerate(evs):
+                            data_2d_tmp.append(evs)
+                            if add_columns:
+                                columns_2d.append(('eigenvalues', i))
+                    if lyapunov_exp:
+                        for i, p in enumerate(self.get_lyapunov_exponent(evs, p)):
+                            data_2d_tmp.append(p)
+                            if add_columns:
+                                columns_2d.append(('lyapunov_exponents', i))
+
+            data_2d.append(data_2d_tmp)
+            data_1d.append(data_1d_tmp)
+            add_columns = False
+
+        # arrange data into DataFrame
+        df = DataFrame(data=data_2d, columns=MultiIndex.from_tuples(columns_2d), index=indices)
+        df2 = DataFrame(data=data_1d, columns=columns_1d, index=indices)
+        for i, key in enumerate(columns_1d):
+            df[key] = df2.loc[:, key]
+        return df
 
     def _call_auto(self, starting_point: Union[str, int], origin: Union[Any, dict], **auto_kwargs) -> Any:
         if starting_point:

diff --git a/pycobi_tests/test_odesystem.py b/pycobi_tests/test_odesystem.py
@@ -44,15 +44,13 @@ def test_1_1_init(auto_dir):
 
     # initialize ODESystem from YAML file
     model = "model_templates.neural_mass_models.qif.qif"
-    ode2, t1 = ODESystem.from_yaml(model, init_cont=False, file_name='qif_eq2', func_name="qif2", auto_dir=auto_dir)
-    ode2.close_session()
-    clear(t1)
+    ode2 = ODESystem.from_yaml(model, init_cont=False, file_name='qif_eq2', func_name="qif2", auto_dir=auto_dir)
+    ode2.close_session(clear_files=True)
 
     # initialize ODESystem from YAML file with different parameters
-    ode3, t2 = ODESystem.from_yaml(model, init_cont=False, file_name='qif_eq3', func_name="qif3", auto_dir=auto_dir,
-                                   node_vars={'p/qif_op/eta': 2.0})
-    ode3.close_session()
-    clear(t2)
+    ode3 = ODESystem.from_yaml(model, init_cont=False, file_name='qif_eq3', func_name="qif3", auto_dir=auto_dir,
+                               node_vars={'p/qif_op/eta': 2.0})
+    ode3.close_session(clear_files=True)
 
     # these tests should pass
     assert isinstance(ode1, ODESystem)
@@ -72,16 +70,14 @@ def test_1_2_run(auto_dir):
 
     # initialize ODESystem from YAML file
     model = "model_templates.neural_mass_models.qif.qif"
-    ode2, t1 = ODESystem.from_yaml(model, init_cont=True, file_name='qif_eq2', func_name="qif2", auto_dir=auto_dir,
-                                   NPR=100, NMX=5000)
-    ode2.close_session()
-    clear(t1)
+    ode2 = ODESystem.from_yaml(model, init_cont=True, file_name='qif_eq2', func_name="qif2", auto_dir=auto_dir,
+                               NPR=100, NMX=5000)
+    ode2.close_session(clear_files=True)
 
     # initialize ODESystem from YAML file with different parameters
-    ode3, t2 = ODESystem.from_yaml(model, init_cont=True, file_name='qif_eq3', func_name="qif3", auto_dir=auto_dir,
-                                   node_vars={'p/qif_op/eta': 2.0}, NPR=100, NMX=5000)
-    ode3.close_session()
-    clear(t2)
+    ode3 = ODESystem.from_yaml(model, init_cont=True, file_name='qif_eq3', func_name="qif3", auto_dir=auto_dir,
+                               node_vars={'p/qif_op/eta': 2.0}, NPR=100, NMX=5000)
+    ode3.close_session(clear_files=True)
 
     # these tests should pass
     assert (ode1[0].loc[:, "U(1)"] - ode2[0].loc[:, "U(1)"]).sum()[0] == approx(0.0, rel=accuracy, abs=accuracy)