diff --git a/climatecritters/model_critters/__init__.py b/climatecritters/model_critters/__init__.py index 9e88a2e..8d709a1 100644 --- a/climatecritters/model_critters/__init__.py +++ b/climatecritters/model_critters/__init__.py @@ -16,4 +16,4 @@ from .stocker2003_bipolar_seesaw import * from .damped_spring import * from .pendulum import * -from .bistable_melcher import * +from .melcher25 import * diff --git a/climatecritters/model_critters/bistable_melcher.py b/climatecritters/model_critters/melcher25.py similarity index 93% rename from climatecritters/model_critters/bistable_melcher.py rename to climatecritters/model_critters/melcher25.py index e4e890f..e1c3651 100644 --- a/climatecritters/model_critters/bistable_melcher.py +++ b/climatecritters/model_critters/melcher25.py @@ -6,7 +6,7 @@ from ..core.ccmodel import CCModel -__all__ = ['BistableMelcherModel', 'classify_bistable_states'] +__all__ = ['Melcher25', 'classify_bistable_states'] def _classify_states(db, stadial_threshold, interstadial_threshold): @@ -38,7 +38,7 @@ def _classify_states(db, stadial_threshold, interstadial_threshold): return states -class BistableMelcherModel(CCModel): +class Melcher25(CCModel): """Two-equation bistable Itô SDE for stochastic Dansgaard-Oeschger transitions. Models the meridional buoyancy gradient Δb and buoyancy flux B as a coupled @@ -121,11 +121,11 @@ class BistableMelcherModel(CCModel): ```python import numpy as np import matplotlib.pyplot as plt - from climatecritters.model_critters.bistable_melcher import ( - BistableMelcherModel, classify_bistable_states + from climatecritters.model_critters.melcher25 import ( + Melcher25, classify_bistable_states ) - model = BistableMelcherModel(sigma=0.2, gamma=1.5, alpha=-0.4) + model = Melcher25(sigma=0.2, gamma=1.5, alpha=-0.4) output = model.integrate( t_span=(0, 599.88), y0=[1.0, 0.0], method='heun_maruyama', dt=0.012, @@ -135,7 +135,7 @@ class BistableMelcherModel(CCModel): # Δb time series as a pyleoclim Series — ready to plot or analyse ts = output.to_pyleo('db') ts.plot() - plt.savefig('docs/reference/figures/BistableMelcher_example.png', + plt.savefig('docs/reference/figures/Melcher25_example.png', dpi=150, bbox_inches='tight') # Stadial/interstadial classification is computed automatically @@ -153,12 +153,12 @@ class BistableMelcherModel(CCModel): ```python import numpy as np import matplotlib.pyplot as plt - from climatecritters.model_critters.bistable_melcher import BistableMelcherModel + from climatecritters.model_critters.melcher25 import Melcher25 t_arr = np.linspace(0, 599.88, 5000) gamma_ramp = np.linspace(0.8, 3.2, 5000) # γ increases over time - model = BistableMelcherModel( + model = Melcher25( sigma=0.2, gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)), alpha=0.0, @@ -170,7 +170,7 @@ class BistableMelcherModel(CCModel): ) ts = output.to_pyleo('db') ts.plot() - plt.savefig('docs/reference/figures/BistableMelcher_ramp_example.png', + plt.savefig('docs/reference/figures/Melcher25_ramp_example.png', dpi=150, bbox_inches='tight') ``` @@ -340,7 +340,7 @@ def compute_stability_thresholds(self, alpha): def classify_bistable_states(signal, alpha, b0=0.625, q0=-9.0, q1=12.0, tau=0.902): """Classify a Δb signal into stadial (1) / interstadial (0) states. - Applies the same hysteresis threshold logic that ``BistableMelcherModel`` + Applies the same hysteresis threshold logic that ``Melcher25`` uses internally after integration. Useful for reclassifying a Δb signal post-hoc — for example after adding proxy noise — without re-running the SDE. @@ -362,10 +362,10 @@ def classify_bistable_states(signal, alpha, b0=0.625, q0=-9.0, q1=12.0, tau=0.90 See also -------- - BistableMelcherModel.compute_stability_thresholds : Threshold computation. - BistableMelcherModel.populate_diagnostics_from_history : Auto-classification + Melcher25.compute_stability_thresholds : Threshold computation. + Melcher25.populate_diagnostics_from_history : Auto-classification after ``integrate()``. """ - model = BistableMelcherModel(b0=b0, q0=q0, q1=q1, tau=tau) + model = Melcher25(b0=b0, q0=q0, q1=q1, tau=tau) stadial_thresh, interstadial_thresh = model.compute_stability_thresholds(alpha) return _classify_states(np.asarray(signal, dtype=float), stadial_thresh, interstadial_thresh) diff --git a/climatecritters/tests/test_signal_models_bistable_melcher.py b/climatecritters/tests/test_signal_models_melcher25.py similarity index 82% rename from climatecritters/tests/test_signal_models_bistable_melcher.py rename to climatecritters/tests/test_signal_models_melcher25.py index 613e2e9..1239c73 100644 --- a/climatecritters/tests/test_signal_models_bistable_melcher.py +++ b/climatecritters/tests/test_signal_models_melcher25.py @@ -1,18 +1,18 @@ -"""Tests for climatecritters.model_critters.bistable_melcher.""" +"""Tests for climatecritters.model_critters.melcher25.""" import numpy as np import pytest -from climatecritters.model_critters.bistable_melcher import ( - BistableMelcherModel, classify_bistable_states, +from climatecritters.model_critters.melcher25 import ( + Melcher25, classify_bistable_states, ) -class TestSignalModelsBistableMelcherIntegrate: +class TestSignalModelsMelcher25Integrate: @pytest.mark.parametrize('y0', [[1.0, 0.0], [0.6, 0.0]]) @pytest.mark.parametrize('method', ['heun_maruyama', 'euler_maruyama', 'milstein']) def test_integrate_t0(self, y0, method): - model = BistableMelcherModel(sigma=0.2, gamma=1.5, alpha=-0.4) + model = Melcher25(sigma=0.2, gamma=1.5, alpha=-0.4) output = model.integrate( t_span=(0, 12), y0=y0, method=method, dt=0.012, kwargs={'random_seed': 0, 'si': 0.12}, @@ -25,15 +25,15 @@ def test_integrate_t0(self, y0, method): def test_integrate_with_deterministic_method_t1(self): """uses_post_history=True models should also work with non-SDE methods.""" - model = BistableMelcherModel(alpha=-0.4) + model = Melcher25(alpha=-0.4) output = model.integrate(t_span=(0, 1.2), y0=[1.0, 0.0], method='euler', dt=0.012) assert np.all(np.isfinite(output.state_variables['db'])) assert model.stadial_threshold is not None -class TestSignalModelsBistableMelcherThresholds: +class TestSignalModelsMelcher25Thresholds: def test_thresholds_set_after_integrate_t0(self): - model = BistableMelcherModel(alpha=-0.4) + model = Melcher25(alpha=-0.4) model.integrate( t_span=(0, 12), y0=[1.0, 0.0], method='heun_maruyama', dt=0.012, kwargs={'random_seed': 1, 'si': 0.12}, @@ -43,7 +43,7 @@ def test_thresholds_set_after_integrate_t0(self): assert model.stadial_threshold < model.interstadial_threshold def test_thresholds_match_compute_stability_thresholds_t1(self): - model = BistableMelcherModel(alpha=-0.4) + model = Melcher25(alpha=-0.4) model.integrate( t_span=(0, 12), y0=[1.0, 0.0], method='heun_maruyama', dt=0.012, kwargs={'random_seed': 1, 'si': 0.12}, @@ -57,8 +57,8 @@ def test_callable_alpha_thresholds_match_scalar_t2(self): constant-valued callable alpha to the same thresholds as passing the same constant directly (exercises the branch reworked alongside compute_stability_thresholds's own mean/float reduction).""" - const_model = BistableMelcherModel(alpha=-0.4) - tv_model = BistableMelcherModel(alpha=lambda t: -0.4) + const_model = Melcher25(alpha=-0.4) + tv_model = Melcher25(alpha=lambda t: -0.4) t_span, y0, dt = (0, 12), [1.0, 0.0], 0.012 @@ -72,7 +72,7 @@ def test_callable_alpha_thresholds_match_scalar_t2(self): def test_array_alpha_thresholds_use_mean_t3(self): """compute_stability_thresholds should reduce an array-like alpha to its mean.""" - model = BistableMelcherModel() + model = Melcher25() alpha_arr = np.array([-0.2, -0.6]) stadial, interstadial = model.compute_stability_thresholds(alpha_arr) expected_stadial, expected_interstadial = model.compute_stability_thresholds(np.mean(alpha_arr)) @@ -80,9 +80,9 @@ def test_array_alpha_thresholds_use_mean_t3(self): assert interstadial == expected_interstadial -class TestSignalModelsBistableMelcherClassifyStandalone: +class TestSignalModelsMelcher25ClassifyStandalone: def test_classify_bistable_states_matches_model_t0(self): - model = BistableMelcherModel(alpha=-0.4) + model = Melcher25(alpha=-0.4) output = model.integrate( t_span=(0, 12), y0=[1.0, 0.0], method='heun_maruyama', dt=0.012, kwargs={'random_seed': 3, 'si': 0.12}, @@ -96,7 +96,7 @@ def test_classify_bistable_states_hysteresis_t1(self): """A signal that dips below the stadial threshold and rises above the interstadial threshold should flip states with hysteresis (no chatter for values between the two thresholds).""" - model = BistableMelcherModel(alpha=-0.4) + model = Melcher25(alpha=-0.4) stadial, interstadial = model.compute_stability_thresholds(-0.4) mid = 0.5 * (stadial + interstadial) signal = np.array([interstadial + 0.1, mid, stadial - 0.1, mid, interstadial + 0.1]) @@ -104,12 +104,12 @@ def test_classify_bistable_states_hysteresis_t1(self): assert list(states) == [0, 0, 1, 1, 0] -class TestSignalModelsBistableMelcherSDENoise: +class TestSignalModelsMelcher25SDENoise: def test_zero_sigma_is_deterministic_t0(self): """sigma=0 should make euler_maruyama reduce to the deterministic drift, independent of the random seed.""" - model_a = BistableMelcherModel(sigma=0.0, alpha=-0.4) - model_b = BistableMelcherModel(sigma=0.0, alpha=-0.4) + model_a = Melcher25(sigma=0.0, alpha=-0.4) + model_b = Melcher25(sigma=0.0, alpha=-0.4) t_span, y0, dt = (0, 12), [1.0, 0.0], 0.012 out_a = model_a.integrate(t_span=t_span, y0=y0, method='euler_maruyama', dt=dt, kwargs={'random_seed': 1, 'si': 0.12}) @@ -118,16 +118,16 @@ def test_zero_sigma_is_deterministic_t0(self): assert np.allclose(out_a.state_variables['db'], out_b.state_variables['db']) def test_sde_noise_shape_and_scale_t1(self): - model = BistableMelcherModel(sigma=0.3) + model = Melcher25(sigma=0.3) diffusion = model.sde_noise(0.0, [1.0, 0.0]) assert diffusion.shape == (2,) assert np.allclose(diffusion, 0.3) -class TestSignalModelsBistableMelcherTimeVaryingParams: +class TestSignalModelsMelcher25TimeVaryingParams: def test_time_varying_params_match_constants_t0(self): - model_const = BistableMelcherModel(gamma=1.5, alpha=-0.4) - model_tv = BistableMelcherModel( + model_const = Melcher25(gamma=1.5, alpha=-0.4) + model_tv = Melcher25( gamma=lambda t: 1.5, alpha=lambda t, x: -0.4, ) diff --git a/docs/_quarto.yml b/docs/_quarto.yml index d77072d..788cd6b 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -126,7 +126,7 @@ website: - text: "Thermohaline Circulation (Stommel, 1961)" href: notebooks/model_demos/stommel.ipynb - text: "Thermohaline Circulation (Melcher et al 2025)" - href: notebooks/model_demos/bistable_melcher.ipynb + href: notebooks/model_demos/melcher25.ipynb - section: Lorenz & Rössler contents: - text: "Roessler" @@ -242,7 +242,7 @@ quartodoc: - Stocker2003BipolarSeesaw - Stocker2003ExtendedSeaIceSeesaw - Stommel - - BistableMelcherModel + - Melcher25 - title: Utils diff --git a/notebooks/model_demos/bistable_melcher.ipynb b/notebooks/model_demos/melcher25.ipynb similarity index 99% rename from notebooks/model_demos/bistable_melcher.ipynb rename to notebooks/model_demos/melcher25.ipynb index e527c84..ef8c164 100644 --- a/notebooks/model_demos/bistable_melcher.ipynb +++ b/notebooks/model_demos/melcher25.ipynb @@ -5,7 +5,7 @@ "id": "cell-title", "metadata": {}, "source": [ - "# BistableMelcherModel — Stochastic Dansgaard-Oeschger Transitions (Melcher et al. 2025)\n", + "# Melcher25 — Stochastic Dansgaard-Oeschger Transitions (Melcher et al. 2025)\n", "\n", "**Author:** Maryam Niati \n", "**Date:** last-modified\n", @@ -14,9 +14,9 @@ "\n", "_Abstract_\n", "\n", - "> `BistableMelcherModel` implements the two-equation bistable Itô SDE from Melcher et al. (2025) that reproduces the statistical properties of Dansgaard-Oeschger (DO) events recorded in the NGRIP ice core. The system switches stochastically between a cold stadial state (weak AMOC) and a warm interstadial state (strong AMOC). This notebook shows how to generate synthetic DO records, automatically detect tipping-point transitions, visualise the stadial/interstadial classification, and explore the effect of noise amplitude and time-varying CO₂ forcing.\n", + "> `Melcher25` implements the two-equation bistable Itô SDE from Melcher et al. (2025) that reproduces the statistical properties of Dansgaard-Oeschger (DO) events recorded in the NGRIP ice core. The system switches stochastically between a cold stadial state (weak AMOC) and a warm interstadial state (strong AMOC). This notebook shows how to generate synthetic DO records, automatically detect tipping-point transitions, visualise the stadial/interstadial classification, and explore the effect of noise amplitude and time-varying CO₂ forcing.\n", "\n", - "**Keywords:** BistableMelcherModel, Dansgaard-Oeschger, bistable, tipping points, stadial, interstadial, AMOC, stochastic SDE, synthetic palaeoclimate, Heun-Maruyama" + "**Keywords:** Melcher25, Dansgaard-Oeschger, bistable, tipping points, stadial, interstadial, AMOC, stochastic SDE, synthetic palaeoclimate, Heun-Maruyama" ] }, { @@ -26,7 +26,7 @@ "source": [ "## Overview\n", "\n", - "`BistableMelcherModel` models the meridional buoyancy gradient $\\Delta b$ and buoyancy flux $B$ as a coupled stochastic system:\n", + "`Melcher25` models the meridional buoyancy gradient $\\Delta b$ and buoyancy flux $B$ as a coupled stochastic system:\n", "\n", "$$\\frac{d(\\Delta b)}{dt} = -B - \\left|q_0 + q_1(\\Delta b - b_0)\\right|(\\Delta b - b_0) + \\sigma\\,dW$$\n", "\n", @@ -63,8 +63,8 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.patches as mpatches\n", "\n", - "from climatecritters.model_critters.bistable_melcher import (\n", - " BistableMelcherModel, classify_bistable_states\n", + "from climatecritters.model_critters.melcher25 import (\n", + " Melcher25, classify_bistable_states\n", ")" ] }, @@ -93,7 +93,7 @@ } ], "source": [ - "model = BistableMelcherModel(sigma=0.2, gamma=1.5, alpha=-0.4)\n", + "model = Melcher25(sigma=0.2, gamma=1.5, alpha=-0.4)\n", "\n", "output = model.integrate(\n", " t_span=(0, 599.88), y0=[1.0, 0.0],\n", @@ -223,7 +223,7 @@ "axes[1].set_ylabel('B')\n", "axes[1].set_xlabel('time (kyr)')\n", "\n", - "fig.suptitle('BistableMelcherModel — state variables')\n", + "fig.suptitle('Melcher25 — state variables')\n", "plt.tight_layout()\n", "plt.show()" ] @@ -255,7 +255,7 @@ "output_type": "stream", "text": [ "\n", - "Parameters — BistableMelcherModel\n", + "Parameters — Melcher25\n", "══════════════════════════════════════════════════════════════════════════\n", "Name Current value Description\n", "──────────────────────────────────────────────────────────────────────────\n", @@ -335,7 +335,7 @@ "sigma_outputs = {}\n", "\n", "for s in sigmas:\n", - " m = BistableMelcherModel(sigma=s, gamma=1.5, alpha=-0.4)\n", + " m = Melcher25(sigma=s, gamma=1.5, alpha=-0.4)\n", " out = m.integrate(\n", " t_span=(0, 599.88), y0=[1.0, 0.0],\n", " method='heun_maruyama', dt=0.012,\n",