.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/automated_continuation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_automated_continuation.py: Automated Codim-2 Search and Period-Doubling Cascades ===================================================== `PyCoBi` ships two convenience helpers for the bookkeeping-heavy parts of a bifurcation analysis: * :func:`codim2_search` — given a list of codim-1 starting points (folds, Hopfs, period-doublings) on an existing continuation, runs a 2-parameter continuation of each, walks the resulting curves for codim-2 points, and recursively continues the codim-1 bifurcations that emerge from them. Supports recursive handling of zero-Hopf (ZH), generalised-Hopf (GH / Bautin), and Bogdanov-Takens (BT) points. * :func:`continue_period_doubling_bf` — chases a cascade of period-doubling bifurcations in 2 parameters, recursing on every new PD encountered. Useful for tracing the boundaries of period-doubling routes to chaos. This example uses the QIF mean-field model with bi-exponential spike-frequency adaptation (``qif_biexp_sfa.yaml`` next to this script). The bi-exponential kernel is a strict generalisation of the alpha-kernel QIF-SFA from :ref:`Hopf Bifurcation and Limit Cycle Continuation`: both kernels satisfy :math:`\tau_a^2 A'' + 2 \tau_a A' + A = \alpha r \tau_a` when ``tau_r == tau_d == tau_a``, so the codim-2 structure in the :math:`(\bar\eta,\, \Delta)` plane (generalised-Hopf, Bogdanov-Takens, cusps) is identical to the alpha-kernel case at that parameter point — verified by side-by-side bifurcation analysis. Picking the bi-exponential model now means we can also explore the period-doubling regime that opens up when ``tau_r`` is taken much smaller than ``tau_d`` (Section 2.1 below). References ^^^^^^^^^^ .. [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. .. GENERATED FROM PYTHON SOURCE LINES 40-47 Step 1: Load the model ^^^^^^^^^^^^^^^^^^^^^^ Bi-exponential QIF-SFA from the co-located YAML, with rise and decay adaptation time constants both set to 10 — the alpha-kernel-equivalent parameter point. Coupling ``J = 15 sqrt(2)``, ``alpha = 1.0``, ``Delta = 2.0`` to match the QIF-SFA example. .. GENERATED FROM PYTHON SOURCE LINES 47-83 .. code-block:: Python from pathlib import Path import numpy as np import matplotlib.pyplot as plt from pycobi import ODESystem from pycobi.automated_continuation import ( codim2_search, continue_period_doubling_bf, ) here = Path(__file__).resolve().parent yaml_path = str(here / 'qif_biexp_sfa' / 'qif_biexp_sfa') ode = ODESystem.from_yaml( yaml_path, auto_dir="~/PycharmProjects/auto-07p", node_vars={ 'p/qif_biexp_sfa_op/Delta': 2.0, 'p/qif_biexp_sfa_op/alpha': 0.8, # Start at tau_r = 11 (slightly above tau_d=10) and continue down in # tau_r before any eta work — see Step 2 below. 'p/qif_biexp_sfa_op/tau_r': 11.0, 'p/qif_biexp_sfa_op/tau_d': 10.0, 'p/qif_biexp_sfa_op/eta': -8.0, }, edge_vars=[('p/qif_biexp_sfa_op/r', 'p/qif_biexp_sfa_op/r_in', {'weight': 15.0 * np.sqrt(2.0)})], init_cont=True, NPR=100, NMX=30000, ) # Auto-07p doesn't ship a built-in style for the cusp (CP) marker; add one # so it shows up in the codim-2 diagram below alongside BT and GH. ode.update_bifurcation_style('CP', marker='d', color='#7F4FBF') .. GENERATED FROM PYTHON SOURCE LINES 84-99 Step 2: Pre-scan in :math:`\\tau_r` to anchor the two regimes ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We'll cover two parameter regimes of this model in one example: 1. :math:`\\tau_r = \\tau_d = 10` (alpha-kernel-equivalent): clean codim-2 structure in :math:`(\\bar\\eta,\\, \\Delta)` — BT / GH / CP, no PD. 2. :math:`\\tau_r = 0.1 \\ll \\tau_d`: same codim-2 backbone plus a period-doubling cascade on the LC born at the upper Hopf. Rather than re-instantiate the model twice, we continue the steady state in :math:`\\tau_r` from the initial :math:`\\tau_r = 11` down to small values, planting user points at :math:`\\tau_r = \\tau_d = 10` (UZ1) and at :math:`\\tau_r = 0.1` (UZ2). Each UZ then serves as the starting point for an :math:`\\bar\\eta` continuation in its own regime. .. GENERATED FROM PYTHON SOURCE LINES 99-110 .. code-block:: Python tau_sols, tau_cont = ode.run( starting_point='EP2', name='tau_r_branch', ICP='p/qif_biexp_sfa_op/tau_r', IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4, NMX=5000, NPR=100, DS=-1e-3, DSMIN=1e-9, DSMAX=5e-2, UZR={'p/qif_biexp_sfa_op/tau_r': [10.0, 0.1]}, UZSTOP={'p/qif_biexp_sfa_op/tau_r': [0.01, 10.0]}, ) print(f"tau_r pre-scan bifurcations: {dict(tau_sols['bifurcation'].value_counts())}") .. GENERATED FROM PYTHON SOURCE LINES 111-122 Shared bookkeeping ^^^^^^^^^^^^^^^^^^ We run the same analysis pipeline at two values of :math:`\\tau_r` (Steps 3-4 and 5-6 below), so it pays to factor the codim-1 / codim-2 bookkeeping into a small helper. ``codim1_points`` lists the four starting points we expect on every :math:`\\bar\\eta` scan (the model has two folds and two Hopfs in the parameter regime we're exploring); ``CURVE_COLORS`` keys colours by bifurcation type so the codim-2 diagrams collapse their legends to two entries even though we draw eight curves. .. GENERATED FROM PYTHON SOURCE LINES 122-260 .. code-block:: Python codim1_points = [ ('LP1', 'fold'), ('LP2', 'fold'), ('HB1', 'Hopf'), ('HB2', 'Hopf'), ] CURVE_COLORS = {'fold': '#1F77B4', 'Hopf': '#FF7F0E'} PD_COLOR = '#D62728' def run_codim2_in_delta(origin_name, name_prefix): """Run the 8-call codim-2 search (LP1/LP2/HB1/HB2 × DS=±1e-3) in :math:`(\\bar\\eta,\\, \\Delta)` and return ``[(key, bif_type)]``. Forward + reverse continuations are kept separate (rather than one bidirectional call) because bidirectional gets stuck bouncing at the cusp adjacent to LP1 / LP2. """ shared_kwargs = dict( pyauto_instance=ode, params=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'], origin=origin_name, max_recursion_depth=0, NMX=1500, NPR=10, DSMIN=1e-9, DSMAX=5e-2, RL0=-15.0, RL1=5.0, bidirectional=False, UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.0, 4.0]}, ) curves = [] for sp, bif_type in codim1_points: for ds in (1e-3, -1e-3): try: result = codim2_search( starting_points=[sp], DS=ds, name=f'{name_prefix}_{sp}_{"pos" if ds > 0 else "neg"}', **shared_kwargs, ) curves.append((list(result.keys())[0], bif_type)) except Exception as exc: print(f"{name_prefix} {sp} DS={ds:+g}: skipped " f"({type(exc).__name__}: {exc})") return curves def plot_2d_diagram(ax, codim2_curves, eta_origin, *, pd_names=(), title, xlim=(-6.0, 2.0), ylim=(0.0, 2.5)): """Plot a 2D codim-2 diagram in :math:`(\\bar\\eta,\\, \\Delta)` from the curves produced by :func:`run_codim2_in_delta`. Optionally overlay one or more PD continuations (rendered in red). Fold curves all share the blue colour, Hopf curves all share orange — so the legend has at most three entries (fold curve, Hopf curve, PD curve).""" labels_used = set() for key, bif_type in codim2_curves: color = CURVE_COLORS[bif_type] label = f'{bif_type} curve' if bif_type not in labels_used else None labels_used.add(bif_type) ode.plot_continuation( 'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta', cont=key, ax=ax, line_color_stable=color, line_color_unstable=color, line_style_stable='solid', line_style_unstable='solid', bifurcation_legend=False, get_stability=False, ignore=['LP', 'HB', 'UZ'], label=label, ) for pd_name in pd_names: label = 'PD curve' if 'PD' not in labels_used else None labels_used.add('PD') ode.plot_continuation( 'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta', cont=pd_name, ax=ax, line_color_stable=PD_COLOR, line_color_unstable=PD_COLOR, line_style_stable='solid', line_style_unstable='solid', bifurcation_legend=False, get_stability=False, ignore=['PD', 'UZ'], label=label, ) for sp, bif_type in codim1_points: try: sol, _, _ = ode.get_solution(point=sp, cont=eta_origin) eta_p = float(sol['eta']) delta_p = float(sol['Delta']) ax.scatter(eta_p, delta_p, marker='*', s=120, c=CURVE_COLORS[bif_type], edgecolor='k', linewidth=0.5, zorder=10) ax.annotate(sp, (eta_p, delta_p), xytext=(5, 5), textcoords='offset points', fontsize=9) except Exception: pass ax.set_xlim(*xlim) ax.set_ylim(*ylim) ax.set_xlabel(r'$\bar\eta$') ax.set_ylabel(r'$\Delta$') ax.set_title(title) ax.legend(loc='best') def plot_1d_diagram(ax, eta_cont_name, lc_cont_name, *, title, xlim=(-8.0, 2.0), cascade_names=()): """Plot a 1D bifurcation diagram in :math:`(\\bar\\eta,\\, r)` with the equilibrium branch (blue), the limit cycle branch (orange), and optionally a list of cascade LC continuations (red — each doubled cycle from a PD branch switch) overlaid. Codim-1 markers (HB / LP / PD) appear automatically; UZ / BP / EP are filtered out as cosmetic noise.""" ode.plot_continuation( 'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r', cont=eta_cont_name, ax=ax, line_color_stable='#1F77B4', line_color_unstable='#1F77B4', bifurcation_legend=False, label='equilibrium branch', ignore=['UZ', 'BP', 'EP'], ) ode.plot_continuation( 'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r', cont=lc_cont_name, ax=ax, line_color_stable='#FF7F0E', line_color_unstable='#FF7F0E', bifurcation_legend=False, label='limit cycle branch', ignore=['UZ', 'BP', 'EP'], ) cascade_label_done = False for cascade_name in cascade_names: label = None if cascade_label_done else 'PD cascade LCs' cascade_label_done = True try: ode.plot_continuation( 'p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/r', cont=cascade_name, ax=ax, line_color_stable=PD_COLOR, line_color_unstable=PD_COLOR, bifurcation_legend=False, label=label, ignore=['UZ', 'BP', 'EP'], ) except Exception as exc: print(f" warning: could not plot cascade LC {cascade_name!r}: {exc}") ax.set_xlim(*xlim) ax.set_xlabel(r'$\bar\eta$') ax.set_ylabel(r'$r$') ax.set_title(title) ax.legend(loc='best') .. GENERATED FROM PYTHON SOURCE LINES 261-269 Step 3: Continuation + codim-2 search at :math:`\\tau_r = \\tau_d = 10` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This is the alpha-kernel-equivalent regime — clean codim-2 structure (BT / GH / CP), no period doublings. A bidirectional :math:`\\bar\\eta` scan from ``UZ1`` exposes the codim-1 starting points; we then branch-switch to the limit cycle born at ``HB2`` and run :func:`codim2_search` from each codim-1 point. .. GENERATED FROM PYTHON SOURCE LINES 269-295 .. code-block:: Python eta_sols, eta_cont = ode.run( origin='tau_r_branch', starting_point='UZ1', name='eta_branch', ICP='p/qif_biexp_sfa_op/eta', bidirectional=True, RL0=-8.0, RL1=2.0, IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4, NMX=2000, NPR=10, DS=1e-4, DSMIN=1e-8, DSMAX=5e-2, ITMX=40, ITNW=40, NWTN=12, ) print("\neta scan (tau_r=10):", dict(eta_sols['bifurcation'].value_counts())) lc_sols, lc_cont = ode.run( origin='eta_branch', starting_point='HB2', name='lc_branch', IPS=2, ISP=2, ISW=-1, ICP=['p/qif_biexp_sfa_op/eta', 11], NMX=2000, NPR=10, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2, bidirectional=True, get_period=True, STOP=["BP3", "LP5"] ) print("LC (tau_r=10):", dict(lc_sols['bifurcation'].value_counts())) codim2_curves = run_codim2_in_delta('eta_branch', name_prefix='delta10') print("codim-2 curves at tau_r=10:") for key, bif_type in codim2_curves: print(f" {key} ({bif_type}): " f"{dict(ode.get_summary(key)['bifurcation'].value_counts())}") .. GENERATED FROM PYTHON SOURCE LINES 296-303 Step 4: Figures at :math:`\\tau_r = 10` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Two figures: a 1D bifurcation diagram in :math:`(\\bar\\eta,\\, r)` with the equilibrium and limit-cycle branches overlaid, and a 2D codim-2 diagram in :math:`(\\bar\\eta,\\, \\Delta)` showing the BT / CP / GH structure. .. GENERATED FROM PYTHON SOURCE LINES 303-316 .. code-block:: Python fig, ax = plt.subplots(figsize=(7, 4.5)) plot_1d_diagram(ax, 'eta_branch', 'lc_branch', title=r'1D bifurcation diagram, $\tau_r = 10$') plt.tight_layout() plt.show() fig, ax = plt.subplots(figsize=(7, 5)) plot_2d_diagram(ax, codim2_curves, 'eta_branch', title=r'codim-2 bifurcation diagram, $\tau_r = 10$') plt.tight_layout() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 317-331 Reading the 2D diagram: * Blue fold curves traced from ``LP1`` / ``LP2``, both directions each. The cusp ``CP`` where these collide marks the closing of the bistable wedge. * Orange Hopf curves traced from ``HB1`` / ``HB2``. The ``BT`` (Bogdanov-Takens) point is where the Hopf curve from ``HB1`` meets the fold manifold — two fold branches pass through BT and one Hopf branch terminates at it (a homoclinic curve also emerges from BT, but tracking it requires auto-07p's HomCont package ``IPS=9`` and is left as a manual follow-up). * The ``GH`` (generalised-Hopf / Bautin) point on the Hopf curve from ``HB2`` is where the first Lyapunov coefficient changes sign — supercritical Hopfs on one side, subcritical on the other. .. GENERATED FROM PYTHON SOURCE LINES 333-340 Step 5: Continuation + codim-2 + PD cascade at :math:`\\tau_r = 0.1` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Same pipeline at the second user point ``UZ2``, with the addition of a :func:`continue_period_doubling_bf` call after the LC continuation — at :math:`\\tau_r = 0.1 \\ll \\tau_d` the LC born at ``HB2`` carries at least one period-doubling point. .. GENERATED FROM PYTHON SOURCE LINES 340-401 .. code-block:: Python eta_lo_sols, eta_lo_cont = ode.run( origin='tau_r_branch', starting_point='UZ2', name='eta_branch_lo', ICP='p/qif_biexp_sfa_op/eta', bidirectional=True, RL0=-8.0, RL1=2.0, IPS=1, ILP=1, ISP=2, ISW=1, NTST=400, NCOL=4, NMX=2000, NPR=10, DS=1e-4, DSMIN=1e-8, DSMAX=5e-2, ITMX=40, ITNW=40, NWTN=12, ) print("\neta scan (tau_r=0.1):", dict(eta_lo_sols['bifurcation'].value_counts())) lc_lo_sols, lc_lo_cont = ode.run( origin='eta_branch_lo', starting_point='HB2', name='lc_branch_lo', IPS=2, ISP=2, ISW=-1, ICP=['p/qif_biexp_sfa_op/eta', 11], NMX=2000, NPR=20, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2, bidirectional=True, get_period=True, STOP=["BP3", "LP5"] ) print("LC (tau_r=0.1):", dict(lc_lo_sols['bifurcation'].value_counts())) codim2_curves_lo = run_codim2_in_delta('eta_branch_lo', name_prefix='delta01') print("codim-2 curves at tau_r=0.1:") for key, bif_type in codim2_curves_lo: print(f" {key} ({bif_type}): " f"{dict(ode.get_summary(key)['bifurcation'].value_counts())}") # Trace the locus of the *first PD point* on the LC in (eta, Delta), the # same way we'd trace any other codim-1 manifold — `codim2_search` treats # PD as a codim-1 starting point when `periodic=True` (so ICP gets PAR(11) # appended for the period). This gives the red PD curve in the 2D diagram. pd_locus = codim2_search( pyauto_instance=ode, params=['p/qif_biexp_sfa_op/eta', 'p/qif_biexp_sfa_op/Delta'], starting_points=['PD1'], origin='lc_branch_lo', max_recursion_depth=0, periodic=True, NMX=1500, NPR=20, DSMIN=1e-9, DSMAX=5e-2, bidirectional=True, UZSTOP={'p/qif_biexp_sfa_op/Delta': [0.0, 4.0]}, name='pd_locus', ) print("PD-locus continuations:", list(pd_locus.keys())) # Chase the period-doubling cascade itself with the dedicated helper. # Unlike `codim2_search` (which traces a locus in 2 params), this one # does the auto-07p "branch-switch onto each new doubled LC and repeat" # workflow — same recipe as the c.lor.2 / c.lor.3 manual cascade in # auto-07p's lor demo. Each cascade step is a separate 1-parameter LC # continuation in `eta` (with PAR(11) for the doubled period). cascade_names, _ = continue_period_doubling_bf( solution=ode.results[ode.get_continuation('lc_branch_lo').key], continuation=lc_lo_cont, pyauto_instance=ode, icp='p/qif_biexp_sfa_op/eta', max_iter=4, NMX=2000, NPR=20, DS=1e-3, DSMIN=1e-9, DSMAX=5e-2, STOP=["BP3", "LP5"], ) print("PD cascade LCs:", cascade_names) .. GENERATED FROM PYTHON SOURCE LINES 402-409 Step 6: Figures at :math:`\\tau_r = 0.1` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Same shape as Step 4. The 1D diagram now overlays each LC born along the PD cascade (each doubled cycle is its own continuation in eta). The 2D diagram overlays the PD locus from :func:`codim2_search` on the usual codim-2 backbone. .. GENERATED FROM PYTHON SOURCE LINES 409-425 .. code-block:: Python fig, ax = plt.subplots(figsize=(7, 4.5)) plot_1d_diagram(ax, 'eta_branch_lo', 'lc_branch_lo', title=r'1D bifurcation diagram, $\tau_r = 0.1$', cascade_names=cascade_names) plt.tight_layout() plt.show() fig, ax = plt.subplots(figsize=(7, 5)) plot_2d_diagram(ax, codim2_curves_lo, 'eta_branch_lo', pd_names=list(pd_locus.keys()), title=r'codim-2 + PD-locus bifurcation diagram, ' r'$\tau_r = 0.1$') plt.tight_layout() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 426-437 Step 7: Failure modes and what to expect on unfamiliar models ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Both helpers wrap a series of nested ``ODESystem.run`` calls. If any individual sub-run raises inside auto-07p (common on low-quality parameter regimes or aggressive step sizes), the failure surfaces as a :class:`UserWarning` and the search continues with the remaining starting points rather than aborting. Read the warnings carefully: they cite the auto-07p exception type and message, the sub-run label that failed, and the kwargs hook (``kwargs_1D_lc_cont``, ``kwargs_2D_cont``, ...) you can use to override the default constants for that path. .. GENERATED FROM PYTHON SOURCE LINES 439-441 Step 8: Clean up ^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 441-443 .. code-block:: Python ode.close_session(clear_files=True) .. _sphx_glr_download_auto_examples_automated_continuation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: automated_continuation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: automated_continuation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: automated_continuation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_