multiscale_joint_y_1_y_2_cause.m

% Code corresponding to the multiscale atmospheric-variability model in Section
% 3.2 - "Impact of Resolving Non-Target Variables in the CIR Metrics of a
% Multiscale Model for Atmospheric Variability" 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 unconditional Assimilative Causal Inference 
% (ACI) and of the associated forward and backward Causal Influence Range (CIR) 
% metrics to a stochastic multiscale conceptual model for studying low-frequency 
% atmospheric variability:
%
%     dx_1/dt = a_1·x_1 - c_1·x_1^3 - x_2·(M + M_1·x_1 + M_2·x_2)
%               + I_11·x_1·y_1 + I_12·x_1·y_2 + L_11·y_1 + L_12·y_2 + f_1^x(t)
%               + σ_x1·dot(W_x1) + Σ_{k=1,2} (σ_yk/γ_k)·(L_1k - I_1k·x_1)·dot(W_yk),
%     dx_2/dt = -c_2·x_2 + x_1·(M + M_1·x_1 + M_2·x_2)
%               + I_21·x_2·y_1 + I_22·x_2·y_2 + L_21·y_1 + L_22·y_2 + f_2^x(t)
%               + σ_x2·dot(W_x2) + Σ_{k=1,2} (σ_yk/γ_k)·(L_2k - I_2k·x_2)·dot(W_yk),
%     dy_1/dt = -(γ_1/ε)·y_1 - L_11·x_1 - L_21·x_2
%               + N·y_2 - I_11·x_1^2 - I_21·x_2^2 + f_1^y(t)
%               + (σ_y1/sqrt(ε))·dot(W_y1),
%     dy_2/dt = -(γ_2/ε)·y_2 - L_12·x_1 - L_22·x_2
%               - N·y_1 - I_12·x_1^2 - I_22·x_2^2 + f_2^y(t)
%               + (σ_y2/sqrt(ε))·dot(W_y2).
%
% This is a reduced-order model for large-scale geophysical flow, often derived
% via a stochastic-mode reduction procedure on general circulation models; see
% works on the Majda--Timofeyev--Vanden-Eijnden (MTV) method and this paper on 
% normal forms for reduced stochastic climate models: 
%   ➤ https://doi.org/10.1073/pnas.0900173106
% It is also a conditional Gaussian nonlinear system (CGNS) for the posterior 
% (y_1,y_2) | (x_1,x_2); see Appendix G of the paper. It governs the nonlinear, 
% coupled evolution between the resolved slow atmospheric/climate/environmental
% variables x = (x_1,x_2)ᵀ and the unresolved fast weather components 
% y = (y_1,y_2)ᵀ. Quadratic nonlinearities preserve energy and produce 
% intermittent instabilities and non-Gaussian marginal/joint PDFs, while the 
% linear-operator dynamics in x encode multiscale structures typical of 
% turbulent flows (skew-symmetric: Coriolis effects, topographic Rossby wave 
% propagation; negative-definite: Stochastic stability via dissipative processes 
% like surface drag), as well as oscillatory feedbacks in y. The small parameter 
% 0 < ε <= 1 separates the slow and fast timescales. The system also includes 
% scale-dependent time-varying forcings, with large-scale feedbacks representing 
% slow external inputs such as decadal oscillations, while small-scale forcings 
% capture faster intraseasonal variabilities and/or interannual phase-locking. 
% Finally, correlated additive and multiplicative (CAM) noise is added to the 
% system; the multiplicative noise stems from the nonlinear advection of the 
% resolved scales by the unresolved scales, whereas the additive noise part 
% stems from the linear operator obtained by linearising the equations of motion 
% around the climatological background state. This CAM noise is a common feature 
% added to conceptual models for efficient and effective description of
% non-Gaussian low-frequency climate variability in the presence of unresolved
% processes (e.g. heavy-tailed PDFs, skewness, and intermittency). 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 unconditional, joint causal relationship (y_1(t),y_2(t)) → (x_1,x_2) 
% over time t∈[0,T]. Therefore, the weather vector y is treated as the latent
% candidate cause, while the atmospheric vector x is the observed effect. Unlike
% the other scripts in this codebase, which employ conditional ACI, no 
% observational channel is resolved or masked in this implementation, with both 
% x_1 and x_2 contributing to the Bayesian update of y(t). This is the 
% unconditional ACI baseline used in the paper before conditioning on the 
% ancillary/non-target variables to resolve and study the marginal links or 
% contributions in the ACI and CIRs, such as the relationships 
% y_1 → x_1 | (x_2,y_2) and y_2 → x_2 | (x_1,y_1).
%
% In the paper, the bivariate link (y_1(t),y_2(t)) → (x_1,x_2) is used to 
% quantify the aggregate or joint weather-to-atmosphere causal influence before
% "decomposing" it to its marginal constituents through the conditional ACI and 
% CIR metrics of:
%     y_1 → x_1 | (x_2,y_2),
%     y_2 → x_2 | (x_1,y_1).
% Note that, based on the chosen model parameter values from the paper, y_2 and 
% x_1 do not causally/dynamically interact with each other explicitly in the 
% system, but only do so through confounding/spurious associations via the other 
% variables (and likewise for y_1 and x_2); therefore, the causal properties of 
% the joint link (ACI and CIRs) are roughly decomposable to those of the two
% aforementioned conditional marginal relationships. Forward CIR analysis of 
% this joint relationship assesses how long the current fast weather state 
% remains predictively informative for future atmospheric evolution, while the 
% backward CIR attributes current atmospheric excursions or bursts to past 
% weather states. This script produces the analyses associated with Panels 
% (e)-(f) of Figure 7 of the paper.
%
% 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('multiscale_joint_y_1_y_2_cause.m');
%
% M-file Scripts:
%   ➤ multiscale_joint_y_1_y_2_cause.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(444)

