From 24392f2e35a84ae62dd64d146087038604e8b0b2 Mon Sep 17 00:00:00 2001 From: Karthik Date: Fri, 8 Aug 2025 00:15:36 +0200 Subject: [PATCH] Corrected code to more accurately depict the distribution of the order parameter --- .../simulateDistribution.m | 332 +++++++++--------- 1 file changed, 171 insertions(+), 161 deletions(-) diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m index 87d40a6..c584b66 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m @@ -1,174 +1,166 @@ -%% Evolve to skewed to single gaussian with lower mean +%% Evolve across a second-order-like transition clear; clc; N_params = 50; -N_reps = 500; +N_reps = 50; 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 +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.8; % low final mean +mu_start = 1.2; % high initial mean +mu_end = 0.2; % low final mean -sigma_start = 0.2; % wide initial std -sigma_end = 0.07; % narrow final std +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 +max_skew = 5; % peak skew strength -% Loop through alpha for i = 1:N_params - alpha = alpha_values(i); + 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]); + % === 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 - data = skewnormrnd(mu, sigma, skew_strength, N_reps); + sigma = sigma_end; % narrow at end end - all_data{i} = data; + % === 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); - kurt_vals(i) = kappa(4); + kappa4_vals(i) = kappa(4); kappa5_vals(i) = kappa(5); kappa6_vals(i) = kappa(6); end -%% Evolve to bimodal +%% Evolve across a first-order-like transition +% First-order-like distribution evolution with significant bimodality clear; clc; -N_params = 50; -N_reps = 500; -alpha_values = linspace(0, 45, N_params); +N_params = 50; +N_reps = 50; +alpha_values = linspace(0, 45, N_params); -all_data = cell(1, N_params); +all_data = cell(1, N_params); -bimodal_start = 20; % start earlier -transition_width = 5; % wider window for smoothness +% 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 < (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))); - + 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 - % After transition: bimodal, but not strongly balanced - w1 = 0.5; - w2 = 0.5; - - 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))); + % Region IV: Final stable low-mean Gaussian + data = normrnd(mu_low, sigma_initial, [N_reps, 1]); 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); + % 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.5, 1.8, 200); % max[g²] values on y-axis +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 @@ -194,9 +186,45 @@ 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 @@ -208,59 +236,41 @@ 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); +nexttile; plot(scan_vals, mean_vals, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Mean of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); +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 -subplot(3,2,2); +nexttile; plot(scan_vals, var_vals, 's-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Variance of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); +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 -subplot(3,2,3); +nexttile; plot(scan_vals, skew_vals, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Skewness of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); +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. Kurtosis -subplot(3,2,4); -plot(scan_vals, kurt_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); -title('Fourth-order Cumulant of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); +% 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; -% 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');