%% Evolve across a second-order-like transition
clear; clc;

N_params     = 50;
N_reps       = 50;
alpha_values = linspace(0, 45, N_params);
all_data     = cell(1, N_params);

% Transition control
alpha_start       = 5;     % where sigma starts changing
alpha_widen_end   = 15;    % when sigma finishes first change
alpha_shift_start = 15;    % when mean starts shifting
alpha_end         = 40;    % when mean finishes shifting and sigma narrows

mu_start    = 1.2;   % high initial mean
mu_end      = 0.2;   % low final mean

sigma_start = 0.25;  % wide std at start
sigma_mid   = 0.15;  % mid-range std in middle
sigma_end   = 0.07;  % narrow std at end

max_skew    = 5;     % peak skew strength

for i = 1:N_params
    alpha = alpha_values(i);

    % === Sigma evolution (variance large -> small) ===
    if alpha < alpha_start
        sigma = sigma_start;  % wide at start
    elseif alpha < alpha_widen_end
        % Smooth transition from wide to mid
        t_sigma = (alpha - alpha_start) / (alpha_widen_end - alpha_start);
        sigma = sigma_start * (1 - t_sigma) + sigma_mid * t_sigma;
    elseif alpha < alpha_end
        % Smooth transition from mid to narrow
        t_sigma = (alpha - alpha_widen_end) / (alpha_end - alpha_widen_end);
        sigma = sigma_mid * (1 - t_sigma) + sigma_end * t_sigma;
    else
        sigma = sigma_end;  % narrow at end
    end

    % === Mean evolution ===
    if alpha < alpha_shift_start
        mu = mu_start;  % fixed at high initially
    elseif alpha <= alpha_end
        % Smooth cosine shift
        t_mu = (alpha - alpha_shift_start) / (alpha_end - alpha_shift_start);
        smooth_t_mu = (1 - cos(pi * t_mu)) / 2;
        mu = mu_start * (1 - smooth_t_mu) + mu_end * smooth_t_mu;
    else
        mu = mu_end;
    end

    % === Skew evolution ===
    if alpha < alpha_end
        t_skew = (alpha - alpha_start) / (alpha_end - alpha_start);
        skew_strength = max_skew * (1 - t_skew);  % fade out
    else
        skew_strength = 0;
    end

    % Generate data
    if abs(skew_strength) < 1e-2
        data = normrnd(mu, sigma, [N_reps, 1]);
    else
        data = skewnormrnd(mu, sigma, skew_strength, N_reps);
    end

    all_data{i} = data;

    % Cumulants
    kappa = computeCumulants(data, 6);
    mean_vals(i)   = kappa(1);
    var_vals(i)    = kappa(2);
    skew_vals(i)   = kappa(3);
    kappa4_vals(i) = kappa(4);
    kappa5_vals(i) = kappa(5);
    kappa6_vals(i) = kappa(6);
end

%% Evolve across a first-order-like transition
% First-order-like distribution evolution with significant bimodality
clear; clc;

N_params       = 50;
N_reps         = 50;
alpha_values   = linspace(0, 45, N_params);

all_data       = cell(1, N_params);

% Define transition regions
skewed_start       = 10;
bimodal_start      = 20;
bimodal_end        = 35;
final_narrow_start = 40;

% Peak positions and widths
mu_high      = 1.2;   % Initial metastable peak
mu_low       = 0.2;   % Final stable peak
mu_new_peak  = 0.8;   % New peak appears slightly lower
sigma_initial = 0.08;

for i = 1:N_params
    alpha = alpha_values(i);

    if alpha < skewed_start
        % Region I: Narrow unimodal at high mean
        data = normrnd(mu_high, sigma_initial, [N_reps, 1]);

    elseif alpha < bimodal_start
        % Region II: Slightly skewed
        t_skew = (alpha - skewed_start) / (bimodal_start - skewed_start);
        mu     = mu_high - 0.15 * t_skew;
        sigma  = sigma_initial + 0.02 * t_skew;
        skew_strength = 3 * t_skew;
        data   = skewnormrnd(mu, sigma, skew_strength, N_reps);

    elseif alpha < bimodal_end
        % Region III: Bimodal with fixed or slowly drifting peak positions
        t = (alpha - bimodal_start) / (bimodal_end - bimodal_start);  % t in [0, 1]

        % Increased separation between peaks
        drift_amount = 0.3;  % larger = more drift toward final mean
        sep_offset   = 0.25; % larger = more initial separation between peaks

        % Peaks start separated and move toward mu_low
        mu1 = mu_high * (1 - t)^drift_amount + mu_low * (1 - (1 - t)^drift_amount);        % Right peak drifts to left
        mu2 = (mu_new_peak - sep_offset) * (1 - t)^drift_amount + mu_low * (1 - (1 - t)^drift_amount);  % Left peak moves slightly

        sigma1 = sigma_initial + 0.02 * (1 - abs(0.5 - t) * 2);
        sigma2 = sigma1;

        % Weight shift: right peak dies out, left peak grows
        w2 = 0.5 + 0.5 * t;  % left peak grows: 0.5 → 1
        w1 = 1 - w2;         % right peak fades: 0.5 → 0

        N1 = round(N_reps * w1);
        N2 = N_reps - N1;

        mode1 = normrnd(mu1, sigma1, [N1, 1]);
        mode2 = normrnd(mu2, sigma2, [N2, 1]);

        data = [mode1; mode2];
        data = data(randperm(length(data)));

    else
        % Region IV: Final stable low-mean Gaussian
        data = normrnd(mu_low, sigma_initial, [N_reps, 1]);
    end

    % Store data and compute cumulants
    all_data{i}    = data;
    kappa          = computeCumulants(data, 6);
    mean_vals(i)   = kappa(1);
    var_vals(i)    = kappa(2);
    skew_vals(i)   = kappa(3);
    kappa4_vals(i) = kappa(4);
    kappa5_vals(i) = kappa(5);
    kappa6_vals(i) = kappa(6);
