API Reference#

Solver#

pyCRO.solver.RO_solver(par, IC, N, NE, NM='EH', dt=0.1, saveat=1.0, savemethod='sampling', EF=None, noise_custom=None, verbose=True)#

Numerical solution of the Recharge Oscillator (RO) model with stochastic forcing.

This function integrates the Recharge Oscillator (RO) system numerically using either the Euler–Maruyama (EM) or Euler–Heun (EH) scheme. It includes deterministic dynamics, stochastic noise (white or red), and optional external forcing. Ensemble simulations are supported, with results returned at user-specified output intervals.

Parameters:

  • par (dict) –

    Dictionary of model parameters with the following keys (as 1-element arrays):

    • 'R'float

      Damping parameter.

    • 'F1'float

      Feedback parameter relating thermocline depth to SST.

    • 'epsilon'float

      Thermocline damping parameter.

    • 'F2'float

      Feedback parameter relating SST to thermocline.

    • 'sigma_T'float

      Noise amplitude for SST.

    • 'sigma_h'float

      Noise amplitude for thermocline depth.

    • 'B'float

      Multiplicative noise coefficient (used when n_g = 0 or 1).

    • 'n_T'int

      Noise type in SST equation (1 = white noise, 0 = red noise).

    • 'm_T'array-like

      Memory kernel for red noise in SST equation (ignored if n_T=1).

    • 'n_h'int

      Noise type in thermocline equation (1 = white noise, 0 = red noise).

    • 'm_h'array-like

      Memory kernel for red noise in thermocline equation (ignored if n_h=1).

    • 'n_g'int

      Noise structure in SST equation: 0 = multiplicative, 1 = multiplicative + Heaviside, 2 = additive.

  • IC (tuple of float) –

    Initial condition (T0, h0):

    • T0 : initial SST anomaly.

    • h0 : initial thermocline depth anomaly.

  • N (float) – Total simulation length (time units).

  • NE (int) – Number of ensemble members.

  • NM ({'EM', 'EH'}, optional) – Numerical integration method. Default is 'EH' (Euler–Heun).

  • dt (float, optional) – Numerical integration time step. Default is 0.1.

  • saveat (float, optional) – Output saving interval. Must be divisible by dt. Default is 1.0.

  • savemethod ({'sampling', 'mean'}, optional) –

    Method for saving results:

    • 'sampling' : store values every saveat steps (default).

    • 'mean' : average values within each saveat interval.

  • EF (dict, optional) –

    External forcing with the following keys (as 1D arrays of length 5):

    • 'E_T' : SST forcing coefficients.

    • 'E_h' : Thermocline forcing coefficients.

    If None (default), no external forcing is applied.

  • noise_custom ({None, int, ndarray}, optional) –

    Specification of stochastic noise:

    • None (default): generate new Gaussian noise for each ensemble.

    • int : random seed for reproducible noise (same across ensembles).

    • ndarray of shape (NT-1, 4, NE) : user-provided noise realizations.

  • verbose (bool, optional) – If True (default), print detailed setup and progress messages.

Returns:

  • T_out (ndarray of shape (N_out, NE)) – SST anomalies from the numerical integration, after applying the saving scheme.

  • h_out (ndarray of shape (N_out, NE)) – Thermocline anomalies from the numerical integration, after applying the saving scheme.

  • noise_out (ndarray of shape (NT-1, 4, NE)) – Realizations of Gaussian noise used in the simulation.

Notes

  • White noise forcing corresponds to n_T = 1 or n_h = 1.

  • Red noise forcing (n_T = 0 or n_h = 0) requires a nonzero memory kernel m_T or m_h.

  • For the Euler–Maruyama method (NM = 'EM') with multiplicative noise, an Ito-to-Stratonovich correction is applied automatically.

Examples

>>> par = {
...     'R':[0.5], 'F1':[1.2], 'epsilon':[0.3], 'F2':[0.8],
...     'sigma_T':[0.2], 'sigma_h':[0.1], 'B':[0.05],
...     'n_T':[1], 'm_T':[0], 'n_h':[1], 'm_h':[0], 'n_g':[0]
... }
>>> IC = (0.1, -0.2)
>>> T, h, noise = RO_solver(par, IC, N=50, NE=5, dt=0.1, saveat=1.0)
>>> T.shape, h.shape
((51, 5), (51, 5))
pyCRO.analytic.RO_analytic_solver(par, IC, N, NE, dt=0.1, saveat=1.0, savemethod='sampling', noise_custom=None)#

Analytical solution of the Recharge Oscillator (RO) model with deterministic and stochastic forcing. ONLY for linear RO model with white noise (annual mean parameters only)

