lorenz84.m

% Code corresponding to the noisy Lorenz-84 model in Section 3.3 - "Atmospheric 
% Blocking Mechanisms in an Equatorial Circulation Model with Seasonality and 
% Flow-Wave Interactions" of the paper "Bridging Prediction and Attribution: 
% Identifying Forward and Backward Causal Influence Ranges Using Assimilative 
% Causal Inference".
%
% Authors: Marios Andreou, Nan Chen.
%
% Code Information: Application of conditional Assimilative Causal Inference
% (ACI) and of the associated forward and backward Causal Influence Range (CIR)
% metrics to the stochastic Lorenz-84 model for studying atmospheric blocking
% and unblocking mechanisms via a conceptual reduced-order equatorial 
% circulation model with seasonality and flow-wave interactions:
%
%     dx/dt = -y^2 - z^2 - a·(x - f(t)) + σ_x·dot(W_x),
%     dy/dt = -b·x·z + x·y - y + g + σ_y·dot(W_y),
%     dz/dt = b·x·y + x·z - z + σ_z·dot(W_z),
%
% This model describes a coarse-grained equatorial atmospheric circulation, 
% obtained from a Galerkin truncation of the two-layer quasi-geostrophic 
% potential vorticity PDEs in a channel on the equatorial β-plane (with the 
% additive noise statistically capturing the effects of unresolved scales or 
% processes). It can capture key features of the Hadley circulation, as well as 
% oceanic variability and seasonality through atmosphere-ocean coupling 
% interactions. It is also a conditional Gaussian nonlinear system (CGNS) for 
% the posterior x | (y,z); see Appendix G of the paper. The unobserved variable 
% x represents the amplitude of the westerly zonal flow (as well as the 
% zonally-averaged meridional temperature gradient after a thermal wind balance 
% argument), while the observed variables y and z represent the cosine and sine 
% phases of a large-scale baroclinic wave, corresponding to a chain of vortices 
% superimposed on the zonal flow via the -(y^2 + z^2) term in the x-equation. 
% The quadratic nonlinearities conserve energy (wave amplification at the 
% expense of the westerly current, with the wave transporting heat poleward via 
% eddies, and vice versa), with wave displacement governed by the advection 
% speed b > 1, while the Prandtl number a < 1 allows the zonal flow to dampen 
% slower than the waves. External forcings include the seasonal thermal contrast 
% f(t)= f_0 + f_1·cos(ωt), proportional to the symmetric cross-latitude heating 
% contrast between solar heating at low and high latitudes and therefore 
% stronger during the winter rather than the summer, which drives the zonal flow 
% at a one-year period (characteristic time scale of L-84 is about 5 days based 
% on the adopted parameter values). The secondary constant forcing g affects the 
% wave asymmetrically through y and mimics asymmetric contrasting thermal
% properties of the topography (zonally-alternating oceans and continents), like 
% mountain torques. This code uses the same model parameter values as those 
% cited in the paper.
%
% In this script, ACI and its forward and backward CIR metrics are employed to
% study the conditional causal relationships x(t) → y | z and x(t) → z | y over 
% time t∈[0,T], where the non-target/resolved observed channel of z or y, 
% respectively, can be chosen freely by the user. Therefore, x is treated as the
% latent candidate cause, while y or z are the target observed effects, and the
% respective remaining variable is assumed to be the resolved conditioning 
% variable (assumed to be observable as well). Operationally, for conditional
% ACI purposes, either z or y (depending on the conditional causal relationship 
% selected by the user) can be inserted into the model coefficients using its 
% observed/simulated time series, with the Bayesian update for the state 
% estimation of x using only the observations from the target observed channel y 
% or z, respectively. This is the explicit model shortcut for conditional ACI: 
% The contribution of the conditioning variable, z or y, is resolved as a 
% prescribed forcing during Bayesian inference, leaving the uncertainty 
% reduction in the state estimation of x(t) attributable to y or z solely. As 
% noted in the original ACI paper, https://doi.org/10.1038/s41467-026-68568-0, 
% this is equivalent to "letting the marginal uncertainty of z/y in the 
% likelihood to grow to infinity", which is equivalent to setting the elements 
% corresponding to z/y in the inverse of the observational noise Grammian matrix 
% to zero during state estimation. This effectively removes the influence of z/y 
% from the Bayesian update of x(t) while still retaining: (a) Its influence in 
% the dynamics and (b) the observational contributions from y/z. This latter 
% "masking" procedure is implemented in this script for calculating the 
% conditional ACI and CIR metrics for explicitness.
%
% In the paper, the conditional causal links x(t) → y | z and x(t) → z | y are
% studied over time to assess the causal structure of atmospheric blocking and
% unblocking. Based on the chosen parameter values, the model exhibits 
% bistability via two stable equilibria (node(s) or limit cycle(s)), which mimic
% alternations between atmospheric blocking and unblocking states: For a 
% dominant westerly wind jet x and small wave amplitudes y and z, we have 
% predominant zonal flow and unblocked conditions, while for weakened westerlies
% and amplified waves, which is consistent with high-pressure systems that 
% disrupt normal flow, we instead have a blocked atmosphere. Based on this, the
% CIRs are utilised to assess the temporal causal extent that the mean zonal 
% flow x has on the large-scale vortices y and z, mirroring real-world 
% teleconnections (i.e. climate anomalies related to each other at large 
% distances via atmospheric pathways). Forward CIR analysis of these 
% relationships can be used to identify when an increasing zonally-averaged 
% meridional temperature gradient (i.e. a strengthening of the westerly jet) can 
% function as a causal precursor to rapidly oscillating waves in the future, 
% while the backward CIR analysis can assess whether an observed wave amplitude 
% amplification can be traced back and causally attributed to a past weak jet 
% stream or atmospheric blocking episode (poleward heat transport). This script 
% produces the analyses associated with Figure 8 of the paper.
%
% As an added note for the Lorenz 84 mode, since the reverse conditional causal 
% relationships are also of high interest, y(t) → x | z and z(t) → x | y, the 
% same conditional ACI and CIR framework can be applied to these relationships
% but with the use of an appropriate ensemble Kalman-Bucy filter and smoother,
% since Lorenz 84 is not a CGNS for the posterior (y,z) | x. The theoretical
% and operational details are outlined in the following paper:
%   ➤ https://arxiv.org/abs/2604.25157
% with code implementations to be used as an example available at:
%   ➤ https://github.com/jiangzh67/EnKBS
% 
% This script also implements additional studies that do not appear in the paper 
% for brevity, but are included in the public repository for transparency and 
% reproducibility:
%
%     (i)   Optional calculation of the exact objective forward CIR using the 
%           ε-average definition of the subjective forward CIR, together with 
%           an exact-vs-approximate objective forward CIR comparison plot.
%     (ii)  Optional calculation of the exact objective backward CIR using the 
%           ε-average definition of the subjective backward CIR, together with 
%           an exact-vs-approximate objective backward CIR comparison plot.
%     (iii) Optional calculation of the normalised forward CIR metric δ^f(T';t) 
%           and its heatmap over natural time t and lagged observational time 
%           after t, T'-t.
%     (iv)  Optional calculation of the normalised complete backward CIR metric,
%           lim_{T'→T^-} δ^b(t;T'), using the discrete-time adaptive-lag online 
%           smoother approximation from Appendix G of the paper, and its heatmap 
%           over observational time T and backward lag T-t. Here T denotes the 
%           current online-smoother observation time, not the terminal 
%           simulation time.
%
% The exact objective forward and backward CIRs are defined in Eqs. (10) and 
% (18) of the paper, respectively, as ε-average integrals of the corresponding 
% subjective CIRs, while their computationally efficient approximations are
% correspondingly defined in Eqs. (12) and (20), as time integrals of the 
% corresponding CIR relative-entropy-based metrics (specifically, their L1-to-L∞ 
% ratio). These approximations become exact under the monotonicity conditions 
% described in the paper for the respective CIR metrics, which hold to a 
% satisfactory degree for general complex dynamical systems. In general, they 
% act as lower- and upper-bounding approximations to the exact CIRs, 
% respectively.
%
% Written and tested in MATLAB R2024b.
%
% MATLAB Toolbox and M-file Requirements:
%
% Code used to obtain the required m-file scripts and MATLAB toolboxes:
% [fList, pList] = matlab.codetools.requiredFilesAndProducts('lorenz84.m');
%
% M-file Scripts:
%   ➤ lorenz84.m
%   ➤ progress_bar.m
%   ➤ simps.m (https://www.mathworks.com/matlabcentral/fileexchange/25754-simpson-s-rule-for-numerical-integration)
%   ➤ shared_helper.m
%
% Data:
%   ➤ N/A - The script utilises synthetic observations simulated by the model
%
% Toolboxes:
%   ➤ Statistics and Machine Learning Toolbox
%
% GitHub Repository: https://github.com/marandmath/FBCIR_code
% MIT License Information: https://github.com/marandmath/FBCIR_code/blob/main/LICENSE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  USER CONTROLS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fixing the random seed for reproducibility across each simulation.
rng(17)

% Total number of time steps within the given time interval.
N = 30000;
% Numerical integration time step.
dt = 0.005;
% Total simulation time.
T = N*dt;
% Time grid (in natural time) for the time-dependent seasonal forcing.
time = (0:N)*dt;

% Plotting time window for the figures.
time_start_plot = 0;
time_end_plot = T;

% This threshold is used to treat small relative-entropy-based CIR-metric maxima 
% at each time as non-causal so as to avoid operationally inflated objective 
% forward and backward CIR lengths (both exact and approximate), as well as 
% optional normalised CIR-metric heatmaps, caused by normalising floating-point 
% noise when dividing by the L∞ norm of the respective relative-entropy-based 
% CIR metric.
RE_metric_threshold = 1e-8;

% Optional repository diagnostics; see the script's header for details. Set any
% flag to false to skip only that extra calculation/plot. The main simulation,
% ACI metric calculation, approximate objective forward and backward CIR 
% calculations, and plotting of the ACI and CIR figures from the paper remain 
% intact.
calculate_and_plot_exact_objective_forward_CIR = true;
calculate_and_plot_exact_objective_backward_CIR = true;
calculate_and_plot_normalised_forward_CIR_metric = true;
calculate_and_plot_normalised_backward_CIR_metric = true;

% Optional normalised CIR metric (forward and complete backward CIR metrics) 
% heatmap display-lag windows for optimal plotting. Defining them before the CIR 
% metric loops so that we retain only the necessary displayed metric rows for 
% this optionally plotted heatmap instead of a full O(N²) cache, without 
% changing the actual mathematical logic or displayed metric values.
% The smaller current-time-centred lagged observational time window for optimal 
% plotting (avoiding O(N²) memory-heavy plots) of the normalised forward CIR 
% metric heatmap, 
%     δ^f(T';t)/max{δ^f(T';t): T'∈[t,T]},
% over the natural time t in the plotting window [time_start_plot,time_end_plot] 
% and lagged observational time after t, T'-t ∈ [0,lag_obs_time_end_plot], 
% where:
%     lag_obs_time_end_plot = T'-t ∈ [0,time_end_plot-time_start_plot].
% Note that, by definition of the forward CIR metric, T'∈[t,T] at each t.
lag_obs_time_end_plot = 6;
% The smaller natural-time-centred, current-observational-time-based backward
% lag window for optimal plotting (avoiding O(N²) memory-heavy plots) of the 
% normalised complete backward CIR metric heatmap, 
%     lim_{T'→T^-} δ^b(t;T')/max_t{lim_{T'→T^-} δ^b(t;T')},
% over the observational time T in the plotting window [time_start_plot,
% time_end_plot] and backward lag before T, T-t ∈ [0,back_lag_time_end_plot],
% where:
%     back_lag_time_end_plot = T-t ∈ [0,time_end_plot-time_start_plot].
% Note that the t-maximum in the normalised complete backward CIR metric is 
% taken over the retained natural-time window (shortened by lag_bound; see 
% "BACKWARD CIR" section of the script below), while still considering the
% plotting window of [time_start_plot,time_end_plot], i.e.
%     t ∈ [T-lag_bound,T] 
%         = [min{time_start_plot, time_end_plot-lag_bound},time_end_plot]
%         = [{earliest retained t closer to t=0},time_end_plot] ⊆ [time_start_plot,time_end_plot].
% Note that, by definition of the (complete) backward CIR metric through a fixed
% but arbitrary observational time T, T-t is the backward lag before the current 
% observational time T.
back_lag_time_end_plot = 6;

% Conditional Lorenz-84 causal relationship to study. Choose exactly one of the
% two conditional links for each run; the same conditional ACI and CIR 
% implementation is used after changing only the target and conditioning 
% observed channels during the conditional ACI and CIR analysis. 
% Available options: "x_to_y_given_z" and "x_to_z_given_y".
conditional_relationship_to_study = "x_to_y_given_z";
% Choosing the target effect observed channel index (and by extension the 
% conditioning/non-target index as well) used by conditional ACI, as well as the
% associated plotting titles, based on the conditional relationship to be 
% studied, with the observed vector ordered as (y,z).
if conditional_relationship_to_study == "x_to_y_given_z"
    % Conditional causal relationship x(t) → y | z.
    selected_target_obs_idx = 1;
    selected_forward_relationship_title = 'x(t) \rightarrow y | z';
    selected_backward_relationship_title = 'x \rightarrow y(T) | z';
elseif conditional_relationship_to_study == "x_to_z_given_y"
    % Conditional causal relationship x(t) → z | y.
    selected_target_obs_idx = 2;
    selected_forward_relationship_title = 'x(t) \rightarrow z | y';
    selected_backward_relationship_title = 'x \rightarrow z(T) | y';
else
    error('Unknown conditional_relationship_to_study value. Use "x_to_y_given_z" or "x_to_z_given_y".')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  MODEL SETUP  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimension of the full observational vector, with the non-target z or y to be
% resolved under conditional ACI, and y or z being the target effect to be 
% studied, respectively. The observed/resolved vector is ordered as (y,z) in 
% this script. The latent candidate cause x is scalar and is therefore stored 
% using single-index arrays throughout the filtering, smoothing, and conditional 
% ACI and CIR calculations.
obs_dim = 2;

% Throughout the script, x is used for the observed variables (y,z) and y for 
% the latent variable x, so as to be consistent with the CGNS notation from 
% Appendix G of the paper. As for conditional ACI, using the notation from 
% Appendix A of the paper, the target effect variable is either x_A = y or 
% x_A = z, while the conditioning non-target variable is x_B = z or x_B = y, 
% respectively.

% Physical state variables of the Lorenz-84 model. The variable x denotes the
% unobserved zonal-flow amplitude candidate cause, while y and z denote the 
% observed wave phases (one considered to be the target effect and the other the
% resolved conditioning variable based on conditional_relationship_to_study).
x = zeros(1, N+1);
y = zeros(1, N+1);
z = zeros(1, N+1);

% The default parameter values used in the original Lorenz-84 works are adopted
% for the deterministic part of the dynamics (but with a slightly stronger 
% seasonal forcing amplitude and variance chosen). Small additive noise is 
% included to capture unresolved processes while still ensuring the 
% deterministic dynamics dominate the system state's evolution.

% Prandtl number (ratio of the zonal-flow damping time scale to the wave damping 
% time scale).
a = 1/4;
% Seasonal forcing parameters (f_0 is the mean forcing and f_1 is the amplitude 
% of the seasonal oscillation with period 2π/ω).
f_0 = 8;
f_1 = 3;
% Since the characteristic time scale of the Lorenz-84 model is about 5 days 
% with these parameter values, the seasonal forcing is chosen to have a one-year 
% period; 73*5 days = 365 days.
omega = 2*pi/73;
% Wave advection strength parameter (b > 1 ensures that the displacement of the
% wave by the zonal current can overcome the wave's amplification due to its
% interaction with the zonal flow).
b = 4;
% Constant forcing parameter that mimics asymmetric contrasting thermal effects 
% on the large-scale wave due to the topography.
g = 1;
% Additive/state-independent noise/diffusion amplitudes in the noisy Lorenz-84 
% model.
sigma_x = 0.2;
sigma_y = 0.2;
sigma_z = 0.2;

% CGNS representation with x latent and (y,z) observed; the "^x" and "^y"
% superscripts used in this script denote the observable (i.e. (y,z)) and
% unobservable (i.e. x) processes, respectively, and are used solely for
% consistency with the notation in Appendix G of the paper. (So they should not
% be confused with the x and y state variables of the Lorenz-84 model!!)
%     d(y,z) = [L^x(t,y,z)·x + f^x(t,y,z)]dt + S^x(t,y,z)·dW_x,
%     dx     = [L^y(t,y,z)·x + f^y(t,y,z)]dt + S^y(t,y,z)·dW_y.

% Observable coefficient matrix: Feedback of x in (y,z).
L_x = zeros(obs_dim, N+1);
% Forcing in the observable process.
f_x = zeros(obs_dim, N+1);
% Noise feedback matrix in the observable process.
S_x = [sigma_y, 0; 0, sigma_z];

% Unobservable coefficient matrix: Feedback of x in x.
L_y = -a;
% Forcing in the unobservable process.
f_y = zeros(1, N+1);
% Noise feedback matrix in the unobservable process.
S_y = sigma_x;

% Same initial conditions as those used in the original Lorenz-84 works.
x(1) = 1;
y(1) = 0;
z(1) = 1;

% Initiating the time-dependent model components.
L_x(:, 1) = [y(1) - b*z(1); z(1) + b*y(1)];
f_x(:, 1) = [g - y(1); -z(1)];
f_y(1) = a*(f_0 + f_1*cos(omega*time(1))) - y(1)^2 - z(1)^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%  GENERATING THE TRUE SIGNALS  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Used for the text-based progress bar.
start_time = tic;

for j = 2:N+1

    % Text-based progress bar.
    progress_bar('Simulation of the Lorenz-84 Model', j-1, N, start_time);

    % Wiener increments.
    dW_x = sqrt(dt)*randn;
    dW_y = sqrt(dt)*randn;
    dW_z = sqrt(dt)*randn;

    % Euler-Maruyama update of the state variables based on the dynamical model.
    x(j) = x(j-1) + (L_y*x(j-1) + f_y(j-1))*dt + S_y*dW_x;
    y(j) = y(j-1) + (L_x(1, j-1)*x(j-1) + f_x(1, j-1))*dt + sigma_y*dW_y;
    z(j) = z(j-1) + (L_x(2, j-1)*x(j-1) + f_x(2, j-1))*dt + sigma_z*dW_z;

    % Updating the time-dependent model components.
    L_x(:, j) = [y(j) - b*z(j); z(j) + b*y(j)];
    f_x(:, j) = [g - y(j); -z(j)];
    f_y(j) = a*(f_0 + f_1*cos(omega*time(j))) - y(j)^2 - z(j)^2;

end

% Observed physical state vector of the Lorenz-84 model, ordered as (y,z), for
% the purposes of conditional ACI and CIR analysis for the latent candidate 
% cause x.
observed_state = [y; z];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%  PLOTTING MODEL DIAGNOSTICS  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot_idx = round(time_start_plot/dt)+1:round(time_end_plot/dt)+1;
plot_time = (plot_idx - 1)*dt;

% 3D phase space and 2D phase-space projection plots illustrating the 
% wave-flow interactions geometry and Lorenz-84 chaotic attractor structure.
figure('WindowState', 'maximized');

subplot(2, 2, 1)
plot(x(plot_idx), y(plot_idx), 'k', LineWidth=1)
xlabel('x')
ylabel('y')
title('Phase Plot (x,y)')
box on
grid on
fontsize(16, 'points')

subplot(2, 2, 2)
plot(x(plot_idx), z(plot_idx), 'k', LineWidth=1)
xlabel('x')
ylabel('z')
title('Phase Plot (x,z)')
box on
grid on
fontsize(16, 'points')

subplot(2, 2, 3)
plot(y(plot_idx), z(plot_idx), 'k', LineWidth=1)
xlabel('y')
ylabel('z')
title('Phase Plot (y,z)')
box on
grid on
fontsize(16, 'points')

subplot(2, 2, 4)
plot3(x(plot_idx), y(plot_idx), z(plot_idx), 'k', LineWidth=1)
xlabel('x')
ylabel('y')
zlabel('z')
title('Phase Plot (x,y,z)')
box on
grid on
view(35, 20)
fontsize(16, 'points')

% Time series and time-averaged marginal PDFs of the latent zonal flow x and 
% the observed large-scale wave phases y and z.
figure('WindowState', 'maximized');

subplot(3, 3, [1, 2])
plot(plot_time, x(plot_idx), 'b', LineWidth=1.5)
ylabel('x')
title('Time Series of x')
box on
grid on
fontsize(16, 'points')

subplot(3, 3, 3)
xx = linspace(min(x(plot_idx)), max(x(plot_idx)), 250);
tavg_pdf_x = ksdensity(x(plot_idx), xx);
plot(xx, tavg_pdf_x, 'b', LineWidth=1.5)
hold on
mean_x = mean(x(plot_idx));
std_x = std(x(plot_idx));
plot(xx, normpdf(xx, mean_x, std_x), 'm--', LineWidth=1.5)
xlabel('x')
ylabel('p(x)')
title('PDF of x')
legend('Empirical', 'Gaussian Fit', Location='best')
box on
grid on
fontsize(16, 'points')

subplot(3, 3, [4, 5])
plot(plot_time, y(plot_idx), 'k', LineWidth=1.5)
ylabel('y')
title('Time Series of y')
box on
grid on
fontsize(16, 'points')

subplot(3, 3, 6)
yy = linspace(min(y(plot_idx)), max(y(plot_idx)), 250);
tavg_pdf_y = ksdensity(y(plot_idx), yy);
plot(yy, tavg_pdf_y, 'k', LineWidth=1.5)
hold on
mean_y = mean(y(plot_idx));
std_y = std(y(plot_idx));
plot(yy, normpdf(yy, mean_y, std_y), 'm--', LineWidth=1.5)
xlabel('y')
ylabel('p(y)')
title('PDF of y')
box on
grid on
fontsize(16, 'points')

subplot(3, 3, [7, 8])
plot(plot_time, z(plot_idx), 'm', LineWidth=1.5)
xlabel('t')
ylabel('z')
title('Time Series of z')
box on
grid on
fontsize(16, 'points')

subplot(3, 3, 9)
zz = linspace(min(z(plot_idx)), max(z(plot_idx)), 250);
tavg_pdf_z = ksdensity(z(plot_idx), zz);
plot(zz, tavg_pdf_z, 'm', LineWidth=1.5)
hold on
mean_z = mean(z(plot_idx));
std_z = std(z(plot_idx));
plot(zz, normpdf(zz, mean_z, std_z), 'm--', LineWidth=1.5)
xlabel('z')
ylabel('p(z)')
title('PDF of z')
box on
grid on
fontsize(16, 'points')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FILTERING  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Inverse of the observational noise Grammian matrix defining the weight of the
% observational influence on the Bayesian update of the unobserved scalar x(t);
% noisier observations have less impact on the posterior uncertainty reduction
% of x(t). Using this fact, for conditional ACI, this matrix is also the
% precision mask that selects which observable wave variable is allowed to 
% reduce the posterior uncertainty of x(t) and therefore be considered the 
% effect in the candidate causal link: Based on the ordering (y,z), 
% selected_target_obs_idx = 1 corresponds to target y, resolved z (x(t) → y | z) 
% and selected_target_obs_idx = 2 corresponds to target z, resolved y 
% (x(t) → z | y). Setting the elements corresponding to the resolved 
% conditioning wave variable in this matrix to zero simulates setting the 
% marginal uncertainty of this variable in the likelihood to infinity (e.g. if
% selected_target_obs_idx = 1, so we resolve z, then S_xoS_x_inv(2,2) = 0, which
% corresponds to 1/(σ_z)^2 → 0 and similarly for the (y,z) correlation in (1,2)
% (2,1) of S_xoS_x_inv, which are already zero due to the independent additive 
% noise in the dynamics, as in this limiting regime the resolved conditioned z 
% can be thought of as a diffuse or flat Bayesian prior). Therefore, we remove 
% the direct observational uncertainty reduction for x(t) from this wave 
% variable while still retaining: (a) That from the target observed wave 
% variable and (b) the conditioning wave variable's influence in the dynamics 
% of x(t).
S_xoS_x_inv = zeros(obs_dim, obs_dim);
S_xoS_x_inv(selected_target_obs_idx, selected_target_obs_idx) = 1/S_x(selected_target_obs_idx, selected_target_obs_idx)^2;

% Grammian of the unobservable process noise feedback.
S_yoS_y = S_y^2;

% Posterior filter mean of the latent variable x.
filter_mean = zeros(1, N+1);
% Initial value of the posterior filter mean.
filter_mean(1) = x(1); mu0 = filter_mean(1);
% Posterior filter covariance matrix of the latent variable x.
filter_cov = zeros(1, N+1);
% Initial value of the posterior filter covariance. Choosing a positive definite 
% matrix to preserve the positive-definiteness of the posterior covariance 
% matrices over time.
filter_cov(1) = 0.1; R0 = filter_cov(1);

% Used for the text-based progress bar.
start_time = tic;

for j = 2:N+1

    % Text-based progress bar.
    progress_bar('Filter Algorithm', j-1, N, start_time);

    dx_obs = observed_state(:, j) - observed_state(:, j-1);
    filter_gain_aux = R0*L_x(:, j-1)';

    % Update the posterior filter mean and posterior filter covariance using the
    % optimal nonlinear filter state estimation equations for CGNSs; see Theorem 
    % 2.1 in Section 2.1.2 of the Supplementary Information in the original ACI 
    % paper, https://doi.org/10.1038/s41467-026-68568-0.
    mu = mu0 + (L_y*mu0 + f_y(j-1))*dt ...
             + filter_gain_aux*S_xoS_x_inv*(dx_obs - (L_x(:, j-1)*mu0 + f_x(:, j-1))*dt);
    R = R0 + (L_y*R0 + R0*L_y' + S_yoS_y)*dt ...
           - (filter_gain_aux*S_xoS_x_inv*filter_gain_aux')*dt;

    filter_mean(j) = mu; mu0 = filter_mean(j);
    filter_cov(j) = R; R0 = filter_cov(j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SMOOTHING  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Posterior smoother mean of the latent variable x.
smoother_mean = zeros(1, N+1);
% Posterior smoother covariance matrix of the latent variable x.
smoother_cov = zeros(1, N+1);
% Smoother runs backwards: "Initial" values of the smoother statistics (i.e. at
% the last time instant T) are the corresponding posterior filter statistics.
smoother_mean(N+1) = filter_mean(N+1); muT = smoother_mean(N+1);
smoother_cov(N+1) = filter_cov(N+1); RT = smoother_cov(N+1);

% Auxiliary matrices used for the calculation of the online smoother for this
% CGNS. The online smoother is required for the calculation of the subjective
% and objective forward and backward CIR lengths of the selected Lorenz-84
% conditional causal relationship at each time t∈[0,T]. Notation used is
% consistent with that of the original CGNS online smoother work:
%   ➤ https://doi.org/10.1007/s00332-026-10271-x
% and the accompanying martingale-free introduction to CGNSs paper:
%   ➤ https://doi.org/10.3390/e27010002
E_j_matrices = zeros(1, N+1);
F_j_matrices = zeros(1, obs_dim, N+1);
G_y_j = L_y + S_yoS_y/filter_cov(N+1);
C_jj = 1 - G_y_j*dt;
E_j_matrices(N+1) = C_jj;
F_j_matrices(:, :, N+1) = G_y_j*filter_cov(N+1)*L_x(:, N+1)'*S_xoS_x_inv*dt;

% Used for the text-based progress bar.
start_time = tic;

for j = N:-1:1

    % Text-based progress bar.
    progress_bar('Smoother Algorithm', N-j+1, N, start_time);

    % Calculation of the online smoother auxiliary matrices.
    G_y_j = L_y + S_yoS_y/filter_cov(j);
    C_jj = 1 - G_y_j*dt;
    E_j_matrices(j) = C_jj;
    F_j_matrices(:, :, j) = G_y_j*filter_cov(j)*L_x(:, j)'*S_xoS_x_inv*dt;

    % Auxiliary matrices for the posterior smoother mean and covariance updates.
    A_j = L_y;
    B_j = S_yoS_y;

    % Update the posterior smoother mean and posterior smoother covariance using 
    % the optimal nonlinear smoother state estimation backward equations for 
    % CGNSs; see Theorem 2.2 in Section 2.1.2 of the Supplementary Information 
    % in the original ACI paper, https://doi.org/10.1038/s41467-026-68568-0.
    mu = muT - (L_y*muT + f_y(j) - B_j/filter_cov(j)*(filter_mean(j) - muT))*dt;
    R = RT - ((A_j + B_j/filter_cov(j))*RT + RT*(A_j + B_j/filter_cov(j))' - B_j)*dt;

    % Alternative and equivalent discrete-time smoother update equations for 
    % CGNSs using the online smoother auxiliary matrices. See the 
    % martingale-free introduction to CGNSs paper for more details:
    %   ➤ https://doi.org/10.3390/e27010002
    % dx_obs = observed_state(:, j+1) - observed_state(:, j);
    % E_j = E_j_matrices(j);
    % F_j = F_j_matrices(:, :, j);
    % mu = filter_mean(j) ...
    %      + E_j*(muT - (1 + L_y*dt)*filter_mean(j) - f_y(j)*dt) ...
    %      + F_j*(dx_obs - (L_x(:, j)*filter_mean(j) + f_x(:, j))*dt);
    % R = filter_cov(j) ...
    %     + E_j*(RT*E_j' - (1 + L_y*dt)*filter_cov(j)) ...
    %     - F_j*L_x(:, j)*filter_cov(j)*dt;

    smoother_mean(j) = mu; muT = smoother_mean(j);
    smoother_cov(j) = R; RT = smoother_cov(j);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%  PLOTTING FILTER AND SMOOTHER RESULTS  %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Time series of the posterior filter and smoother Gaussian statistics of the
% latent candidate cause x after the conditional ACI precision masking (based
% on the chosen conditional causal relationship).
figure('WindowState', 'maximized');

plot(plot_time, x(plot_idx), 'b', LineWidth=1.5, DisplayName='Truth')
hold on
plot(plot_time, filter_mean(plot_idx), 'g', LineWidth=1.5, DisplayName='Filter')
plot(plot_time, smoother_mean(plot_idx), 'r', LineWidth=1.5, DisplayName='Smoother')
filter_upper = filter_mean(plot_idx) + 2*sqrt(filter_cov(plot_idx));
filter_lower = filter_mean(plot_idx) - 2*sqrt(filter_cov(plot_idx));
smoother_upper = smoother_mean(plot_idx) + 2*sqrt(smoother_cov(plot_idx));
smoother_lower = smoother_mean(plot_idx) - 2*sqrt(smoother_cov(plot_idx));
patch([plot_time, fliplr(plot_time)], [filter_lower, fliplr(filter_upper)], ...
      'g', FaceAlpha=0.2, LineStyle='none', DisplayName='±2Std Filter')
patch([plot_time, fliplr(plot_time)], [smoother_lower, fliplr(smoother_upper)], ...
      'r', FaceAlpha=0.2, LineStyle='none', DisplayName='±2Std Smoother')
xlabel('t')
ylabel('x')
title(sprintf('Time Series of x and of its Filter and Smoother Posterior Statistics (%s)', selected_forward_relationship_title))
legend(Location='best')
box on
grid on
fontsize(16, 'points')

% Comparison of the time-averaged low-order statistical diagnostics for the 
% posterior state estimates of x (filter and smoother after the conditional ACI 
% masking based on the chosen conditional causal relationship); normalised RMSE 
% (NRMSE) and Pearson correlation coefficient.
nrmse_filter = rmse(x(plot_idx), filter_mean(plot_idx))/std(x(plot_idx));
nrmse_smoother = rmse(x(plot_idx), smoother_mean(plot_idx))/std(x(plot_idx));
relative_nrmse_change = 100*(nrmse_filter - nrmse_smoother)/nrmse_filter;
fprintf('Relative percentage change in smoother NRMSE from filter NRMSE for x in %s = %0.2f%% (+ = decreased, better)\n', strrep(selected_forward_relationship_title, '\rightarrow', '→'), relative_nrmse_change);
corr_filter_matrix = corrcoef(x(plot_idx), filter_mean(plot_idx));
corr_smoother_matrix = corrcoef(x(plot_idx), smoother_mean(plot_idx));
corr_filter = corr_filter_matrix(1, 2);
corr_smoother = corr_smoother_matrix(1, 2);
relative_corr_change = corr_smoother - corr_filter;
fprintf('Difference/Δ in smoother correlation from filter correlation for x in %s = %0.2f (+ = increased, better)\n', strrep(selected_forward_relationship_title, '\rightarrow', '→'), relative_corr_change);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ACI ANALYSIS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Theoretically, since the partially-observed stochastic dynamical system 
% studied in this script is a real-valued CGNS, then the posterior filter and 
% smoother covariances should remain, almost surely over the space of 
% observable sample paths, finite, real-valued, symmetric, and positive 
% semi-definite (PSD) at all times (or positive definite (PD) at all times if 
% the filter covariance at the initial time t=0 is positive definite as well): 
%   ➤ https://doi.org/10.3390/e27010002
%   ➤ https://doi.org/10.1007/978-3-662-10028-8
% However, due to finite numerical floating-point precision and/or because of 
% discretisation errors (since we solve the continuous-time CGNS posterior 
% filter/smoother state estimation equations in discrete time using the 
% Euler-Maruyama method and a discretisation time step of dt=Δt), then the 
% produced posterior filter/smoother covariances may violate these mathematical
% properties at some time instants. This is especially true for stiff update 
% equations which might require very small time steps dt=Δt and for the smoother
% covariance, which is calculated backwards in time using the inverse of the
% filter covariance, and is therefore more prone to numerical instabilities. 
% Still, these violations are extremely rare and usually occur only at one or
% two time instants, while also being very negligible. Still, to avoid script
% errors of this type to terminate the execution of the code, we regularise the
% posterior filter and/or smoother covariances only when such a finite numerical
% instability makes them non-real, non-symmetric, or non-PSD, and before any of
% the relative-entropy-based ACI and CIR calculations. This is achieved through 
% a finite real-symmetric PD projection via a machine-scale eigenvalue or 
% variance lift on the posterior covariance only whenever it is necessary. Valid 
% posterior covariances over time remain unchanged. Also, for transparency and 
% as a diagnostic summary, a brief report of the instances of these corrections 
% and of any non-finite values is printed to the console for the user's 
% information; a non-finite covariance naturally raises a clear error because it 
% cannot be regularised without an intrusive intervention to the underlying
% mathematical, theoretical, and operational logic of the ACI and CIR framework.
filter_cov = shared_helper('regularise_posterior_covariance_history', filter_cov, dt, 'Filter');
smoother_cov = shared_helper('regularise_posterior_covariance_history', smoother_cov, dt, 'Smoother');

% Calculating the ACI metric for the selected conditional relationship at each
% time in the plotting window, i.e. the relative entropy from the posterior
% smoother to the filter of x(t) after the conditional ACI masking of the
% non-target observed wave variable: 
%     P(p_t^{s|z}(x|y), p_t^{f|z}(x|y)) (x(t) → y | z),
%     P(p_t^{s|y}(x|z), p_t^{f|y}(x|z)) (x(t) → z | y),
% which, due to the Gaussianity of these modified posterior distributions, both 
% divergences are given by the signal-dispersion formula; see Appendix G of the 
% paper.
signal_smoother_filter = 0.5*(smoother_mean(plot_idx) - filter_mean(plot_idx)).^2./filter_cov(plot_idx);
cov_ratio_smoother_filter = smoother_cov(plot_idx)./filter_cov(plot_idx);
dispersion_smoother_filter = 0.5*(cov_ratio_smoother_filter - log(cov_ratio_smoother_filter) - 1);
ACI_metric = signal_smoother_filter + dispersion_smoother_filter;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FORWARD CIR  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implementation of the fixed-lag online smoother for CGNSs; see Appendix G of 
% the paper and original CGNS online smoother paper for the theoretical and 
% implementation details:
%   ➤ https://doi.org/10.1007/s00332-026-10271-x

% The fixed-lag online smoother distributions, 
% p(x(t_j)|target(s<=t_n); conditioning), after the conditional ACI masking of 
% the non-target wave variable where target = y or z and conditioning = z or y, 
% respectively, based on the chosen conditional relationship, are needed to 
% evaluate the forward CIR metric (and by extension the subjective and
% objective forward CIRs) for the selected Lorenz-84 conditional relationship as
% a function of the lagged observational time t_n<=T at each fixed natural time
% t_j. In the following implementation, the relative entropy from the complete
% smoother distribution p(x(t_j)|target(s<=t_N); conditioning) to the lagged
% smoother distribution p(x(t_j)|target(s<=t_n); conditioning) is calculated at
% each natural time t_j for all lagged observational times t_j<=t_n<=t_N=T,
%     P^j_n := P(p(x(t_j)|target(s<=t_N); conditioning), p(x(t_j)|target(s<=t_n); conditioning)),  
% where both posterior distributions are Gaussian for this CGNS. See Appendix G 
% of the paper for the background, notations, and details.

% The fixed-lag parameter is set equal to the total number of observations, 
% N = ⌈T/Δt⌉, such that at each time instant the full backward algorithm is 
% carried out. This is because each conditional ACI masked online smoother 
% distribution, p(x(t_j)|target(s<=t_n); conditioning), is needed for the 
% calculation of the subjective and objective forward CIRs; see Appendix G of 
% the paper. (N+1 due to the initial condition.)
fixed_lag = N+1;

% Saving the online smoother mean, online smoother covariance matrices, and 
% update matrices in a cell array where each row is another cell array with as 
% many columns as the cardinal number of the current row. Using such nested cell
% arrays efficiently simulates staggered arrays in MATLAB. This approach 
% efficiently preserves space in memory without defining unnecessarily large 
% high-order tensors to store the online smoother estimations and update 
% matrices. In these nested cell arrays, the first/parent index corresponds to 
% n∈{j,j+1,...,N}, for the current target effect observation y^n = y(t_n) or
% z^n = z(t_n), while the second/child index corresponds to j∈{0,1,...,N}, for 
% the time instant t_j at which we carry out the online smoother state 
% estimation for x^j = x(t_j).
online_fixed_mean = cell(N+1, 1);
online_fixed_cov = cell(N+1, 1);
update_matrices_fixed = cell(N-1, 1);
for n = 1:(N-1)
    update_matrices_fixed{n} = zeros(1, n+1);
    online_fixed_mean{n} = zeros(1, n);
    online_fixed_cov{n} = zeros(1, n);
end
for n = N:N+1
    online_fixed_mean{n} = zeros(1, n);
    online_fixed_cov{n} = zeros(1, n);
end

% Need to do the first two observations manually.

% A single observation (n=1).
online_fixed_mean{1}(1) = filter_mean(1);
online_fixed_cov{1}(1) = filter_cov(1);

% Two observations (n=2).
online_fixed_mean{2}(2) = filter_mean(2);
online_fixed_cov{2}(2) = filter_cov(2);
if fixed_lag == 0

    online_fixed_mean{2}(1) = online_fixed_mean{1}(1);
    online_fixed_cov{2}(1) = online_fixed_cov{1}(1);

else

    dx_obs = observed_state(:, 2) - observed_state(:, 1);
    aux_vec = filter_mean(1) ...
              - E_j_matrices(1)*((1 + L_y*dt)*filter_mean(1) + f_y(1)*dt) ...
              + F_j_matrices(:, :, 1)*(dx_obs - (L_x(:, 1)*filter_mean(1) + f_x(:, 1))*dt);
    online_fixed_mean{2}(1) = E_j_matrices(1)*filter_mean(2) + aux_vec;
    aux_mat = filter_cov(1) ...
              - E_j_matrices(1)*(1 + L_y*dt)*filter_cov(1) ...
              - F_j_matrices(:, :, 1)*L_x(:, 1)*filter_cov(1)*dt;
    online_fixed_cov{2}(1) = E_j_matrices(1)*filter_cov(2)*E_j_matrices(1)' + aux_mat;

end

% Used for the text-based progress bar.
start_time = tic;

for n = 3:N+1

    % Text-based progress bar.
    progress_bar('Fixed-Lag Online Smoother Algorithm for Forward CIR', n-2, N-1, start_time);

    online_fixed_mean{n}(n) = filter_mean(n);
    online_fixed_cov{n}(n) = filter_cov(n);

    if fixed_lag == 0

        online_fixed_mean{n}(n-1) = online_fixed_mean{n-1}(n-1);
        online_fixed_cov{n}(n-1) = online_fixed_cov{n-1}(n-1);

    else

        dx_obs = observed_state(:, n) - observed_state(:, n-1);
        aux_vec = filter_mean(n-1) ...
                  - E_j_matrices(n-1)*((1 + L_y*dt)*filter_mean(n-1) + f_y(n-1)*dt) ...
                  + F_j_matrices(:, :, n-1)*(dx_obs - (L_x(:, n-1)*filter_mean(n-1) + f_x(:, n-1))*dt);
        online_fixed_mean{n}(n-1) = E_j_matrices(n-1)*filter_mean(n) + aux_vec;
        aux_mat = filter_cov(n-1) ...
                  - E_j_matrices(n-1)*(1 + L_y*dt)*filter_cov(n-1) ...
                  - F_j_matrices(:, :, n-1)*L_x(:, n-1)*filter_cov(n-1)*dt;
        online_fixed_cov{n}(n-1) = E_j_matrices(n-1)*filter_cov(n)*E_j_matrices(n-1)' + aux_mat;

    end

    for j = (n-1):-1:1

        if (1 <= j) && (j <= n-1-fixed_lag)

            online_fixed_mean{n}(j) = online_fixed_mean{n-1}(j);
            online_fixed_cov{n}(j) = online_fixed_cov{n-1}(j);

        elseif (n-fixed_lag <= j) && (j <= n-1)

            if j == n-1
                update_matrices_fixed{n-2}(n-1) = 1;
            elseif j == n-2
                update_matrices_fixed{n-2}(n-2) = E_j_matrices(n-2);
            else
                update_matrices_fixed{n-2}(j) = update_matrices_fixed{n-3}(j)*E_j_matrices(n-2);
            end
            online_mean_innovation = online_fixed_mean{n}(n-1) - filter_mean(n-1);
            online_cov_innovation = online_fixed_cov{n}(n-1) - filter_cov(n-1);
            online_fixed_mean{n}(j) = online_fixed_mean{n-1}(j) + update_matrices_fixed{n-2}(j)*online_mean_innovation;
            online_fixed_cov{n}(j) = online_fixed_cov{n-1}(j) + update_matrices_fixed{n-2}(j)*online_cov_innovation*update_matrices_fixed{n-2}(j)';

        end

    end

end

% Letting 10⁻⁶ <= ε <= 10¹ with a resolution of 513 points.
epsilon_resolution = 513;
lowest_order = -6;
highest_order = 1;
eps_ord_values = flip(linspace(lowest_order, highest_order, epsilon_resolution));
% Threshold values used when evaluating the subjective CIR lengths.
CIR_epsilon_values = 10.^eps_ord_values;
epsilon_axis_tick_values = unique([lowest_order:1:floor(highest_order), highest_order]);

% The following snippet uses a more adaptive mesh of ε instead of the purely 
% logarithmic one used in the above implementation, for a more realistic
% integration of the subjective CIRs over ε for obtaining the exact
% corresponding objective CIR length at each time.
% adaptive_point = -2;
% half_resolution = 250;
% eps_ord_values = unique([
%     flip(log10(linspace(10^adaptive_point, 10^highest_order, half_resolution))), ...
%     flip(linspace(lowest_order, adaptive_point, half_resolution))
% ], 'stable');

% Calculating the subjective and objective CIRs over the plotting time interval 
% of choice. We add a lookahead tolerance for the lagged observational time: 
%     T'∈[t,time_end_plot+lookahead_tolerance],
% to avoid observational saturation as t approaches time_end_plot. If the
% plotting interval reaches the final simulation time, the available lookahead
% is set automatically to zero.
lookahead_tolerance = 0.6;
first_idx = round(time_start_plot/dt)+1;
target_last_idx = round(time_end_plot/dt)+1;
if time_end_plot+lookahead_tolerance < T
    last_idx = round((time_end_plot+lookahead_tolerance)/dt)+1;
else
    last_idx = target_last_idx;
end
lookahead_steps = last_idx - target_last_idx;
calc_len = length(first_idx:last_idx);
plot_len = calc_len - lookahead_steps;

% Calculating the subjective forward CIR length for the selected conditional
% relationship at each time t in the plotting interval and for various orders
% O(10⁻ᵏ) of ε values. The associated objective forward CIR length is also
% calculated using: (a) Its definition by integrating the respective subjective
% forward CIR length over ε and (b) its computationally efficient
% lower-bounding/underestimating approximation through the time integral formula
% of the forward CIR metric (specifically its L1-to-L∞ ratio). The theory behind
% the subjective and objective forward CIR length and the latter's approximation
% is given in Sections 2.3.2-2.3.3 of the paper, while the computational details
% for CGNSs are given in Appendix G.
forw_subjective_CIR = zeros(length(eps_ord_values), calc_len);
forw_approx_objective_CIR = zeros(1, calc_len);

% Forward CIR relative-entropy metric δ^f(T';t) (see Section 2.3.1 and 2.3.2 of
% the paper):
%     δ^f(T';t) := δ(t,T') = P(p(x(t)|target(s<=T); conditioning), p(x(t)|target(s<=T'); conditioning)),
% where the non-target wave variable has been conditionally masked in the
% posterior distributions under the conditional ACI framework. This metric is 
% used to calculate the approximate objective forward CIR via a time integral 
% over the lagged observational time T' instead of integrating the associated 
% subjective forward CIR over ε as in the definition to get the exact objective
% forward CIR length (as well as to calculate the subjective forward CIR 
% itself). Using the notation from Appendix G of the paper, RE_n below stores:
%     P^j_n = δ(t_j,t_n) = δ^f(t_n;t_j).
% The columns of the 1D RE_n vector correspond to the lagged observational time 
% T' indexed by t_n (i.e. n∈{j,j+1,...,last_idx}) and calculated at each natural 
% time t indexed by t_j (i.e. j∈{first_idx,first_idx+1,...,last_idx}, where 
% first_idx and last_idx are determined by the plotting interval and the 
% lookahead tolerance), iterated over in the loop below. Note that this forward 
% CIR metric vector is calculated regardless of whether 
% calculate_and_plot_normalised_forward_CIR_metric = false, since its values are 
% required for computing the subjective forward CIRs, the definition-based 
% (exact) objective forward CIR, and its computationally efficient objective 
% approximation.
if calculate_and_plot_normalised_forward_CIR_metric
    % The optional normalised forward CIR metric heatmap cache retains only the 
    % requested displayed observational lag rows, with the rows indexed by the 
    % current-natural-time-centred lagged observational time T'-t, induced by 
    % the moving natural time t=t_j in the plotting interval (+ lookahead 
    % tolerance at the right end-point; j∈{first_idx, ..., last_idx}) and the 
    % lagged observational time T' ([obs]∈{j,...,last_idx}), and the columns 
    % indexed by the moving natural time t in the plotting interval. This 
    % approach avoids storing a full memory-heavy O(N²) CIR metric cache; 
    % specifically this retains only the portion of the forward CIR metric 
    % calculated at each natural time t=t_j in the loop below, under RE_n, that
    % is necessary for plotting the normalised heatmap, leaving its calculated
    % values unchanged and omitted lag rows as NaN. The computations of the
    % displayed CIR metric values are unchanged by this convention. Cache 
    % columns retain j=first_idx:target_last_idx only; the right end-point 
    % lookahead tolerance remains in RE_n and all forward CIR calculations, but 
    % it is not displayed.
    forw_lag_steps = min(round(lag_obs_time_end_plot/dt)+1, plot_len);
    % Not to be confused with the RE_n variable, which is the 1D vector of the 
    % forward CIR metric values at each natural time t=t_j in the plotting 
    % interval, iterated over in the loop below; see the comments above for more 
    % details. Essentially, this cache stores the RE_n values at each natural 
    % time t=t_j in the plotting interval and only for the retained lagged 
    % observational time rows chosen by the user through lag_obs_time_end_plot, 
    % with the omitted lag rows set to NaN.
    forw_RE_metric = nan(forw_lag_steps, plot_len);
end
max_forw_RE_metric = zeros(1, calc_len);

% Used for the text-based progress bar.
start_time = tic;

for j = first_idx:last_idx

    % Text-based progress bar.
    progress_bar('Calculation of the Forward CIRs', j-first_idx+1, calc_len, start_time);

    % Calculating and storing P^j_n over n∈{j,j+1,...,last_idx} for a fixed
    % j∈{first_idx,first_idx+1,...,last_idx}.
    RE_n = zeros(1, length(j:last_idx));
    for obs = j:last_idx
        cov_ratio = online_fixed_cov{end}(j)/online_fixed_cov{obs}(j);
        if ~isreal(cov_ratio) || ~isfinite(cov_ratio) || cov_ratio <= 0
            error('FBCIR:InvalidCovarianceRatio', 'The covariance ratio for calculating the forward CIR metric and lengths is invalid at natural time index %d and lagged observational time index %d.', j, obs)
        end
        RE_n(obs-j+1) = 0.5*(online_fixed_mean{end}(j) - online_fixed_mean{obs}(j))^2/online_fixed_cov{obs}(j) ...
                        + 0.5*(cov_ratio - log(cov_ratio) - 1);
    end
    % Store the forward CIR metric row and its maximum for the current natural
    % time t_j; the maximum is reused for objective forward CIR normalisation 
    % and for the optional normalised forward CIR metric heatmap.
    if calculate_and_plot_normalised_forward_CIR_metric && j <= target_last_idx
        forw_lag_cache_len = min(forw_lag_steps, length(RE_n));
        forw_RE_metric(1:forw_lag_cache_len, j-first_idx+1) = RE_n(1:forw_lag_cache_len);
    end
    max_forw_RE_metric(j-first_idx+1) = max(RE_n);

    % If the maximum is essentially zero, corresponding to a non-causal 
    % relationship at this time, then do not calculate the CIR and set
    % equal to 0 instead (to avoid operationally inflated CIR values
    % due to the normalisation in the objective CIR approximation and
    % to numerical precision errors).
    if max(RE_n) > RE_metric_threshold
        % Efficient objective forward CIR approximation. It is a lower-bound
        % approximation in general and becomes exact under the monotonicity
        % condition for the forward CIR metric described in the paper.

        % Estimation of the forward CIR metric's integral using the 
        % L1 grid-function norm from Appendix G of the paper since we deal with
        % discrete-time quantities with the fixed-lag CGNS online smoother.
        forw_approx_objective_CIR(j-first_idx+1) = sum(RE_n)*dt/max(RE_n);

        % Smooth quadrature: Estimation of the forward CIR metric's integral 
        % using a composite trapezoidal rule.
        % forw_approx_objective_CIR(j-first_idx+1) = trapz(RE_n)*dt/max(RE_n);

        % Smoother quadrature: Estimation of the forward CIR metric's integral 
        % using a composite Simpson's 1/3 rule.
        % forw_approx_objective_CIR(j-first_idx+1) = simps(RE_n)*dt/max(RE_n);
    else
        forw_approx_objective_CIR(j-first_idx+1) = 0;
    end

    % Calculation of the subjective forward CIR lengths for all ε values. The
    % reverse cumulative maximum makes each threshold crossing equivalent to
    % find(RE_n > ε, 1, 'last') without repeatedly scanning the full row.

    % Indexing convention: RE_n(1) corresponds to T'=t_j, so subj_CIR_idx counts
    % the retained grid samples up to the last ε-threshold exceedance in 
    % (t_j,T]. Multiplying by Δt (dt) follows the L1 grid-function norm 
    % convention outlined in Appendix G of the paper due to the discrete-time 
    % implementation adopted in the code, matching the dt*sum(...) L1 
    % grid-function norm in the computationally efficient objective 
    % approximation above. An interval-endpoint convention ([t_j,T]) could be 
    % used instead with max(subj_CIR_idx-1, 0)*dt, such that an exceedance only 
    % at T'=t_j has zero length; still, we keep the grid-function norm 
    % convention for consistency with the discrete-time objective forward CIR 
    % approximation.
    tail_max_RE = cummax(RE_n, 'reverse');
    subj_CIR_idx = zeros(length(eps_ord_values), 1);
    last_exceedance_idx = 0;
    for eps_idx = 1:length(eps_ord_values)
        while (last_exceedance_idx < length(tail_max_RE)) && ...
              (tail_max_RE(last_exceedance_idx+1) > CIR_epsilon_values(eps_idx))
            last_exceedance_idx = last_exceedance_idx + 1;
        end
        subj_CIR_idx(eps_idx) = last_exceedance_idx;
    end
    forw_subjective_CIR(:, j-first_idx+1) = subj_CIR_idx*dt;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  NOTE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% When choosing a sufficiently large epsilon_resolution and a sufficiently large 
% k where 10⁻ᵏ ≤ ε, then the objective CIR at each time tⱼ can be calculated 
% through the definition by averaging the corresponding subjective CIR over ε 
% via a numerical quadrature method using the following command (the flip 
% operations are needed to put the ε (see eps_ord_values variable) interval in 
% ascending order):

% Estimation of the integral using a composite trapezoidal rule (a mask using 
% RE_metric_threshold, like the one utilised below and in the approximate 
% objective CIR calculation, can be considered to avoid operationally inflated 
% CIR values in the following due to the normalisation by the maxima).
% forw_defn_objective_CIR = trapz(10.^flip(eps_ord_values), flipud(forw_subjective_CIR(:, 1:plot_len)), 1)./max_forw_RE_metric(1:plot_len);

% Estimation of the integral using a composite Simpson's 1/3 rule for better
% accuracy. The same RE_metric_threshold convention used in the approximate 
% objective CIR calculation is applied below before normalisation of this 
% definition-based/exact length by the maxima, so as to avoid operationally 
% inflated CIR values.
forw_RE_metric_mask = max_forw_RE_metric(1:plot_len) > RE_metric_threshold;
if calculate_and_plot_exact_objective_forward_CIR
    forw_defn_objective_CIR = simps(10.^flip(eps_ord_values), flipud(forw_subjective_CIR(:, 1:plot_len)), 1);
    forw_defn_objective_CIR(~forw_RE_metric_mask) = 0;
    forw_defn_objective_CIR(forw_RE_metric_mask) = forw_defn_objective_CIR(forw_RE_metric_mask)./max_forw_RE_metric(forw_RE_metric_mask);
end

% Here, by "exact", we mean "definition-based" as opposed to 
% "approximation-based" (the latter being the computationally efficient 
% objective forward CIR approximation), since numerical quadrature with a 
% truncated ε interval is still used for the integral, which is sufficiently 
% accurate for a large enough ε resolution and coverage.

% The fixed-lag online smoother history is no longer needed after all forward
% CIR metric and objective/subjective CIR-length calculations are complete.
% Clear it before allocating the adaptive-lag history to reduce peak memory use.
clear online_fixed_mean online_fixed_cov update_matrices_fixed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  BACKWARD CIR  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Implementation of the adaptive-lag online smoother for CGNSs; see Appendix G 
% of the paper and original CGNS online smoother paper for the theoretical and 
% implementation details:
%   ➤ https://doi.org/10.1007/s00332-026-10271-x

% The adaptive-lag online smoother distributions, 
% p^{upd}(x(t_j)|target(s<=t_N); conditioning), obtained from updating 
% p^{lag}(x(t_j)|target(s<=t_{N-1}); conditioning) due to the newly observed 
% effect target(t_N) = target(T) using the online smoother's information gain 
% criterion (see the original CGNS online smoother paper and Appendix G of the 
% paper for context here), as well as the latter lagged distributions, after the 
% conditional ACI masking of the non-target wave variable where target = y or z 
% and conditioning = z or y, respectively, based on the chosen conditional 
% relationship, are needed to evaluate the (complete) backward CIR metric (and 
% by extension the subjective and objective backward CIRs) as a function of 
% natural time t_j<=T at each observational time t_N=T. In the following 
% implementation, the relative entropy from the updated adaptive-lag online 
% smoother distribution p^{upd}(x(t_j)|target(s<=t_N); conditioning) to the 
% lagged online smoother distribution 
% p^{lag}(x(t_j)|target(s<=t_{N-1}); conditioning) is calculated at each 
% observational time t_N=T for all natural times 0<=t_j<=T,
%     P^j_T := P(p^{upd}(x(t_j)|target(s<=t_N); conditioning), p^{lag}(x(t_j)|target(s<=t_{N-1}); conditioning)), 
% where both posterior distributions are Gaussian for this CGNS; note how using
% the adaptive-lag online smoother distributions going from target(t_{N-1}) to 
% target(t_N) is approximating the limit of complete information incorporation 
% for the complete backward CIR metric, from the theory of the backward CIR in 
% Section 2.3.4 of the paper. See Appendix G of the paper for the background, 
% notations, and details.

% While theoretically an upper-lag bound of N = ⌈T/Δt⌉ (for each observational 
% time T; not to be confused with the terminal simulation time), lag_bound = N, 
% should be used to make sure the adaptive lag (and by extension the backward 
% CIR lengths) are fully captured and not prematurely truncated, practically a 
% much smaller lag bound is sufficient and can be used due to the chaotic nature 
% of the system to reduce the computational burden of the backward CIR algorithm 
% (since new observations affect the Bayesian state estimation of the unobserved 
% variables for a finite and relatively short interval in the past). When the 
% t=0 baseline falls outside this shortened lag window, the earliest retained 
% lag (i.e. closest to t=0) is used as the finite-memory approximation of the 
% baseline term in the centred (complete) backward CIR metric. Still, for 
% consistency with the results of the paper, we use the full lag_bound of N 
% (although lag_bound = round(3/dt) to round(4/dt) has been shown to be 
% sufficient during testing).
lag_bound = N;

% Saving the online smoother mean, online smoother covariance matrices, and 
% update matrices in a cell array where each row is another cell array with as 
% many columns as the cardinal number of the current row. Using such nested cell
% arrays efficiently simulates staggered arrays in MATLAB. This approach 
% efficiently preserves space in memory without defining unnecessarily large 
% high-order tensors to store the online smoother estimations and update 
% matrices. In these nested cell arrays, the first/parent index corresponds to 
% n∈{j,j+1,...,N}, for the current effect observation y^N = y(t_N) = y(T) or
% z^N = z(t_N) = z(T), while the second/child index corresponds to 
% j∈{0,1,...,N}, for the time instant t_j at which we carry out the online 
% smoother state estimation for x^j = x(t_j).
online_adapt_mean = cell(N+1, 1);
online_adapt_cov = cell(N+1, 1);
update_matrices_adapt = cell(N-1, 1);
for n = 1:(N-1)
    update_matrices_adapt{n} = zeros(1, n+1);
    online_adapt_mean{n} = zeros(1, n);
    online_adapt_cov{n} = zeros(1, n);
end
for n = N:N+1
    online_adapt_mean{n} = zeros(1, n);
    online_adapt_cov{n} = zeros(1, n);
end

% Need to do the first two observations manually.

% A single observation (n=1).
online_adapt_mean{1}(1) = filter_mean(1);
online_adapt_cov{1}(1) = filter_cov(1);

% Two observations (n=2).
online_adapt_mean{2}(2) = filter_mean(2);
online_adapt_cov{2}(2) = filter_cov(2);
dx_obs = observed_state(:, 2) - observed_state(:, 1);
aux_vec = filter_mean(1) ...
          - E_j_matrices(1)*((1 + L_y*dt)*filter_mean(1) + f_y(1)*dt) ...
          + F_j_matrices(:, :, 1)*(dx_obs - (L_x(:, 1)*filter_mean(1) + f_x(:, 1))*dt);
online_adapt_mean{2}(1) = E_j_matrices(1)*filter_mean(2) + aux_vec;
aux_mat = filter_cov(1) ...
          - E_j_matrices(1)*(1 + L_y*dt)*filter_cov(1) ...
          - F_j_matrices(:, :, 1)*L_x(:, 1)*filter_cov(1)*dt;
online_adapt_cov{2}(1) = E_j_matrices(1)*filter_cov(2)*E_j_matrices(1)' + aux_mat;

% Calculating the subjective and objective backward CIRs over the plotting time 
% interval of choice. Importantly, due to the nature of the backward CIR theory, 
% unlike its forward counterpart, everything can be computed efficiently in the
% same loop in an online/real-time manner: The adaptive-lag CGNS online 
% smoother, the (complete) backward CIR metric, and the subjective and objective 
% backward CIR lengths (both exact and approximate for the latter). Also, unlike 
% in the forward CIR calculations, no lookahead in observational time is added 
% here: At each observation time T=t_N, the backward metric only uses the 
% retained natural-time window before T, with its maximum length controlled by 
% lag_bound. So:
% back_first_idx = first_idx;
% back_last_idx = target_last_idx;
back_plot_len = length(first_idx:target_last_idx);

% Calculating the subjective backward CIR length for the selected conditional
% relationship at each observational time T in the plotting interval, for the
% same O(10⁻ᵏ) orders of ε values as in the forward CIR implementation. The
% associated objective backward CIR length is also calculated using: (a) Its
% definition by integrating the respective subjective backward CIR length over ε
% and (b) its computationally efficient upper-bounding/overestimating 
% approximation through the time integral formula of the complete backward CIR 
% metric (specifically its L1-to-L∞ ratio). The theory behind the subjective and 
% objective backward CIR length and the latter's approximation is given in 
% Sections 2.3.4-2.3.6 of the paper, while the computational details for CGNSs 
% are given in Appendix G.
back_subjective_CIR = zeros(length(eps_ord_values), back_plot_len);
back_approx_objective_CIR = zeros(1, back_plot_len);

% Complete backward CIR relative-entropy metric δ^b(t;T) (see Sections
% 2.3.4-2.3.5 of the paper):
%     δ^b(t;T) := |δ(t,T') - δ(0,T')| as T' → T^-,
%               = |P(p^{upd}(x(t)|target(s<=T); conditioning), p^{lag}(x(t)|target(s<=T-Δt); conditioning))
%                  - P(p^{upd}(x(0)|target(s<=T); conditioning), p^{lag}(x(0)|target(s<=T-Δt); conditioning))| as Δt → 0^+,
% where the non-target wave variable has been conditionally masked in the
% posterior distributions under the conditional ACI framework. As discussed in 
% Appendix G of the paper, in practice, this quantity is approximated in 
% discrete time by the adaptive-lag online smoother update from t_{N-1} to t_N=T 
% for each arbitrary observational time T in the simulated interval. This metric 
% is used to calculate the approximate objective backward CIR via a time 
% integral over the natural time t instead of integrating the associated 
% subjective backward CIR over ε as in the definition to get the exact objective 
% backward CIR length (as well as to calculate the subjective backward CIR 
% itself). Using the notation from Appendix G of the paper, back_RE_metric below 
% stores the uncentred increments P^j_T,
%     P^j_T ≈ lim_{T'→T^-} δ(t_j,T'),
% over the retained natural-time window (shortened by lag_bound), approximating
% the limit of complete information incorporation lim_{T'→T^-} δ(t_j,T') for
% small enough Δt (observational time-step difference between target(t_{N-1})
% and target(t_N), assumed uniform throughout for simplicity). Centring
% back_RE_metric against its earliest retained lag gives complete_back_RE_metric
% (based on lag_bound), we get:
%     |P^j_T - P^{earliest retained j closer to t=0}_T| ≈ |P^j_T - P^0_T|
%                                                       ≈ lim_{T'→T^-} δ^b(t_j;T')
%                                                       = lim_{Δt→0^+} δ^b(t_j;T-Δt).
% The columns of the 1D complete_back_RE_metric vector correspond to the natural 
% time t indexed by t_j (i.e. j∈{T-lag_bound,...,T}) and calculated at each 
% observational time T=t_N in the simulated interval, iterated over in the loop 
% below. This lag-truncated approximation to |P^j_T - P^0_T| approaches the 
% exact value at each observational time T whenever the retained natural-time 
% window still reaches t=0 based on the lag_bound. Note that this complete
% backward CIR metric vector is calculated regardless of whether
% calculate_and_plot_normalised_backward_CIR_metric = false, since its values 
% are required for computing the subjective backward CIRs, the definition-based 
% (exact) objective backward CIR, and its computationally efficient objective 
% approximation.
if calculate_and_plot_normalised_backward_CIR_metric
    % The optional normalised complete backward CIR metric heatmap cache retains 
    % only the requested displayed backward lag rows, with the rows indexed by 
    % the natural-time-centred, current-observational-time-based backward lag 
    % T-t, induced by the moving observation time T=t_N in the plotting interval 
    % (T∈[time_start_plot,time_end_plot]) and the retained natural time t=t_j 
    % (j∈{T-lag_bound,...,T} or j∈{first_idx,first_idx+1,...,last_idx}, whatever 
    % is smaller, where first_idx and last_idx are determined by the plotting 
    % interval), and the columns indexed by the moving observational time t_N=T 
    % in the plotting window. This approach avoids storing a full memory-heavy 
    % O(N²) CIR metric cache; specifically, this retains only the portion of the 
    % complete backward CIR metric calculated at each observational time T=t_N 
    % in the loop below, under complete_back_RE_metric, that is necessary for 
    % plotting the normalised heatmap, leaving its calculated values unchanged 
    % and omitted lag rows as NaN. The computations of the displayed CIR metric 
    % values are unchanged by this convention. Cache columns preserve the 
    % original available backward-lag bound while also retaining only the 
    % amount of user-requested lag rows for optimal display.
    back_lag_steps = min(round(back_lag_time_end_plot/dt)+1, min(lag_bound+1, back_plot_len));
    % Not to be confused with the complete_back_RE_metric variable, which is the 
    % 1D vector of the complete backward CIR metric values at each observational 
    % time T=t_N in the simulated interval, iterated over in the loop below; see 
    % the comments above for more details. Essentially, this cache stores the
    % complete_back_RE_metric values at each observational time T=t_N in the 
    % plotting interval and only for the retained backward lag rows chosen by 
    % the user through back_lag_time_end_plot, with the omitted lag rows set 
    % to NaN.
    back_complete_RE_metric = nan(back_lag_steps, plot_len);
end
back_max_RE_metric = zeros(1, back_plot_len);

% Used for the text-based progress bar.
start_time = tic;

for n = 3:N+1

    % Text-based progress bar.
    progress_bar('Adaptive-Lag Online Smoother Algorithm and Calculation of the Backward CIRs', n-2, N-1, start_time);

    online_adapt_mean{n}(n) = filter_mean(n);
    online_adapt_cov{n}(n) = filter_cov(n);

    dx_obs = observed_state(:, n) - observed_state(:, n-1);
    aux_vec = filter_mean(n-1) ...
              - E_j_matrices(n-1)*((1 + L_y*dt)*filter_mean(n-1) + f_y(n-1)*dt) ...
              + F_j_matrices(:, :, n-1)*(dx_obs - (L_x(:, n-1)*filter_mean(n-1) + f_x(:, n-1))*dt);
    online_adapt_mean{n}(n-1) = E_j_matrices(n-1)*filter_mean(n) + aux_vec;
    aux_mat = filter_cov(n-1) ...
              - E_j_matrices(n-1)*(1 + L_y*dt)*filter_cov(n-1) ...
              - F_j_matrices(:, :, n-1)*L_x(:, n-1)*filter_cov(n-1)*dt;
    online_adapt_cov{n}(n-1) = E_j_matrices(n-1)*filter_cov(n)*E_j_matrices(n-1)' + aux_mat;

    online_mean_innovation = online_adapt_mean{n}(n-1) - filter_mean(n-1);
    online_cov_innovation = online_adapt_cov{n}(n-1) - filter_cov(n-1);

    first_lagged_idx = max((n-1)-lag_bound, 1);
    lag_window_len = (n-1) - first_lagged_idx + 1;
    back_RE_metric = zeros(1, lag_window_len);
    for t = (n-1):-1:first_lagged_idx
        if t == n-1
            update_matrices_adapt{n-2}(n-1) = 1;
        elseif t == n-2
            update_matrices_adapt{n-2}(n-2) = E_j_matrices(n-2);
        else
            update_matrices_adapt{n-2}(t) = update_matrices_adapt{n-3}(t)*E_j_matrices(n-2);
        end

        % Calculating the information gain criterion for the adaptive-lag online 
        % smoother update from t_{N-1} to t_N=T for each retained natural time 
        % t_j∈{t_{N-1}-lag_bound,...,t_{N-1}}.
        mean_diff = update_matrices_adapt{n-2}(t)*online_mean_innovation;
        cov_ratio_temp = 1 ... 
                         + (update_matrices_adapt{n-2}(t)*online_cov_innovation*update_matrices_adapt{n-2}(t)')/online_adapt_cov{n-1}(t);
        signal_temp = 0.5*mean_diff'*(online_adapt_cov{n-1}(t) \ mean_diff);
        if ~isreal(cov_ratio_temp) || ~isfinite(cov_ratio_temp) || cov_ratio_temp <= 0
            error('FBCIR:InvalidCovarianceRatio', 'The covariance ratio for calculating the complete backward CIR metric and lengths is invalid at natural time index %d and observational time index %d.', t, n)
        end
        dispersion_temp = 0.5*(cov_ratio_temp - log(cov_ratio_temp) - 1);
        lag_local_idx = t - first_lagged_idx + 1;
        back_RE_metric(lag_local_idx) = signal_temp + dispersion_temp;
    end

    % Complete backward CIR metric approximation at observation time t_N=T. The
    % subtraction by the earliest retained lag approximates the t=0 baseline
    % information deficit in the shortened-lag implementation above; see
    % Appendix G of the paper and the original online smoother paper
    % (https://doi.org/10.1007/s00332-026-10271-x) for more information.
    complete_back_RE_metric = abs(back_RE_metric - back_RE_metric(1));
    max_back_metric = max(complete_back_RE_metric);

    % If the maximum is essentially zero, corresponding to a non-causal
    % relationship at this time, then do not calculate the CIR and set
    % equal to 0 instead (to avoid operationally inflated CIR values
    % due to the normalisation in the objective CIR approximation and
    % to numerical precision errors).
    if max_back_metric > RE_metric_threshold
        % Efficient objective backward CIR approximation. It is an upper-bound
        % approximation in general and becomes exact under the monotonicity
        % condition for the complete backward CIR metric described in the paper.

        % Estimation of the complete backward CIR metric's integral using the
        % L1 grid-function norm from Appendix G of the paper since we deal with
        % discrete-time quantities with the adaptive-lag CGNS online smoother.
        adaptive_approx_objective_backward_CIR_n = sum(complete_back_RE_metric)*dt/max_back_metric;

        % Smooth quadrature: Estimation of the complete backward CIR metric's
        % integral using a composite trapezoidal rule.
        % adaptive_approx_objective_backward_CIR_n = trapz(complete_back_RE_metric)*dt/max_back_metric;

        % Smoother quadrature: Estimation of the complete backward CIR metric's
        % integral using a composite Simpson's 1/3 rule.
        % adaptive_approx_objective_backward_CIR_n = simps(complete_back_RE_metric)*dt/max_back_metric;

        % Number of steps to update back in time for the adaptive-lag CGNS
        % online smoother due to the new target effect observation 
        % target(t_N) = target(T).
        adaptive_lag_steps = min(ceil(adaptive_approx_objective_backward_CIR_n/dt), n-1);
    else
        adaptive_approx_objective_backward_CIR_n = 0;
        adaptive_lag_steps = 0;
    end

    if (first_idx <= n) && (n <= target_last_idx)
        % Map the global observation index n (corresponding to an arbitrary
        % t_N = T) to the plotting-window column and store the complete backward
        % CIR maximum and efficient objective backward CIR approximation for
        % normalisation and plotting later on.
        back_plot_idx = n - first_idx + 1;
        back_max_RE_metric(back_plot_idx) = max_back_metric;
        back_approx_objective_CIR(back_plot_idx) = adaptive_approx_objective_backward_CIR_n;

        % Calculation of the subjective backward CIR lengths for all ε values.
        % The centred backward CIR metric is ordered from the earliest retained
        % natural time up to T-Δt = T-dt; the first threshold exceedance
        % determines how far back from T the attributable interval extends and,
        % by extension, which state estimation updates in the adaptive-lag
        % online smoother possess significant information gain due to observing
        % target(t_N) = target(T). The prefix cumulative maximum makes each
        % threshold crossing equivalent to find(complete_back_RE_metric > ε, 1, 'first') 
        % without repeatedly scanning the full retained lag window.
        
        % Indexing convention: complete_back_RE_metric(end) corresponds to T-Δt,
        % the lagged online smoother endpoint used to approximate the limit of
        % complete information, T'→T^-, through the adaptive-lag online smoother
        % update from T-Δt to T due to the newly obtained target observation
        % target(t_N) = target(T). Therefore, subj_CIR_idx counts grid intervals 
        % from the first ε-threshold exceedance to [0,T-Δt]. A 
        % terminal-observation ([0,t_N=T]) convention could be used instead with 
        % (subj_CIR_idx+1)*dt, such that we measure the lag directly to T; 
        % still, we keep the T-Δt/grid-function convention because Appendix G of 
        % the paper defines the complete backward CIR metric through this lagged 
        % online smoother update and uses an L1 grid-function norm-based sum for 
        % the computationally efficient objective backward CIR approximation due 
        % to the discrete-time formulation adopted in the code.
        prefix_max_back_RE = cummax(complete_back_RE_metric);
        subj_CIR_idx = zeros(length(eps_ord_values), 1);
        first_exceedance_idx = 1;
        for eps_idx = length(eps_ord_values):-1:1
            while (first_exceedance_idx <= length(prefix_max_back_RE)) && ...
                  (prefix_max_back_RE(first_exceedance_idx) <= CIR_epsilon_values(eps_idx))
                first_exceedance_idx = first_exceedance_idx + 1;
            end

            if first_exceedance_idx <= length(prefix_max_back_RE)
                subj_CIR_idx(eps_idx) = length(prefix_max_back_RE) - first_exceedance_idx;
            end
        end
        back_subjective_CIR(:, back_plot_idx) = subj_CIR_idx*dt;

        if calculate_and_plot_normalised_backward_CIR_metric
            % Store the normalised complete backward CIR metric's numerator on
            % a backward-lag grid for the optional heatmap over (T,T-t).
            natural_idx = first_lagged_idx:(n-1);
            plot_nat_mask = (natural_idx >= first_idx) & (natural_idx <= target_last_idx);
            lag_plot_rows = n - natural_idx(plot_nat_mask);
            metric_values = complete_back_RE_metric(plot_nat_mask);
            lag_plot_mask = lag_plot_rows <= back_lag_steps;
            back_complete_RE_metric(lag_plot_rows(lag_plot_mask), back_plot_idx) = metric_values(lag_plot_mask);
        end
    end

    for j = (n-1):-1:1
        if (1 <= j) && (j <= n-1-adaptive_lag_steps)
            online_adapt_mean{n}(j) = online_adapt_mean{n-1}(j);
            online_adapt_cov{n}(j) = online_adapt_cov{n-1}(j);
        elseif (n-adaptive_lag_steps <= j) && (j <= n-1)
            online_adapt_mean{n}(j) = online_adapt_mean{n-1}(j) + update_matrices_adapt{n-2}(j)*online_mean_innovation;
            online_adapt_cov{n}(j) = online_adapt_cov{n-1}(j) + update_matrices_adapt{n-2}(j)*online_cov_innovation*update_matrices_adapt{n-2}(j);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  NOTE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% As in the forward case, the exact objective backward CIR is calculated
% through the definition by averaging the corresponding subjective backward CIR
% over ε. Here, both the definition-based value and the efficient approximation
% are computed for the same lag-truncated backward CIR metric. Here, by "exact",
% we mean "definition-based" as opposed to "approximation-based" (the latter 
% being the computationally efficient objective backward CIR approximation), 
% since numerical quadrature with a truncated ε interval is still used for the
% integral, which is sufficiently accurate for a large enough ε resolution and 
% coverage.

% When choosing a sufficiently large epsilon_resolution and a sufficiently large 
% k where 10⁻ᵏ ≤ ε, then the objective CIR at each time tⱼ can be calculated 
% through the definition by averaging the corresponding subjective CIR over ε 
% via a numerical quadrature method using the following command (the flip 
% operations are needed to put the ε (see eps_ord_values variable) interval in 
% ascending order):

% Estimation of the integral using a composite trapezoidal rule (a mask using 
% RE_metric_threshold, like the one utilised below and in the approximate 
% objective CIR calculation, can be considered to avoid operationally inflated 
% CIR values in the following due to the normalisation by the maxima).
% back_defn_objective_CIR = trapz(10.^flip(eps_ord_values), flipud(back_subjective_CIR), 1)./back_max_RE_metric;

% Estimation of the integral using a composite Simpson's 1/3 rule for better
% accuracy. The same RE_metric_threshold convention used in the approximate 
% objective CIR calculation is applied below before normalisation of this 
% definition-based/exact length by the maxima, so as to avoid operationally 
% inflated CIR values.
back_RE_metric_mask = back_max_RE_metric > RE_metric_threshold;
if calculate_and_plot_exact_objective_backward_CIR
    back_defn_objective_CIR = simps(10.^flip(eps_ord_values), flipud(back_subjective_CIR), 1);
    back_defn_objective_CIR(~back_RE_metric_mask) = 0;
    back_defn_objective_CIR(back_RE_metric_mask) = back_defn_objective_CIR(back_RE_metric_mask)./back_max_RE_metric(back_RE_metric_mask);
end

% The adaptive-lag online smoother history is no longer needed after all 
% backward CIR metric and objective/subjective CIR-length calculations are 
% complete. Clear it before plotting to reduce memory use.
clear online_adapt_mean online_adapt_cov update_matrices_adapt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%  PLOTTING ACI AND CIR ANALYSIS RESULTS  %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Time series of the ACI metric, heatmap of the subjective forward CIR length
% over time and ε, and time series of the computationally efficient approximate 
% objective forward CIR length for the selected Lorenz-84 conditional causal 
% relationship.
figure('WindowState', 'maximized');

subplot(3, 1, 1)
plot(plot_time, ACI_metric, 'k', LineWidth=1.5)
title(sprintf('ACI Metric for %s', selected_forward_relationship_title))
xlabel('t')
ylabel('Relative Entropy')
box on
grid on
fontsize(16, 'points')

subplot(3, 1, 2)
xx = plot_time;
yy = eps_ord_values;
[X, Y] = meshgrid(xx, yy);
pcolor(X, Y, log10(forw_subjective_CIR(:, 1:plot_len)));
shading interp
colormap("jet")
clim([-2.5, 1.5])
cb = colorbar('eastoutside', Ticks=-2.5:0.5:1.5, TickLabels=arrayfun(@(x) sprintf('10^{%.1f}', x), -2.5:0.5:1.5, 'UniformOutput', false));
cb.Position(1) = cb.Position(1) + 0.1;
cb.Position(3) = 0.015;
ylim([lowest_order, highest_order])
set(gca, 'YDir','reverse')
title(sprintf('Subjective Forward CIR Length for %s (Logarithmic Scale)', selected_forward_relationship_title))
xlabel('t (Natural Time)')
ylabel('ε (Relative Entropy Threshold)')
yticks(epsilon_axis_tick_values)
yticklabels(arrayfun(@(x) sprintf('10^{%g}', x), epsilon_axis_tick_values, 'UniformOutput', false))
fontsize(16, 'points')

subplot(3, 1, 3)
plot(plot_time, forw_approx_objective_CIR(1:plot_len), 'Color', [0.850 0.325 0.098], LineWidth=1.5)
xlim([time_start_plot, time_end_plot])
title(sprintf('Approximate Objective Forward CIR Length for %s', selected_forward_relationship_title))
xlabel('t')
ylabel('Forward CIR Length')
set(gca,'XGrid', 'on', 'YGrid', 'off')
box on
grid on
fontsize(16, 'points')

% Comparison of the time series of the exact objective forward CIR, calculated 
% via the ε-average definition of the subjective forward CIR, and its 
% computationally efficient lower-bounding approximation.
if calculate_and_plot_exact_objective_forward_CIR

figure('WindowState', 'maximized');
plot(plot_time, forw_approx_objective_CIR(1:plot_len), 'Color', [0.850 0.325 0.098], LineWidth=1.5)
hold on
plot(plot_time, forw_defn_objective_CIR, 'k--', LineWidth=1)
title(sprintf('Comparison of the Objective Forward CIR Length Algorithms for %s: Definition vs Approximation', selected_forward_relationship_title))
xlabel('t')
ylabel('Objective Forward CIR Length')
legend('Efficient Lower-Bound Approximation', 'Definition', NumColumns=2)
box on
grid on
fontsize(16, 'points')

end

% Heatmap of the normalised forward CIR metric for the selected Lorenz-84 
% conditional causal relationship:
%     δ^f(T';t)/max{δ^f(T';t): T'∈[t,T]},
% over natural time t and lagged observational time after t, T'-t. This is the
% CIR metric from Section 2.3.1 of the paper with t fixed and T' varying:
%     δ^f(T';t) = δ(t,T')
% normalised by its maximum over T'∈[t,T].
if calculate_and_plot_normalised_forward_CIR_metric

figure('WindowState', 'maximized');

% The plotting on the y-axis is done in terms of the lagged observational time 
% after t, T'-t, and in the smaller current-natural-time-centred lagged 
% observational time window defined by lag_obs_time_end_plot for optimal 
% plotting.
forw_metric_max = max_forw_RE_metric(1:plot_len);
% We apply the same RE_metric_threshold mask convention before normalising the 
% optional CIR metric heatmap values, as it was used for calculating the 
% objective CIRs while avoiding operationally inflated values.
forw_metric_max(~forw_RE_metric_mask) = NaN;
forw_data = forw_RE_metric;
clear forw_RE_metric
forw_data = forw_data./forw_metric_max;
% To avoid script errors during this optional heatmap plotting, we set the 
% non-finite normalised CIR metric values to NaN and, since we plot on a 
% logarithmic scale, we take the log10 of the absolute value and real part of 
% the normalised CIR metric values. Nevertheless, for transparency and as a 
% diagnostic summary, a brief report of the display-only numerical adjustments 
% made to the normalised CIR metric values for the optional logarithmic heatmap 
% plotting, i.e. number of non-real, non-finite, and non-positive values, etc., 
% is printed to the console for the user's information. The computed CIR metrics 
% and lengths remain unchanged.
shared_helper('report_heatmap_display_diagnostic', forw_data, plot_time, 'Forward')
forw_data = log10(abs(real(forw_data)));
forw_data(~isfinite(forw_data)) = NaN;
h = imagesc([time_start_plot time_end_plot], [0, (forw_lag_steps-1)*dt], forw_data(1:forw_lag_steps,:));
set(h, 'AlphaData', ~isnan(forw_data(1:forw_lag_steps,:)))
colormap("jet")
clim([-10, 0])
cb = colorbar('eastoutside', Ticks=-10:0, TickLabels=arrayfun(@(x) sprintf('10^{%.1f}', x), -10:0, 'UniformOutput', false));
cb.Position(1) = cb.Position(1) + 0.1;
cb.Position(3) = 0.015;
set(gca, 'Color', [0.7, 0.7, 0.7])
set(gca, 'YDir','reverse')
ylim([0, (forw_lag_steps-1)*dt])
title({['Normalised Forward CIR Metric for ', selected_forward_relationship_title, ':'], ...
       '\delta^f(T'';t)/max_{T''\in[t,T]}\{\delta^f(T'';t)\}, \delta^f(T'';t):=\delta(t,T'') (Logarithmic Scale)'}, ...
      'Interpreter', 'tex')
xlabel('t (Natural Time)')
ylabel("T'-t (Lagged Observational Time After t)")
fontsize(16, 'points')

end

% Time series of the ACI metric, heatmap of the subjective backward CIR length
% over time and ε, and time series of the computationally efficient approximate 
% objective backward CIR length for the selected Lorenz-84 conditional causal 
% relationship.
figure('WindowState', 'maximized');

subplot(3, 1, 1)
plot(plot_time, ACI_metric, 'k', LineWidth=1.5)
title(sprintf('ACI Metric for %s', selected_forward_relationship_title))
xlabel('t')
ylabel('Relative Entropy')
box on
grid on
fontsize(16, 'points')

subplot(3, 1, 2)
xx = plot_time;
yy = eps_ord_values;
[X, Y] = meshgrid(xx, yy);
pcolor(X, Y, log10(back_subjective_CIR));
shading interp
colormap("jet")
clim([-2.5, 1.5])
cb = colorbar('eastoutside', Ticks=-2.5:0.5:1.5, TickLabels=arrayfun(@(x) sprintf('10^{%.1f}', x), -2.5:0.5:1.5, 'UniformOutput', false));
cb.Position(1) = cb.Position(1) + 0.1;
cb.Position(3) = 0.015;
ylim([lowest_order, highest_order])
set(gca, 'YDir','reverse')
title(sprintf('Subjective Backward CIR Length for %s (Logarithmic Scale)', selected_backward_relationship_title))
xlabel('T (Observational Time)')
ylabel('ε (Relative Entropy Threshold)')
yticks(epsilon_axis_tick_values)
yticklabels(arrayfun(@(x) sprintf('10^{%g}', x), epsilon_axis_tick_values, 'UniformOutput', false))
fontsize(16, 'points')

subplot(3, 1, 3)
plot(plot_time, back_approx_objective_CIR, 'Color', [0 0.247 0.541], LineWidth=1.5)
xlim([time_start_plot, time_end_plot])
title(sprintf('Approximate Objective Backward CIR Length for %s', selected_backward_relationship_title))
xlabel('T')
ylabel('Backward CIR Length')
set(gca,'XGrid', 'on', 'YGrid', 'off')
box on
grid on
fontsize(16, 'points')

% Comparison of the time series of the exact objective backward CIR, calculated
% via the ε-average definition of the subjective backward CIR, and its
% computationally efficient upper-bounding approximation.
if calculate_and_plot_exact_objective_backward_CIR

figure('WindowState', 'maximized');
plot(plot_time, back_approx_objective_CIR, 'Color', [0 0.247 0.541], LineWidth=1.5)
hold on
plot(plot_time, back_defn_objective_CIR, 'k--', LineWidth=1)
title(sprintf('Comparison of the Objective Backward CIR Length Algorithms for %s: Definition vs Approximation', selected_backward_relationship_title))
xlabel('T')
ylabel('Objective Backward CIR Length')
legend('Efficient Upper-Bound Approximation', 'Definition', NumColumns=2)
box on
grid on
fontsize(16, 'points')

end

% Heatmap of the normalised complete backward CIR metric for the selected 
% Lorenz-84 conditional causal relationship:
%     lim_{T'→T^-} δ^b(t;T')/max_t{lim_{T'→T^-} δ^b(t;T')},
% over the current online-smoother observation time T and backward lag T-t. Here
% T denotes the moving observation time t_N, not the terminal simulation time.
% This is the centred complete backward form of the CIR metric from Section 
% 2.3.1:
%     lim_{T'→T^-} δ^b(t;T') := lim_{T'→T^-} |δ(t,T') - δ(0,T')|,
% normalised by its maximum over the retained natural-time window. Note that 
% this t-maximum in the normalised complete backward CIR metric is taken over 
% the retained natural-time window (shortened by lag_bound), while still 
% considering the plotting window of [time_start_plot,time_end_plot], i.e.
%     t ∈ [T-lag_bound,T] 
%         = [min{time_start_plot, time_end_plot-lag_bound},time_end_plot]
%         = [{earliest retained t closer to t=0},time_end_plot] ⊆ [time_start_plot,time_end_plot].
if calculate_and_plot_normalised_backward_CIR_metric

figure('WindowState', 'maximized');

% The plotting on the y-axis is done in terms of the backward lag, T-t, and in 
% the smaller natural-time-centred, current-observational-time-based backward
% lag window defined by back_lag_time_end_plot for optimal plotting.
back_metric_max = back_max_RE_metric;
% We apply the same RE_metric_threshold mask convention before normalising the 
% optional CIR metric heatmap values, as it was used for calculating the 
% objective CIRs while avoiding operationally inflated values.
back_metric_max(~back_RE_metric_mask) = NaN;
back_data = back_complete_RE_metric;
clear back_complete_RE_metric
back_data = back_data./back_metric_max;
% To avoid script errors during this optional heatmap plotting, we set the 
% non-finite normalised CIR metric values to NaN and, since we plot on a 
% logarithmic scale, we take the log10 of the absolute value and real part of 
% the normalised CIR metric values. Nevertheless, for transparency and as a 
% diagnostic summary, a brief report of the display-only numerical adjustments 
% made to the normalised CIR metric values for the optional logarithmic heatmap 
% plotting, i.e. number of non-real, non-finite, and non-positive values, etc., 
% is printed to the console for the user's information. The computed CIR metrics 
% and lengths remain unchanged.
shared_helper('report_heatmap_display_diagnostic', back_data, plot_time, 'Complete Backward')
back_data = log10(abs(real(back_data)));
back_data(~isfinite(back_data)) = NaN;
h = imagesc([time_start_plot time_end_plot], [0, (back_lag_steps-1)*dt], back_data(1:back_lag_steps,:));
set(h, 'AlphaData', ~isnan(back_data(1:back_lag_steps,:)))
colormap("jet")
clim([-10, 0])
cb = colorbar('eastoutside', Ticks=-10:0, TickLabels=arrayfun(@(x) sprintf('10^{%.1f}', x), -10:0, 'UniformOutput', false));
cb.Position(1) = cb.Position(1) + 0.1;
cb.Position(3) = 0.015;
set(gca, 'Color', [0.7, 0.7, 0.7])
set(gca, 'YDir','reverse')
xlim([time_start_plot, time_end_plot])
ylim([0, (back_lag_steps-1)*dt])
title({['Normalised Complete Backward CIR Metric for ', selected_backward_relationship_title, ':'], ...
       'lim_{T''\rightarrow T^-}\delta^b(t;T'')/max_{t\in[0,T]}\{lim_{T''\rightarrow T^-}\delta^b(t;T'')\},', ...
       'lim_{T''\rightarrow T^-}\delta^b(t;T''):=lim_{T''\rightarrow T^-}|\delta(t,T'')-\delta(0,T'')| (Logarithmic Scale)'}, ...
      'Interpreter', 'tex')
xlabel('T (Observational Time)')
ylabel('T-t (Backward Lag Before T)')
fontsize(16, 'points')

end