From b74e3a0cbddc9f1caae01af504c1ef6e140857fc Mon Sep 17 00:00:00 2001 From: Richert Date: Mon, 17 Oct 2022 12:09:40 -0500 Subject: [PATCH] added new use example, updated some of the plotting functions --- CHANGELOG.rst | 8 ++ documentation/doc/source/index.rst | 7 ++ documentation/doc/source/simple_example.rst | 9 +- documentation/doc/source/use_examples.rst | 6 ++ documentation/use_examples/qif_sfa.py | 59 ++++++++--- pycobi/pycobi.py | 111 ++++++++++++++------ pycobi_tests/test_odesystem.py | 26 ++--- 7 files changed, 163 insertions(+), 63 deletions(-) create mode 100644 documentation/doc/source/use_examples.rst diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d1fd431..5c76132 100644 --- 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 index 94b3d47..1c34dd3 100755 --- 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 index 58b2ded..a3dfb57 100644 --- 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 `_, which allows us to generate the Fortran files +To this end, we make use of the code generation software `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 new file mode 100644 index 0000000..c4b43c2 --- /dev/null +++ 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 index f1ec461..4cb8d4e 100644 --- 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 `_ 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 `_) # and grey triangles (`fold bifurcations `_), # 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 index aa3689b..fce4160 100644 --- 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 index 3b50eac..eed4456 100644 --- 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)