Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 58 additions & 75 deletions climatecritters/model_critters/bistable_melcher.py
Original file line number Diff line number Diff line change
Expand Up @@ -116,84 +116,67 @@ class BistableMelcherModel(CCModel):

Examples
--------
**Generate a synthetic DO record and inspect states**

.. code-block:: python

import numpy as np
from climatecritters.model_critters.bistable_melcher import (
BistableMelcherModel, classify_bistable_states
)

model = BistableMelcherModel(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,
kwargs={'random_seed': 42, 'si': 0.12},
)

# Δb time series as a pyleoclim Series — ready to plot or analyse
ts = output.to_pyleo('db')
ts.plot()

# Stadial/interstadial classification is computed automatically
states = output.diagnostic_variables['states'] # 1 = cold, 0 = warm
print('stadial threshold:', model.stadial_threshold)
print('interstadial threshold:', model.interstadial_threshold)

# Count DO events (interstadial onsets)
n_events = int(np.sum(np.diff(states) == -1))
print('number of DO events:', n_events)

**Inspect all parameters interactively**

.. code-block:: python

model.doc('parameters') # pretty-prints names, current values, descriptions

**Time-varying forcing** (e.g. a slow CO2 ramp across a glacial cycle)

.. code-block:: python

import numpy as np
from climatecritters.model_critters.bistable_melcher import BistableMelcherModel

t_arr = np.linspace(0, 599.88, 5000)
gamma_ramp = np.linspace(0.8, 3.2, 5000) # γ increases over time

model = BistableMelcherModel(
sigma=0.2,
gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)),
alpha=0.0,
)
output = model.integrate(
t_span=(t_arr[0], t_arr[-1]), y0=[1.0, 0.0],
method='heun_maruyama', dt=0.012,
kwargs={'random_seed': 0, 'si': 0.12},
)
ts = output.to_pyleo('db')

**Reclassify after adding proxy noise**

.. code-block:: python

import numpy as np
from climatecritters.model_critters.bistable_melcher import classify_bistable_states

noisy_db = ts.value + np.random.default_rng(7).normal(0, 0.1, len(ts.value))
states_noisy = classify_bistable_states(noisy_db, alpha=-0.4)

**Sensitivity to calibration constants**

.. code-block:: python

model = BistableMelcherModel(b0=0.6, q0=-8.5, q1=11.0, tau=0.95)
model.doc('parameters') # verify the new values are registered
Generate a synthetic DO record and inspect states

```python
import numpy as np
import matplotlib.pyplot as plt
from climatecritters.model_critters.bistable_melcher import (
BistableMelcherModel, classify_bistable_states
)

model = BistableMelcherModel(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,
kwargs={'random_seed': 42, 'si': 0.12},
)

# Δ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',
dpi=150, bbox_inches='tight')

# Stadial/interstadial classification is computed automatically
states = output.diagnostic_variables['states'] # 1 = cold, 0 = warm
print('stadial threshold:', model.stadial_threshold)
print('interstadial threshold:', model.interstadial_threshold)

# Count DO events (interstadial onsets)
n_events = int(np.sum(np.diff(states) == -1))
print('number of DO events:', n_events)
```

Add a time-varying forcing (e.g. a slow CO2 ramp across a glacial cycle):

```python
import numpy as np
import matplotlib.pyplot as plt
from climatecritters.model_critters.bistable_melcher import BistableMelcherModel

t_arr = np.linspace(0, 599.88, 5000)
gamma_ramp = np.linspace(0.8, 3.2, 5000) # γ increases over time

model = BistableMelcherModel(
sigma=0.2,
gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)),
alpha=0.0,
)
output = model.integrate(
t_span=(t_arr[0], t_arr[-1]), y0=[1.0, 0.0],
method='heun_maruyama', dt=0.012,
kwargs={'random_seed': 0, 'si': 0.12},
)
ts = output.to_pyleo('db')
ts.plot()
plt.savefig('docs/reference/figures/BistableMelcher_ramp_example.png',
dpi=150, bbox_inches='tight')
```

References
----------
Melcher et al. (2025), Clim. Past, 21, 115-132.
https://cp.copernicus.org/articles/21/115/2025/
Melcher et al. (2025), Clim. Past, 21, 115-132. https://cp.copernicus.org/articles/21/115/2025/
"""

uses_post_history = True
Expand Down
3 changes: 3 additions & 0 deletions docs/_quarto.yml
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,8 @@ website:
href: notebooks/model_demos/g24.ipynb
- text: "Stommel"
href: notebooks/model_demos/stommel.ipynb
- text: "Melcher 2025, Bistable"
href: notebooks/model_demos/bistable_melcher.ipynb
- section: Attractors
contents:
- text: "Roessler"
Expand Down Expand Up @@ -238,6 +240,7 @@ quartodoc:
- Stocker2003BipolarSeesaw
- Stocker2003ExtendedSeaIceSeesaw
- Stommel
- BistableMelcherModel

- title: Utils

Expand Down
Binary file added docs/figures/BistableMelcher_example.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/figures/BistableMelcher_ramp_example.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified docs/figures/Lorenz96_example.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
10 changes: 10 additions & 0 deletions docs/objects.txt
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,16 @@ climatecritters.Stocker2003ExtendedSeaIceSeesaw py:class 1 api/Stocker2003Extend
climatecritters.model_critters.stocker2003_bipolar_seesaw.Stocker2003ExtendedSeaIceSeesaw py:class 1 api/Stocker2003ExtendedSeaIceSeesaw.html#climatecritters.Stocker2003ExtendedSeaIceSeesaw climatecritters.Stocker2003ExtendedSeaIceSeesaw
climatecritters.Stommel py:class 1 api/Stommel.html#climatecritters.Stommel -
climatecritters.model_critters.stommel.Stommel py:class 1 api/Stommel.html#climatecritters.Stommel climatecritters.Stommel
climatecritters.BistableMelcherModel.compute_stability_thresholds py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.compute_stability_thresholds -
climatecritters.model_critters.bistable_melcher.BistableMelcherModel.compute_stability_thresholds py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.compute_stability_thresholds climatecritters.BistableMelcherModel.compute_stability_thresholds
climatecritters.BistableMelcherModel.dydt py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.dydt -
climatecritters.model_critters.bistable_melcher.BistableMelcherModel.dydt py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.dydt climatecritters.BistableMelcherModel.dydt
climatecritters.BistableMelcherModel.populate_diagnostics_from_history py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.populate_diagnostics_from_history -
climatecritters.model_critters.bistable_melcher.BistableMelcherModel.populate_diagnostics_from_history py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.populate_diagnostics_from_history climatecritters.BistableMelcherModel.populate_diagnostics_from_history
climatecritters.BistableMelcherModel.sde_noise py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.sde_noise -
climatecritters.model_critters.bistable_melcher.BistableMelcherModel.sde_noise py:function 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel.sde_noise climatecritters.BistableMelcherModel.sde_noise
climatecritters.BistableMelcherModel py:class 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel -
climatecritters.model_critters.bistable_melcher.BistableMelcherModel py:class 1 api/BistableMelcherModel.html#climatecritters.BistableMelcherModel climatecritters.BistableMelcherModel
climatecritters.utils.noise.from_series py:function 1 api/utils.noise.from_series.html#climatecritters.utils.noise.from_series -
climatecritters.utils.noise.from_param py:function 1 api/utils.noise.from_param.html#climatecritters.utils.noise.from_param -
climatecritters.utils.resample.downsample py:function 1 api/utils.resample.downsample.html#climatecritters.utils.resample.downsample -
Expand Down
Loading