Parameters:

  • par (dict) –

    Dictionary of model parameters, each as at least a one-element array:

    • 'R'float

      Damping parameter.

    • 'F1'float

      Feedback parameter relating thermocline depth to SST.

    • 'epsilon'float

      Thermocline damping parameter.

    • 'F2'float

      Feedback parameter relating SST to thermocline.

    • 'sigma_T'float

      Noise amplitude for SST.

    • 'sigma_h'float

      Noise amplitude for thermocline depth.

  • IC (tuple of float) –

    Initial condition (To, ho):

    • To : initial SST anomaly.

    • ho : initial thermocline depth anomaly.

  • N (float) – Total simulation length (time units).

  • NE (int) – Number of ensemble members.

  • dt (float, optional) – Numerical integration time step. Default is 0.1.

  • saveat (float, optional) – Output saving interval. Must be divisible by dt. Default is 1.0.

  • savemethod ({'sampling', 'mean'}, optional) –

    Method for saving results:

    • 'sampling' : take samples every saveat.

    • 'mean' : average values over each saveat interval.

  • noise_custom ({None, int, ndarray}, optional) –

    Specification of stochastic noise:

    • None (default): new Gaussian noise is generated for each ensemble.

    • int : seed for reproducible noise (same across ensembles).

    • ndarray of shape (NT-1, 4) : user-provided noise, repeated across ensembles.

    • ndarray of shape (NT-1, 4, NE) : user-provided noise for each ensemble.

Returns:

  • T_anal_out (ndarray of shape (N_out, NE)) – SST anomalies including deterministic and stochastic forcing, after applying the saving scheme. The final row is NaN-padded for dimensional consistency.

  • h_anal_out (ndarray of shape (N_out, NE)) – Thermocline anomalies including deterministic and stochastic forcing, after applying the saving scheme. The final row is NaN-padded.

  • noise_array (ndarray of shape (NT-1, 4, NE)) – Realizations of Gaussian noise used in the simulation.

Notes

The analytical solution consists of:

  • Deterministic part: exponential-sinusoidal functions of time.

  • Stochastic part: weighted sums (discrete convolutions) of noise with exponential and sinusoidal kernels.

Ensemble members are evaluated simultaneously using vectorized NumPy broadcasting, avoiding explicit Python loops.

Examples

>>> par = {'R':[0.5], 'F1':[1.2], 'epsilon':[0.3], 'F2':[0.8],
...        'sigma_T':[0.2], 'sigma_h':[0.1]}
>>> IC = (0.1, -0.2)
>>> T, h, noise = RO_solver_analytic(par, IC, N=50, NE=10, dt=0.1, saveat=1.0)
>>> T.shape, h.shape
((51, 10), (51, 10))
pyCRO.analytic.RO_analytic_std(par)#

Compute analytical standard deviation of T and h for the Recharge Oscillator (RO) model. ONLY for linear RO model with white noise (annual mean parameters only)

Parameters:

par (dict) – Dictionary containing parameter arrays: ‘R’, ‘F1’, ‘epsilon’, ‘F2’, ‘sigma_T’, ‘sigma_h’.

Returns:

  • T_std (float) – Standard deviation of T.

  • h_std (float) – Standard deviation of h.

Fitting#

pyCRO.fitting.RO_fitting(T, h, par_option_T, par_option_h, par_option_noise, method_fitting=None, dt=None, verbose=True)#

Fit the Recharge Oscillator (RO) model parameters to T and h time series.

This function performs parameter fitting for the RO system, using the specified prescribed parameter options and noise characteristics. Fitting can be done via linear regression (LR) or maximum likelihood estimation (MLE), with optional handling for multiplicative or additive noise. It can automatically select a default fitting method based on parameter options.

Parameters:

  • T (ndarray) – Time series of SST anomalies (1D array).

  • h (ndarray) – Time series of thermocline anomalies (1D array).

  • par_option_T (dict) – Prescribed options for SST equation parameters.

  • par_option_h (dict) – Prescribed options for thermocline equation parameters.

  • par_option_noise (dict) – Noise settings, e.g., {‘T’: ‘red’, ‘h’: ‘red’, ‘T_type’: ‘multiplicative’}.

  • method_fitting (str, optional) – Fitting method to use: ‘LR-F’, ‘LR-C’, ‘LR-F-MAC’, or ‘MLE’. If None, the method is determined automatically.

  • dt (float, optional) – Time step of the input time series. Default is 1.0 (months).

  • verbose (bool, optional) – If True, prints fitting information and progress.

Returns:

Dictionary of fitted CRO parameters: Keys include: ‘R’, ‘F1’, ‘F2’, ‘epsilon’, ‘b_T’, ‘c_T’, ‘d_T’, ‘b_h’, ‘sigma_T’, ‘sigma_h’, ‘B’, ‘m_T’, ‘m_h’, ‘n_T’, ‘n_h’, ‘n_g’.

Return type:

dict

Notes

  • Automatically removes mean from input time series.

  • Applies Ito-to-Stratonovich correction for multiplicative noise.

  • Negligible seasonal amplitudes are rounded to zero.

  • Ensures consistent noise settings for T and h equations.

  • Prints a summary of the fitting process if verbose is True.