% Total number of time steps within the given time interval; not to be confused
% with the feedback parameter N in the model equations, corresponding to the
% N_cross variable in the code, that controls the oscillatory cross-coupling 
% between the fast weather variables.
N = 20000;
% Numerical integration time step.
dt = 0.005;
% Total simulation time.
T = N*dt;
% Time grid (in natural time) for the time-dependent forcings.
time = (0:N)*dt;

% Main timescale separation parameter. The value used in the paper is ε = 0.1.
epsilon = 0.1;

% Plotting time window for the figures.
time_start_plot = 50;
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 = 5;
% 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 = 5;

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

% Dimension of the observed effect and unobserved candidate-cause vectors.
obs_dim = 2;
state_dim = 2;
% state_dim-dimensional identity matrix for quick reference in the algorithms.
eye_state = eye(state_dim);

% Observed slow atmospheric/effect variables: x = (x_1,x_2)ᵀ.
x = zeros(obs_dim, N+1);
% Unobserved fast weather/candidate-cause variables: y = (y_1,y_2)ᵀ.
y = zeros(state_dim, N+1);

% Slow atmospheric self-interaction and energy-conserving feedback parameters.
a_1 = 1;
c_1 = 1/3;
M = 0.5;
M_1 = 0.5;
c_2 = 0.4;
M_2 = -1.5;

% Fast weather damping and oscillatory cross-coupling parameters.
gamma_1 = 0.5;
N_cross = 4;
gamma_2 = 1.2;

% Quadratic interactions-operator feedback matrix I.
I = zeros(2, 2);
I(1, 1) = 0.6;
I(1, 2) = 0;
I(2, 1) = 0;
I(2, 2) = 2;

% Linear-operator feedback matrix L.
L = zeros(2, 2);
L(1, 1) = 1;
L(1, 2) = 0;
L(2, 1) = 0;
L(2, 2) = 1.5;

% Large-scale forcing in x_1.
fx_1 = 0;
% Large-scale forcing in x_2.
fx2_mean = 4;
fx2_amplitude = 0.5;
fx2_period = 36;
fx_2 = zeros(1, N+1);
fx_2(1) = fx2_mean + fx2_amplitude*sin(2*pi*time(1)/fx2_period);

% Small-scale weather forcings in y_1 and y_2.
fy_1 = 1;
fy_2 = -1;

% Additive/state-independent noise/diffusion amplitudes in x and y.
sigma_x1 = 0.15;
sigma_x2 = 0.3;
sigma_y1 = 1;
sigma_y2 = 2;

% CGNS representation with y latent and x observed; the "^x" and "^y"
% superscripts used in the comments denote the observable and unobservable
% processes, respectively, and are used solely for consistency with the notation
% in Appendix G of the paper.
%     dx = [L^x(t,x)·y + f^x(t,x)]dt + S_1^x(t,x)·dW_x + S_2^x(t,x)·dW_y,
%     dy = [L^y(t,x)·y + f^y(t,x)]dt + S_1^y(t,x)·dW_x + S_2^y(t,x)·dW_y.

% Observable coefficient matrix: Feedback of y in x.
L_x = zeros(obs_dim, state_dim, N+1);
% Forcing in the observable process.
f_x = zeros(obs_dim, N+1);
% Noise feedback matrix corresponding to the independent x-noise channel in the 
% observable process x, i.e. the (1,1)-diagonal block of the full diffusion 
% matrix of the joint process (x,y).
Sx_1 = [
    sigma_x1, 0;
    0, sigma_x2
];
% Noise feedback matrix corresponding to the y-noise channel in the observable 
% process x, producing the CAM noise in the (1,2) and (2,1) off-diagonal blocks 
% of the full diffusion matrix of the joint process (x,y).
Sx_2 = zeros(obs_dim, state_dim, N+1);

