%% Evolve to skewed to single gaussian with lower mean
clear; clc;

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

% Transition control
alpha_start = 5;     % where transition begins
alpha_end   = 40;    % where transition ends

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

sigma_start = 0.2;   % wide initial std
sigma_end   = 0.07;  % narrow final std

max_skew    = 5;     % peak skew strength

% Loop through alpha
for i = 1:N_params
    alpha         = alpha_values(i);

    % Normalized transition variable t ∈ [0, 1]
    t             = min(max((alpha - alpha_start) / (alpha_end - alpha_start), 0), 1);

    % Use cosine-based smooth interpolation
    smooth_t      = (1 - cos(pi * t)) / 2;  % ease-in-out

    % Mean and sigma interpolation (smoothly decrease)
    mu            = mu_start    * (1 - smooth_t) + mu_end    * smooth_t;
    sigma         = sigma_start * (1 - smooth_t) + sigma_end * smooth_t;

    % Skewness: sinusoidal profile, max at middle
    skew_strength = max_skew * sin(t * pi);

    % Generate data
    if abs(skew_strength) < 1e-2  % near-zero skew, use normal
        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);
    kurt_vals(i)   = kappa(4);
    kappa5_vals(i) = kappa(5);
    kappa6_vals(i) = kappa(6);
end

%% Evolve to bimodal
clear; clc;

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

all_data              = cell(1, N_params);

bimodal_start         = 20;     % start earlier
transition_width      = 5;      % wider window for smoothness

for i = 1:N_params
    alpha             = alpha_values(i);
    
    if alpha < (bimodal_start - transition_width)
        % Pure skewed unimodal, smaller max skewness for subtlety
        skew_strength = 3 * (alpha / bimodal_start);
        data          = skewnormrnd(1, 0.1, skew_strength, N_reps);
        
    elseif alpha <= (bimodal_start + transition_width)
        % Smooth transition window
        t             = (alpha - (bimodal_start - transition_width)) / (2 * transition_width);
        
        
        % Weights transition 1 -> 0.7
        w1            = 1 - 0.3 * t;
        w2            = 1 - w1;
        
        % Peak separation smaller (max delta = 0.3)
        delta_max     = 0.3;
        delta         = delta_max * t;
        
        mu1           = 1 - delta;
        mu2           = 1 + delta;
        
        sigma1        = 0.1;
        sigma2        = 0.1;
        
        N1            = round(N_reps * w1);
        N2            = N_reps - N1;
        
        mode1_samples = normrnd(mu1, sigma1, [N1, 1]);
        mode2_samples = normrnd(mu2, sigma2, [N2, 1]);
        
        data          = [mode1_samples; mode2_samples];
        data          = data(randperm(length(data)));
        
    else
        % After transition: bimodal, but not strongly balanced
        w1            = 0.7;
        w2            = 0.3;
        
        mu1           = 1 - 0.3;
        mu2           = 1 + 0.3;
        
        sigma1        = 0.1;
        sigma2        = 0.1;
        
        N1            = round(N_reps * w1);
        N2            = N_reps - N1;
        
        mode1_samples = normrnd(mu1, sigma1, [N1, 1]);
        mode2_samples = normrnd(mu2, sigma2, [N2, 1]);
        
        data          = [mode1_samples; mode2_samples];
        data          = data(randperm(length(data)));
    end
    
    all_data{i}    = data;
    kappa          = computeCumulants(data,6);
    mean_vals(i)   = kappa(1);
    var_vals(i)    = kappa(2);
    skew_vals(i)   = kappa(3);
    kurt_vals(i)   = kappa(4);
    kappa5_vals(i) = kappa(5);
    kappa6_vals(i) = kappa(6);
end

%% 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.5, 1.8]);


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

    pause(0.3);
end

%% === Compute 2D PDF heatmap: f(x, alpha) ===
x_grid = linspace(0.5, 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);

%% === Plotting ===
figure(1)
set(gcf, 'Color', 'w', 'Position', [100 100 950 750])

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
subplot(3,2,1);
plot(scan_vals, mean_vals, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Mean of Distribution', '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
subplot(3,2,2);
plot(scan_vals, var_vals, 's-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Variance of Distribution', '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
subplot(3,2,3);
plot(scan_vals, skew_vals, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Skewness of Distribution', '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. Kurtosis
subplot(3,2,4);
plot(scan_vals, kurt_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Kurtosis of Distribution', '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;

% 5. 5th-order cumulant
subplot(3,2,5);
plot(scan_vals, kappa5_vals, 'v-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Fifth-order Cumulant of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_5$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize);
grid on;

% 6. 6th-order cumulant
subplot(3,2,6);
plot(scan_vals, kappa6_vals, '>-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Sixth-order Cumulant of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_6$', '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