Examples

>>> par_option_T = {'mean': 1, 'seasonal': 3, 'semi_annual': 0, ...}
>>> par_option_h = {'mean': 1, 'seasonal': 3, 'semi_annual': 0, ...}
>>> par_option_noise = {'T': 'red', 'h': 'red', 'T_type': 'multiplicative'}
>>> fitted_params = RO_fitting(T, h, par_option_T, par_option_h, par_option_noise)
>>> print(fitted_params['sigma_T'])
pyCRO.fit_LR.fit_LR(T, h, par_option_T, par_option_h, par_option_noise, dt, tend_option, fitting_option_B, fitting_option_red)#

Estimates Recharge Oscillator (RO) parameters with linear regression (LR) given ENSO SST (T) and thermocline depth (h) anomaly time series.

Parameters:

  • T (array_like) – Time series of SST anomalies (1D).

  • h (array_like) – Time series of thermocline depth anomalies (1D).

  • par_option_T (array_like) –

    Seasonal fitting options for the T-equation terms. Values must be one of:

    • 0 : constant zero

    • 1 : constant non-zero

    • 3 : annual variation

    • 5 : annual + semi-annual variation

  • par_option_h (array_like) – Seasonal fitting options for the h-equation terms (same coding as above).

  • par_option_noise (array_like of length 3) –

    Noise fitting options [n_T, n_h, n_g]:

    • n_T : 0 = red, 1 = white

    • n_h : 0 = red, 1 = white

    • n_g0 = multiplicative (full),

      1 = multiplicative (Heaviside), 2 = additive

  • dt (float) – Time step (in months).

  • tend_option ({"F", "C"}) –

    Option for numerical derivative:

    • ”F” : forward difference

    • ”C” : centered difference

  • fitting_option_B ({"MAC", "LR"}) –

    Method for estimating multiplicative noise amplitude B:

    • ”MAC” : Moment Approximation Calibration

    • ”LR” : Linear regression

  • fitting_option_red ({"LR", "AR1", "ARn"}) –

    Method for fitting red noise memory:

    • ”LR” : regression of residual derivative

    • ”AR1” : lag-1 autocorrelation

    • ”ARn” : regression using autocorrelation decay

Returns:

final_par –

Stacked parameter array including:

  • Deterministic T-equation parameters (seasonalized)

  • Deterministic h-equation parameters (seasonalized)

  • Noise parameters:

    [sigma_T, sigma_h, B, m_T, m_h, n_T, n_h, n_g]

Return type:

ndarray, shape (N, 5)

Notes

  • Seasonal terms are expanded into constant, annual, and semi-annual components depending on the fitting options.

  • Noise fitting can handle white/red and additive/multiplicative cases.

  • If multiplicative noise is selected and fitting_option_B=”MAC”, the MAC method overwrites sigma_T and B.

pyCRO.fit_MLE.fit_MLE(T, h, par_option_T, par_option_h, par_option_noise, dt)#

Estimate recharge oscillator (RO) parameters using Maximum Likelihood Estimation (MLE)

Parameters:

  • T (array_like) – Time series of SST anomalies (1D).

  • h (array_like) – Time series of thermocline depth anomalies (1D).

  • par_option_T (array_like) – Switches for which SST-related parameters to include (0=off, 1=on, 3=annual, 5=annual+semiannual).

  • par_option_h (array_like) – Switches for which h-related parameters to include (same coding).

  • par_option_noise (array_like of length 3) –

    Noise structure options [n_T, n_h, n_g]: - n_T : noise option for T-equation - n_h : noise option for h-equation - n_g : type of multiplicative noise term

    • 0 = linear multiplicative (B*T)

    • 1 = Heaviside multiplicative (B*H(T)*T)

    • 2 = additive (B=0)

  • dt (float) – Time step (months).

Returns:

par – Final estimated parameter set including deterministic RO parameters, seasonal amplitudes/phases, and red noise variances.

Return type:

ndarray

Notes

  • Implements iterative forward filtering and backward smoothing to estimate hidden noise states.

  • This function may take a long time to converge (default 10 iterations, could be increased for accuracy).

Visualization#

pyCRO.visual.plot_RO_par(par, ax=None, keys=None, ncol=4, label=None)#

Plot deterministic RO parameter seasonal cycles.

Parameters:

  • par (dict) – Dictionary of fitted CRO parameters. Each key (e.g., ‘R’, ‘F1’, etc.) maps to an array of length 1, 3, or 5.

  • ax (matplotlib axis or array-like, optional) – Axis or array of axes to plot on. If None, creates a new figure/axes.

  • keys (list of str, optional) – Which parameters to plot. If None, uses a default set.

  • ncol (int, default=4) – Number of columns for subplot layout.