% Unobservable coefficient matrix: Feedback of y in y.
L_y = [
    -gamma_1/epsilon, N_cross;
    -N_cross, -gamma_2/epsilon
];
% Forcing in the unobservable process.
f_y = zeros(state_dim, N+1);
% Noise feedback matrix corresponding to the x-noise channel in the unobservable 
% process y, producing the CAM noise in the (1,2) and (2,1) off-diagonal blocks 
% of the full diffusion matrix of the joint process (x,y).
Sy_1 = zeros(state_dim, obs_dim);
% Noise feedback matrix corresponding to the independent y-noise channel in the 
% unobservable process y, i.e. the (2,2)-diagonal block of the full diffusion 
% matrix of the joint process (x,y).
Sy_2 = [
    sigma_y1/sqrt(epsilon), 0;
    0, sigma_y2/sqrt(epsilon)
];

% Initial conditions for x = (x_1,x_2)ᵀ and y = (y_1,y_2)ᵀ; zero initial 
% conditions are also a valid choice for this model.
x(:, 1) = [
    0;
    fx_2(1)/c_2
];
y(:, 1) = [
    fy_1/(gamma_1/epsilon);
    fy_2/(gamma_2/epsilon)
];

% Initiating the time-dependent model components.
L_x(:, :, 1) = [
    I(1, 1)*x(1, 1) + L(1, 1), I(1, 2)*x(1, 1) + L(1, 2);
    I(2, 1)*x(2, 1) + L(2, 1), I(2, 2)*x(2, 1) + L(2, 2)
];
f_x(:, 1) = [
    a_1*x(1, 1) - c_1*x(1, 1)^3 - x(2, 1)*(M + M_1*x(1, 1) + M_2*x(2, 1)) + fx_1;
    -c_2*x(2, 1) + x(1, 1)*(M + M_1*x(1, 1) + M_2*x(2, 1)) + fx_2(1)
];
Sx_2(:, :, 1) = [
    sigma_y1*(L(1, 1) - I(1, 1)*x(1, 1))/gamma_1, sigma_y2*(L(1, 2) - I(1, 2)*x(1, 1))/gamma_2;
    sigma_y1*(L(2, 1) - I(2, 1)*x(2, 1))/gamma_1, sigma_y2*(L(2, 2) - I(2, 2)*x(2, 1))/gamma_2
];
f_y(:, 1) = [
    -L(1, 1)*x(1, 1) - L(2, 1)*x(2, 1) - I(1, 1)*x(1, 1)^2 - I(2, 1)*x(2, 1)^2 + fy_1;
    -L(1, 2)*x(1, 1) - L(2, 2)*x(2, 1) - I(1, 2)*x(1, 1)^2 - I(2, 2)*x(2, 1)^2 + fy_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 Multiscale Model', j-1, N, start_time);

    % Wiener increments. For consistency with the other multiscale scripts in
    % terms of the random Gaussian noise sampling during simulation, note the
    % sampling order of x_1 -> x_2 -> y_1 -> y_2, which we conform to 
    % throughout. This is important for a consistent comparison of the ACI and
    % CIR metrics on the same set of simulated true signals and across the 
    % multiscale assimilative causal relationships studied in the paper. 
    dW_x = sqrt(dt)*randn(obs_dim, 1);
    dW_y = sqrt(dt)*randn(state_dim, 1);

    % Euler-Maruyama update of the state variables based on the dynamical model.
    x(:, j) = x(:, j-1) + (L_x(:, :, j-1)*y(:, j-1) + f_x(:, j-1))*dt ...
              + Sx_1*dW_x + Sx_2(:, :, j-1)*dW_y;
    y(:, j) = y(:, j-1) + (L_y*y(:, j-1) + f_y(:, j-1))*dt ...
              + Sy_1*dW_x + Sy_2*dW_y;

    % Updating the time-dependent model components.
    L_x(:, :, j) = [
        I(1, 1)*x(1, j) + L(1, 1), I(1, 2)*x(1, j) + L(1, 2);
        I(2, 1)*x(2, j) + L(2, 1), I(2, 2)*x(2, j) + L(2, 2)
    ];
    fx_2(j) = fx2_mean + fx2_amplitude*sin(2*pi*time(j)/fx2_period);
    f_x(:, j) = [
        a_1*x(1, j) - c_1*x(1, j)^3 - x(2, j)*(M + M_1*x(1, j) + M_2*x(2, j)) + fx_1;
        -c_2*x(2, j) + x(1, j)*(M + M_1*x(1, j) + M_2*x(2, j)) + fx_2(j)
    ];
    Sx_2(:, :, j) = [
        sigma_y1*(L(1, 1) - I(1, 1)*x(1, j))/gamma_1, sigma_y2*(L(1, 2) - I(1, 2)*x(1, j))/gamma_2;
        sigma_y1*(L(2, 1) - I(2, 1)*x(2, j))/gamma_1, sigma_y2*(L(2, 2) - I(2, 2)*x(2, j))/gamma_2
    ];
    f_y(:, j) = [
        -L(1, 1)*x(1, j) - L(2, 1)*x(2, j) - I(1, 1)*x(1, j)^2 - I(2, 1)*x(2, j)^2 + fy_1;
        -L(1, 2)*x(1, j) - L(2, 2)*x(2, j) - I(1, 2)*x(1, j)^2 - I(2, 2)*x(2, j)^2 + fy_2
    ];

end

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

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

% (x,y)-phase-space 2D projection plots; these plots are indicators of how the 
% multiscale, nonlinear interactions between x and y produce non-Gaussian 
% marginal and joint PDFs. 
figure('WindowState', 'maximized');

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

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

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

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

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

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

% Time series and time-averaged marginal PDFs of x_1, x_2, y_1, and y_2.
figure('WindowState', 'maximized');

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

subplot(4, 3, 3)
xx1 = linspace(min(x(1, plot_idx)), max(x(1, plot_idx)), 250);
tavg_pdf_x1 = ksdensity(x(1, plot_idx), xx1);
plot(xx1, tavg_pdf_x1, 'k', LineWidth=1.5)
hold on
mean_x1 = mean(x(1, plot_idx));
std_x1 = std(x(1, plot_idx));
plot(xx1, normpdf(xx1, mean_x1, std_x1), 'm--', LineWidth=1.5)
xlabel('x_1')
ylabel('p(x_1)')
title('PDF of x_1')
legend('Empirical', 'Gaussian Fit', Location='best')
box on
grid on
fontsize(16, 'points')

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

subplot(4, 3, 6)
xx2 = linspace(min(x(2, plot_idx)), max(x(2, plot_idx)), 250);
tavg_pdf_x2 = ksdensity(x(2, plot_idx), xx2);
plot(xx2, tavg_pdf_x2, 'k', LineWidth=1.5)
hold on
mean_x2 = mean(x(2, plot_idx));
std_x2 = std(x(2, plot_idx));
plot(xx2, normpdf(xx2, mean_x2, std_x2), 'm--', LineWidth=1.5)
xlabel('x_2')
ylabel('p(x_2)')
title('PDF of x_2')
box on
grid on
fontsize(16, 'points')

subplot(4, 3, [7, 8])
plot(plot_time, y(1, plot_idx), 'b', LineWidth=1.5)
ylabel('y_1')
title('Time Series of y_1')
box on
grid on
fontsize(16, 'points')

subplot(4, 3, 9)
yy1 = linspace(min(y(1, plot_idx)), max(y(1, plot_idx)), 250);
tavg_pdf_y1 = ksdensity(y(1, plot_idx), yy1);
plot(yy1, tavg_pdf_y1, 'b', LineWidth=1.5)
hold on
mean_y1 = mean(y(1, plot_idx));
std_y1 = std(y(1, plot_idx));
plot(yy1, normpdf(yy1, mean_y1, std_y1), 'm--', LineWidth=1.5)
xlabel('y_1')
ylabel('p(y_1)')
title('PDF of y_1')
box on
grid on
fontsize(16, 'points')

subplot(4, 3, [10, 11])
plot(plot_time, y(2, plot_idx), 'b', LineWidth=1.5)
xlabel('t')
ylabel('y_2')
title('Time Series of y_2')
box on
grid on
fontsize(16, 'points')

subplot(4, 3, 12)
yy2 = linspace(min(y(2, plot_idx)), max(y(2, plot_idx)), 250);
tavg_pdf_y2 = ksdensity(y(2, plot_idx), yy2);
plot(yy2, tavg_pdf_y2, 'b', LineWidth=1.5)
hold on
mean_y2 = mean(y(2, plot_idx));
std_y2 = std(y(2, plot_idx));
plot(yy2, normpdf(yy2, mean_y2, std_y2), 'm--', LineWidth=1.5)
xlabel('y_2')
ylabel('p(y_2)')
title('PDF of y_2')
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 vector y(t).
% Unlike the conditional ACI and CIR scripts, this unconditional joint-cause
% analysis uses the full precision matrix, so both x_1 and x_2 directly reduce
% the posterior uncertainty of y(t).
S_xoS_x_inv = zeros(obs_dim, obs_dim, N+1);
S_xoS_x_inv(:, :, 1) = inv(Sx_1*Sx_1' + Sx_2(:, :, 1)*Sx_2(:, :, 1)');

% Cross-Grammian (noise cross-interaction) between the latent and observed 
% process noise feedbacks.
S_yoS_x = zeros(state_dim, obs_dim, N+1);
S_yoS_x(:, :, 1) = Sy_1*Sx_1' + Sy_2*Sx_2(:, :, 1)';

% Grammian of the unobservable process noise feedback.
S_yoS_y = Sy_1*Sy_1' + Sy_2*Sy_2';

% Posterior filter mean of the latent vector y.
filter_mean = zeros(state_dim, N+1);
% Initial value of the posterior filter mean.
filter_mean(:, 1) = y(:, 1); mu0 = filter_mean(:, 1);
% Posterior filter covariance matrix of the latent vector y.
filter_cov = zeros(state_dim, state_dim, 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*eye_state; 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 = x(:, j) - x(:, j-1);
    filter_gain_aux = S_yoS_x(:, :, j-1) + 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(:, :, j-1)*(dx - (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(:, :, j-1)*filter_gain_aux')*dt;

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

    S_xoS_x_inv(:, :, j) = inv(Sx_1*Sx_1' + Sx_2(:, :, j)*Sx_2(:, :, j)');
    S_yoS_x(:, :, j) = Sy_1*Sx_1' + Sy_2*Sx_2(:, :, j)';

end

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

% Posterior smoother mean of the latent vector y.
smoother_mean = zeros(state_dim, N+1);
% Posterior smoother covariance matrix of the latent vector y.
smoother_cov = zeros(state_dim, state_dim, 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 (y_1(t),y_2(t)) → (x_1,x_2) 
% 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(state_dim, state_dim, N+1);
F_j_matrices = zeros(state_dim, obs_dim, N+1);
G_x_j = L_x(:, :, N+1) + S_yoS_x(:, :, N+1)'/filter_cov(:, :, N+1);
G_y_j = L_y + S_yoS_y/filter_cov(:, :, N+1);
C_jj = eye_state - G_y_j*dt;
H_j = filter_cov(:, :, N+1)\(L_y*filter_cov(:, :, N+1) + filter_cov(:, :, N+1)*L_y' + S_yoS_y);
K_j = S_xoS_x_inv(:, :, N+1)*G_x_j;
E_j_matrices(:, :, N+1) = C_jj + S_yoS_x(:, :, N+1)*K_j*dt;
F_j_matrices(:, :, N+1) = -filter_cov(:, :, N+1) ...
                          *(K_j' + (G_x_j'*K_j*filter_cov(:, :, N+1)*K_j' - filter_cov(:, :, N+1)\H_j'*filter_cov(:, :, N+1)*K_j' + L_y'*K_j')*dt ...
                            - L_x(:, :, N+1)'*(S_xoS_x_inv(:, :, N+1) + K_j*filter_cov(:, :, N+1)*K_j'*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_x_j = L_x(:, :, j) + S_yoS_x(:, :, j)'/filter_cov(:, :, j);
    G_y_j = L_y + S_yoS_y/filter_cov(:, :, j);
    C_jj = eye_state - G_y_j*dt;
    H_j = filter_cov(:, :, j)\(L_y*filter_cov(:, :, j) + filter_cov(:, :, j)*L_y' + S_yoS_y);
    K_j = S_xoS_x_inv(:, :, j)*G_x_j;
    E_j_matrices(:, :, j) = C_jj + S_yoS_x(:, :, j)*K_j*dt;
    F_j_matrices(:, :, j) = -filter_cov(:, :, j) ...
                             *(K_j' + (G_x_j'*K_j*filter_cov(:, :, j)*K_j' - filter_cov(:, :, j)\H_j'*filter_cov(:, :, j)*K_j' + L_y'*K_j')*dt ...
                             - L_x(:, :, j)'*(S_xoS_x_inv(:, :, j) + K_j*filter_cov(:, :, j)*K_j'*dt));

    % Auxiliary matrices for the posterior smoother mean and covariance updates.
    A_j = L_y - S_yoS_x(:, :, j)*S_xoS_x_inv(:, :, j)*L_x(:, :, j);
    B_j = S_yoS_y - S_yoS_x(:, :, j)*S_xoS_x_inv(:, :, j)*S_yoS_x(:, :, j)';
    dx = x(:, j+1) - x(:, j);

    % 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 ...
              + S_yoS_x(:, :, j)*S_xoS_x_inv(:, :, j)*(-dx + (L_x(:, :, j)*muT + f_x(:, j))*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 = x(:, j+1) - x(:, j);
    % E_j = E_j_matrices(:, :, j);
    % F_j = F_j_matrices(:, :, j);
    % mu = filter_mean(:, j) ...
    %      + E_j*(muT - (eye_state + L_y*dt)*filter_mean(:, j) - f_y(:, j)*dt) ...
    %      + F_j*(dx - (L_x(:, :, j)*filter_mean(:, j) + f_x(:, j))*dt);
    % R = filter_cov(:, :, j) ...
    %     + E_j*(RT*E_j' - (eye_state + 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 (unconditional) Gaussian 
% statistics of y_1 and y_2.
figure('WindowState', 'maximized');

subplot(2, 1, 1)
hold on
plot(plot_time, y(1, plot_idx), 'b', LineWidth=1.5, DisplayName='Truth')
plot(plot_time, filter_mean(1, plot_idx), 'g', LineWidth=1.5, DisplayName='Filter')
plot(plot_time, smoother_mean(1, plot_idx), 'r', LineWidth=1.5, DisplayName='Smoother')
filter_upper = filter_mean(1, plot_idx) + 2*real(sqrt(squeeze(filter_cov(1, 1, plot_idx))))';
filter_lower = filter_mean(1, plot_idx) - 2*real(sqrt(squeeze(filter_cov(1, 1, plot_idx))))';
smoother_upper = smoother_mean(1, plot_idx) + 2*real(sqrt(squeeze(smoother_cov(1, 1, plot_idx))))';
smoother_lower = smoother_mean(1, plot_idx) - 2*real(sqrt(squeeze(smoother_cov(1, 1, 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('y_1')
title('Time Series of y_1 and of its Filter and Smoother Posterior Statistics ((y_1(t),y_2(t)) \rightarrow (x_1,x_2))')
legend(Location='best')
box on
grid on
fontsize(16, 'points')

subplot(2, 1, 2)
hold on
plot(plot_time, y(2, plot_idx), 'b', LineWidth=1.5, DisplayName='Truth')
plot(plot_time, filter_mean(2, plot_idx), 'g', LineWidth=1.5, DisplayName='Filter')
plot(plot_time, smoother_mean(2, plot_idx), 'r', LineWidth=1.5, DisplayName='Smoother')
filter_upper = filter_mean(2, plot_idx) + 2*real(sqrt(squeeze(filter_cov(2, 2, plot_idx))))';
filter_lower = filter_mean(2, plot_idx) - 2*real(sqrt(squeeze(filter_cov(2, 2, plot_idx))))';
smoother_upper = smoother_mean(2, plot_idx) + 2*real(sqrt(squeeze(smoother_cov(2, 2, plot_idx))))';
smoother_lower = smoother_mean(2, plot_idx) - 2*real(sqrt(squeeze(smoother_cov(2, 2, 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('y_2')
title('Time Series of y_2 and of its Filter and Smoother Posterior Statistics ((y_1(t),y_2(t)) \rightarrow (x_1,x_2))')
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 y_1 and y_2 (unconditional filter and smoother); 
% normalised RMSE (NRMSE) and Pearson correlation coefficient.
nrmse_filter_y1 = rmse(y(1, plot_idx), filter_mean(1, plot_idx))/std(y(1, plot_idx));
nrmse_smoother_y1 = rmse(y(1, plot_idx), smoother_mean(1, plot_idx))/std(y(1, plot_idx));
relative_nrmse_change_y1 = 100*(nrmse_filter_y1 - nrmse_smoother_y1)/nrmse_filter_y1;
fprintf('Relative percentage change in smoother NRMSE from filter NRMSE for y_1 in (y_1(t),y_2(t)) \rightarrow (x_1,x_2) = %0.2f%% (+ = decreased, better)\n', relative_nrmse_change_y1);
nrmse_filter_y2 = rmse(y(2, plot_idx), filter_mean(2, plot_idx))/std(y(2, plot_idx));
nrmse_smoother_y2 = rmse(y(2, plot_idx), smoother_mean(2, plot_idx))/std(y(2, plot_idx));
relative_nrmse_change_y2 = 100*(nrmse_filter_y2 - nrmse_smoother_y2)/nrmse_filter_y2;
fprintf('Relative percentage change in smoother NRMSE from filter NRMSE for y_2 in (y_1(t),y_2(t)) \rightarrow (x_1,x_2) = %0.2f%% (+ = decreased, better)\n', relative_nrmse_change_y2);
corr_filter_y1_matrix = corrcoef(y(1, plot_idx), filter_mean(1, plot_idx));
corr_smoother_y1_matrix = corrcoef(y(1, plot_idx), smoother_mean(1, plot_idx));
corr_filter_y1 = corr_filter_y1_matrix(1, 2);
corr_smoother_y1 = corr_smoother_y1_matrix(1, 2);
relative_corr_change_y1 = corr_smoother_y1 - corr_filter_y1;
fprintf('Difference/Δ in smoother correlation from filter correlation for y_1 in (y_1(t),y_2(t)) \rightarrow (x_1,x_2) = %0.2f (+ = increased, better)\n', relative_corr_change_y1);
corr_filter_y2_matrix = corrcoef(y(2, plot_idx), filter_mean(2, plot_idx));
corr_smoother_y2_matrix = corrcoef(y(2, plot_idx), smoother_mean(2, plot_idx));
corr_filter_y2 = corr_filter_y2_matrix(1, 2);
corr_smoother_y2 = corr_smoother_y2_matrix(1, 2);
relative_corr_change_y2 = corr_smoother_y2 - corr_filter_y2;
fprintf('Difference/Δ in smoother correlation from filter correlation for y_2 in (y_1(t),y_2(t)) \rightarrow (x_1,x_2) = %0.2f (+ = increased, better)\n', relative_corr_change_y2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  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 (y_1(t),y_2(t)) → (x_1,x_2) at each time in the
% plotting window, i.e. the relative entropy from the posterior smoother to the
% filter of y(t) = (y_1(t),y_2(t))ᵀ, P(p_t^s(y|x), p_t^f(y|x)), which, due to 
% the Gaussianity of these posterior distributions, it is given by the 
% signal-dispersion formula; see Appendix G of the paper.
signal_smoother_filter = zeros(1, length(plot_idx));
dispersion_smoother_filter = zeros(1, length(plot_idx));
for idx = 1:length(plot_idx)
    j = plot_idx(idx);
    mean_diff = smoother_mean(:, j) - filter_mean(:, j);
    signal_smoother_filter(idx) = 0.5*mean_diff'*(filter_cov(:, :, j) \ mean_diff);
    cov_ratio_smoother_filter = smoother_cov(:, :, j)/filter_cov(:, :, j);
    % Alternative numerical safeguards due to the multidimensional latent causal
    % space: If covariance roundoff makes the log-det term unstable use
    % log(abs(det(...))) and max(real(...), 0) clamping on the final
    % signal-dispersion value. These are numerical stabilisers only, not part
    % of the mathematical Gaussian relative entropy formula used here.
    dispersion_smoother_filter(idx) = 0.5*(trace(cov_ratio_smoother_filter) - log(det(cov_ratio_smoother_filter)) - state_dim);
end
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(y(t_j)|x(s<=t_n)), are needed
% to evaluate the forward CIR metric (and by extension the subjective and
% objective forward CIRs) for (y_1(t),y_2(t)) → (x_1,x_2) 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(y(t_j)|x(s<=t_N)) to the lagged smoother distribution
% p(y(t_j)|x(s<=t_n)) 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(y(t_j)|x(s<=t_N)), p(y(t_j)|x(s<=t_n))),
% 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 unconditional online smoother distribution,
% p(y(t_j)|x(s<=t_n)), 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 effect observation x^n = x(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 y^j = y(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(state_dim, state_dim, n+1);
    online_fixed_mean{n} = zeros(state_dim, n);
    online_fixed_cov{n} = zeros(state_dim, state_dim, n);
end
for n = N:N+1
    online_fixed_mean{n} = zeros(state_dim, n);
    online_fixed_cov{n} = zeros(state_dim, state_dim, 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 = x(:, 2) - x(:, 1);
    aux_vec = filter_mean(:, 1) ...
              - E_j_matrices(:, :, 1)*((eye_state + L_y*dt)*filter_mean(:, 1) + f_y(:, 1)*dt) ...
              + F_j_matrices(:, :, 1)*(dx - (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)*(eye_state + 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 = x(:, n) - x(:, n-1);
        aux_vec = filter_mean(:, n-1) ...
                  - E_j_matrices(:, :, n-1)*((eye_state + L_y*dt)*filter_mean(:, n-1) + f_y(:, n-1)*dt) ...
                  + F_j_matrices(:, :, n-1)*(dx - (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)*(eye_state + 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) = eye_state;
            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 (y_1(t),y_2(t)) → (x_1,x_2) 
% 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(y(t)|x(s<=T)), p(y(t)|x(s<=T'))).
% 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);
        mean_diff = online_fixed_mean{end}(:, j) - online_fixed_mean{obs}(:, j);
        % Alternative numerical safeguards due to the multidimensional latent 
        % causal space: If covariance roundoff makes the log-det term unstable 
        % use log(abs(det(...))) and max(real(...), 0) clamping on the final
        % signal-dispersion value. These are numerical stabilisers only, not 
        % part of the mathematical Gaussian relative entropy formula used here.
        RE_n(obs-j+1) = 0.5*mean_diff'*(online_fixed_cov{obs}(:, :, j) \ mean_diff) ...
                        + 0.5*(trace(cov_ratio) - log(det(cov_ratio)) - state_dim);
    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}(y(t_j)|x(s<=t_N)),
% obtained from updating p^{lag}(y(t_j)|x(s<=t_{N-1})) due to the newly observed
% x(t_N) = x(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, are needed to 
% evaluate the (complete) backward CIR metric (and by extension the subjective 
% and objective backward CIRs) for (y_1,y_2) → (x_1(T),x_2(T)) 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}(y(t_j)|x(s<=t_N)) to the lagged online smoother
% distribution p^{lag}(y(t_j)|x(s<=t_{N-1})) is calculated at each observational 
% time t_N=T for all natural times 0<=t_j<=T,
%     P^j_T := P(p^{upd}(y(t_j)|x(s<=t_N)), p^{lag}(y(t_j)|x(s<=t_{N-1}))),
% where both posterior distributions are Gaussian for this CGNS; note how using
% the adaptive-lag online smoother distributions going from x(t_{N-1}) to x(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 x^N = x(t_N) = x(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 y^j = y(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(state_dim, state_dim, n+1);
    online_adapt_mean{n} = zeros(state_dim, n);
    online_adapt_cov{n} = zeros(state_dim, state_dim, n);
end
for n = N:N+1
    online_adapt_mean{n} = zeros(state_dim, n);
    online_adapt_cov{n} = zeros(state_dim, state_dim, 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 = x(:, 2) - x(:, 1);
aux_vec = filter_mean(:, 1) ...
          - E_j_matrices(:, :, 1)*((eye_state + L_y*dt)*filter_mean(:, 1) + f_y(:, 1)*dt) ...
          + F_j_matrices(:, :, 1)*(dx - (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)*(eye_state + 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 (y_1,y_2) → (x_1(T),x_2(T)) 
% at each observational time T in the plotting interval, for the same O(10⁻ᵏ) 
% orders of ε values as in the forward CIR implementations. 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}(y(t)|x(s<=T)), p^{lag}(y(t)|x(s<=T-Δt)))
%                  - P(p^{upd}(y(0)|x(s<=T)), p^{lag}(y(0)|x(s<=T-Δt)))| as Δt → 0^+.
% As discussed in Appendix G of the paper, in practice, this quantity is
% approximated in discrete time % 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 x(t_{N-1}) and 
% x(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 = x(:, n) - x(:, n-1);
    aux_vec = filter_mean(:, n-1) ...
              - E_j_matrices(:, :, n-1)*((eye_state + L_y*dt)*filter_mean(:, n-1) + f_y(:, n-1)*dt) ...
              + F_j_matrices(:, :, n-1)*(dx - (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)*(eye_state + 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) = eye_state;
        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 = eye_state ...
                         + (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);
        % Alternative numerical safeguards due to the multidimensional latent 
        % causal space: If covariance roundoff makes the log-det term unstable 
        % use log(abs(det(...))) and max(real(...), 0) clamping on the final
        % signal-dispersion value. These are numerical stabilisers only, not 
        % part of the mathematical Gaussian relative entropy formula used here.
        dispersion_temp = 0.5*(trace(cov_ratio_temp) - log(det(cov_ratio_temp)) - state_dim);
        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 effect observation x(t_N) = x(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
        % x(t_N) = x(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 observation 
        % x(t_N) = x(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 (y_1(t),y_2(t)) → (x_1,x_2).
figure('WindowState', 'maximized');

subplot(3, 1, 1)
plot(plot_time, ACI_metric, 'k', LineWidth=1.5)
title('ACI Metric for (y_1(t),y_2(t)) \rightarrow (x_1,x_2)')
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('Subjective Forward CIR Length for (y_1(t),y_2(t)) \rightarrow (x_1,x_2) (Logarithmic Scale)')
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('Approximate Objective Forward CIR Length for (y_1(t),y_2(t)) \rightarrow (x_1,x_2)')
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('Comparison of the Objective Forward CIR Length Algorithms for (y_1(t),y_2(t)) \rightarrow (x_1,x_2): Definition vs Approximation')
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 (y_1(t),y_2(t)) → (x_1,x_2):
%     δ^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 (y_1(t),y_2(t)) \rightarrow (x_1,x_2):', ...
       'δ^f(T'';t)/max_{T''\in[t,T]}\{δ^f(T'';t)\}, δ^f(T'';t):=δ(t,T'') (Logarithmic Scale)'})
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 (y_1,y_2) → (x_1(T),x_2(T)).
figure('WindowState', 'maximized');

subplot(3, 1, 1)
plot(plot_time, ACI_metric, 'k', LineWidth=1.5)
title('ACI Metric for (y_1(t),y_2(t)) \rightarrow (x_1,x_2)')
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('Subjective Backward CIR Length for (y_1,y_2) \rightarrow (x_1(T),x_2(T)) (Logarithmic Scale)')
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('Approximate Objective Backward CIR Length for (y_1,y_2) \rightarrow (x_1(T),x_2(T))')
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('Comparison of the Objective Backward CIR Length Algorithms for (y_1,y_2) \rightarrow (x_1(T),x_2(T)): Definition vs Approximation')
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
% (y_1,y_2) → (x_1(T),x_2(T)):
%     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 (y_1,y_2) \rightarrow (x_1(T),x_2(T)):', ...
       'lim_{T''\rightarrow T^-}δ^b(t;T'')/max_{t\in[0,T]}\{lim_{T''\rightarrow T^-}δ^b(t;T'')\},', ...
       'lim_{T''\rightarrow T^-}δ^b(t;T''):=lim_{T''\rightarrow T^-}|δ(t,T'')-δ(0,T'')| (Logarithmic Scale)'})
xlabel('T (Observational Time)')
ylabel('T-t (Backward Lag Before T)')
fontsize(16, 'points')

end