ntsa.bifurcation

At a glance

Function	One-liner
bifurcation_sweep(model, param, values, ...)	Parameter sweep collecting observable extrema: one ensemble forecast — member k at values[k] — by default, or serial branch-following with continuation=True.
plot_bifurcation(values, peaks, ...)	Bifurcation diagram from bifurcation_sweep output.

Produced by the module demo: python -m ntsa.bifurcation --model lorenz96 --save-figs.

Full reference

ntsa.bifurcation

Bifurcation sweeps and diagrams.

Run as a script for a bifurcation-diagram demo, e.g.::

python -m ntsa.bifurcation --model lorenz63 --param rho --pmin 20 --pmax 100 --save-figs

bifurcation_sweep(model, param, values, dt=None, t_transient=None, t_sample=None, extrema=('max',), continuation=False, seed=0, **overrides)

Sweep param over values, collecting local extrema of every observable.

By default the whole sweep is ONE ensemble forecast: every parameter value shares the same initial state, param is augmented into the state (init_ensemble(est_alpha=[param])) and member k integrates at values[k] — the integrator parallelizes across parameter values in a single run. With continuation=True it falls back to a serial m=1 loop where each run starts from the previous endpoint plus a small perturbation — the classic way to follow an attractor branch (e.g. through a hysteresis), and inherently sequential.

Returns:

Name	Type	Description
values	ndarray
peaks	dict	peaks[ext][iq][k] = extrema (ragged 1-D arrays) of observable iq at values[k], for each ext in extrema ('max' and/or 'min').

Source code in ntsa/bifurcation.py

def bifurcation_sweep(model, param, values, dt=None, t_transient=None, t_sample=None,
                      extrema=('max',), continuation=False, seed=0, **overrides):
    """Sweep `param` over `values`, collecting local extrema of every observable.

    By default the whole sweep is ONE ensemble forecast: every parameter value
    shares the same initial state, `param` is augmented into the state
    (``init_ensemble(est_alpha=[param])``) and member `k` integrates at
    ``values[k]`` — the integrator parallelizes across parameter values in a
    single run. With ``continuation=True`` it falls back to a serial m=1 loop
    where each run starts from the previous endpoint plus a small perturbation —
    the classic way to follow an attractor branch (e.g. through a hysteresis),
    and inherently sequential.

    Returns
    -------
    values : ndarray
    peaks : dict
        ``peaks[ext][iq][k]`` = extrema (ragged 1-D arrays) of observable `iq`
        at ``values[k]``, for each `ext` in `extrema` ('max' and/or 'min').
    """
    if param not in model.params:
        raise ValueError(f"'{param}' is not in {type(model).__name__}.params = {model.params}")
    rng = np.random.default_rng(seed)
    values = np.asarray(values, dtype=float)
    psi0 = model.psi0[:, 0] + 1e-3 * rng.standard_normal(model.Nphi)
    peaks = {ext: [[] for _ in range(model.Nq)] for ext in extrema}

    def collect(y_k, trim=False):
        i0 = stationary_start(y_k[:, 0]) if trim else 0
        for ext in extrema:
            sign = 1.0 if ext == 'max' else -1.0
            for iq in range(y_k.shape[1]):
                idx, _ = find_peaks(sign * y_k[i0:, iq])
                peaks[ext][iq].append(y_k[i0 + idx, iq])

    if continuation:
        for val in tqdm(values, desc=f'{param} sweep'):
            mi = respawn(model, psi0=psi0, dt=dt, **{param: val}, **overrides)
            _, y, psi = run_long(mi, t_run=t_sample or 20 * mi.t_CR, t_transient=t_transient)
            collect(y)   # run_long already trimmed the residual transient
            psi0 = psi[-1] + 1e-6 * rng.standard_normal(model.Nphi)
        return values, peaks

    # identical ICs: one ensemble run, member k carrying values[k] as its parameter
    mi = respawn(model, psi0=psi0, dt=dt, **overrides)
    aug = np.vstack([np.tile(psi0[:, None], (1, values.size)), values[None, :]])
    mi.init_ensemble(m=values.size, est_alpha=[param], ensemble_psi0=aug)
    if t_transient is None:
        t_transient = mi.t_transient
    Nt_tr = int(round(t_transient / mi.dt))
    if Nt_tr > 0:
        psi, _ = mi.time_integrate(Nt=Nt_tr)
        mi.update_history(psi[-1:], t=np.array([0.]), reset=True)
    psi, tt = mi.time_integrate(int(round((t_sample or 20 * mi.t_CR) / mi.dt)))
    mi.update_history(psi, tt)
    mi.close()
    y = mi.get_observable_hist()                     # (Nt, Nq, m)
    for k in range(values.size):
        collect(y[:, :, k], trim=True)               # per-member residual-transient trim
    return values, peaks

plot_bifurcation(values, peaks, param_label, obs_labels, filename=None)

Bifurcation diagram: one axis per observable, extrema vs parameter.

Source code in ntsa/bifurcation.py

def plot_bifurcation(values, peaks, param_label, obs_labels, filename=None):
    """Bifurcation diagram: one axis per observable, extrema vs parameter."""
    Nq = len(obs_labels)
    fig, axs = plt.subplots(Nq, 1, figsize=(8, 2.5 * Nq), sharex=True,
                            layout='constrained', squeeze=False)
    axs = axs.ravel()
    for iq, (ax, label) in enumerate(zip(axs, obs_labels)):
        for ext, per_obs in peaks.items():
            color = EXTREMA_COLORS.get(ext, 'k')
            for p, pk in zip(values, per_obs[iq]):
                ax.plot(np.full_like(pk, p), pk, '.', color=color, ms=1.5, alpha=0.4)
        ax.set_ylabel(f'local extrema of {label}')
    axs[-1].set_xlabel(param_label)
    if len(peaks) > 1:
        handles = [Line2D([], [], ls='', marker='.', color=EXTREMA_COLORS.get(e, 'k'), label=f'{e}ima')
                   for e in peaks]
        axs[0].legend(handles=handles, loc='best')
    if filename:
        fig.savefig(filename, dpi=200)
    return fig, axs