BistableMelcher -->Melcher25

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 *

climatecritters/model_critters/bistable_melcher.py → 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)

climatecritters/tests/test_signal_models_bistable_melcher.py → 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,17 +72,17 @@ 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))
         assert stadial == expected_stadial
         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,20 +96,20 @@ 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])
         states = classify_bistable_states(signal, alpha=-0.4)
         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,
         )

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
 

notebooks/model_demos/bistable_melcher.ipynb → 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",