end

%% === Compute 2D PDF heatmap: f(x, alpha) ===
x_grid = linspace(0.0, 1.8, 200);  % max[g²] values on y-axis
pdf_matrix = zeros(numel(x_grid), N_params);  % Now: rows = y, columns = alpha

for i = 1:N_params
    data = all_data{i};
    f = ksdensity(data, x_grid, 'Bandwidth', 0.025);
    pdf_matrix(:, i) = f;  % Transpose for y-axis to be vertical
end

% === Plot PDF vs. alpha heatmap ===
figure(2); clf;
set(gcf, 'Color', 'w', 'Position',[100 100 950 750])

imagesc(alpha_values, x_grid, pdf_matrix);
set(gca, 'YDir', 'normal');  % Flip y-axis to normal orientation
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$\mathrm{max}[g^{(2)}]$', 'Interpreter', 'latex', 'FontSize', 14);
title('Evolving PDF of $\mathrm{max}[g^{(2)}]$', ...
      'Interpreter', 'latex', 'FontSize', 16);

colormap(Colormaps.coolwarm());  % More aesthetic than default
colorbar;
c = colorbar;
ylabel(c, 'PDF', 'FontSize', 14, 'Interpreter', 'latex');
set(gca, 'FontSize', 14);

%% Animate evolving distribution and cumulant value
figure(1); clf;
set(gcf, 'Color', 'w', 'Position',[100 100 1300 750])

for i = 1:N_params
    clf;

    % PDF
    subplot(1,2,1); cla; hold on;
    data = all_data{i};
    
    % Plot histogram with normalized PDF
    histogram(data, 'Normalization', 'pdf', 'BinWidth', 0.03, ...
        'FaceColor', [0.3 0.5 0.8], 'EdgeColor', 'k', 'FaceAlpha', 0.6);
    
    title(sprintf('Histogram at $\\alpha = %.1f^\\circ$', alpha_values(i)), ...
        'Interpreter', 'latex', 'FontSize', 16);
    xlabel('$\mathrm{max}[g^{(2)}]$', 'Interpreter', 'latex', 'FontSize', 14);
    ylabel('PDF', 'FontSize', 14);
    set(gca, 'FontSize', 12); grid on;
    xlim([0.0, 2.0]);


    % Cumulant evolution
    subplot(1,2,2); hold on;
    plot(alpha_values(1:i), kappa4_vals(1:i), 'bo-', 'LineWidth', 2);
    title('Binder Cumulant Tracking', 'Interpreter', 'latex', 'FontSize', 16);
    xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', 14);
    ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', 14);
    xlim([0, 45]); grid on;
    set(gca, 'FontSize', 12);

    pause(0.3);
end

%% === Plotting ===
figure(1)
set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact');

scan_vals = alpha_values;  % your parameter sweep values

% Define font style for consistency
axis_fontsize  = 14;
label_fontsize = 16;
title_fontsize = 16;

% 1. Mean with error bars (if you have error data, else just plot)
% If no error, replace errorbar with plot or omit error data
% For now, no error bars assumed
nexttile;
plot(scan_vals, mean_vals, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Mean', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_1$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;

% 2. Variance
nexttile;
plot(scan_vals, var_vals, 's-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Variance', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_2$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;

% 3. Skewness
nexttile;
plot(scan_vals, skew_vals, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Skewness', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_3$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;

% 4. Binder Cumulant
nexttile;
plot(scan_vals, kappa4_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Binder Cumulant', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;

% Super title (you can customize the string)
sgtitle('Cumulants of a simulated evolving distribution', ...
    'FontWeight', 'bold', 'FontSize', 18, 'Interpreter', 'latex');

%% === Helper: Cumulant Calculation ===
function kappa = computeCumulants(data, max_order)
    data = data(:);
    mu = mean(data);
    c = zeros(1, max_order);
    centered = data - mu;
    for n = 1:max_order
        c(n) = mean(centered.^n);
    end
    kappa = zeros(1, max_order);
    kappa(1) = mu;
    kappa(2) = c(2);
    kappa(3) = c(3);
    kappa(4) = c(4) - 3*c(2)^2;
    kappa(5) = c(5) - 10*c(3)*c(2);
    kappa(6) = c(6) - 15*c(4)*c(2) - 10*c(3)^2 + 30*c(2)^3;
end

%% === Helper: Skewed Normal Distribution ===
function x = skewnormrnd(mu, sigma, alpha, n)
    % Skew-normal using Azzalini's method
    delta = alpha / sqrt(1 + alpha^2);
    u0    = randn(n,1);
    v     = randn(n,1);
    u1    = delta * u0 + sqrt(1 - delta^2) * v;
    x     = mu + sigma * u1 .* sign(u0);
end