karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Corrected code to more accurately depict the distribution of the order parameter

Browse Source
This commit is contained in:
Karthik 2025-08-08 00:15:36 +02:00
parent 4a4ddf0e9c
commit 24392f2e35
1 changed files with 171 additions and 161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

266
Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/simulateDistribution.m
View File

@ -1,42 +1,66 @@
%% Evolve to skewed to single gaussian with lower mean %% Evolve across a second-order-like transition
clear; clc; clear; clc;
N_params = 50; N_params = 50;
N_reps = 500; N_reps = 50;
alpha_values = linspace(0, 45, N_params); alpha_values = linspace(0, 45, N_params);
all_data = cell(1, N_params); all_data = cell(1, N_params);
% Transition control % Transition control
alpha_start = 5; % where transition begins alpha_start = 5; % where sigma starts changing
alpha_end = 40; % where transition ends 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_start = 1.2; % high initial mean
mu_end = 0.8; % low final mean mu_end = 0.2; % low final mean
sigma_start = 0.2; % wide initial std sigma_start = 0.25; % wide std at start
sigma_end = 0.07; % narrow final std 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 for i = 1:N_params
alpha = alpha_values(i); alpha = alpha_values(i);
% Normalized transition variable t ∈ [0, 1] % === Sigma evolution (variance large -> small) ===
t = min(max((alpha - alpha_start) / (alpha_end - alpha_start), 0), 1); 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
% Use cosine-based smooth interpolation % === Mean evolution ===
smooth_t = (1 - cos(pi * t)) / 2; % ease-in-out 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
% Mean and sigma interpolation (smoothly decrease) % === Skew evolution ===
mu = mu_start * (1 - smooth_t) + mu_end * smooth_t; if alpha < alpha_end
sigma = sigma_start * (1 - smooth_t) + sigma_end * smooth_t; t_skew = (alpha - alpha_start) / (alpha_end - alpha_start);
skew_strength = max_skew * (1 - t_skew); % fade out
% Skewness: sinusoidal profile, max at middle else
skew_strength = max_skew * sin(t * pi); skew_strength = 0;
end
% Generate data % Generate data
if abs(skew_strength) < 1e-2 % near-zero skew, use normal if abs(skew_strength) < 1e-2
data = normrnd(mu, sigma, [N_reps, 1]); data = normrnd(mu, sigma, [N_reps, 1]);
else else
data = skewnormrnd(mu, sigma, skew_strength, N_reps); data = skewnormrnd(mu, sigma, skew_strength, N_reps);
@ -49,126 +73,94 @@ for i = 1:N_params
mean_vals(i) = kappa(1); mean_vals(i) = kappa(1);
var_vals(i) = kappa(2); var_vals(i) = kappa(2);
skew_vals(i) = kappa(3); skew_vals(i) = kappa(3);
kurt_vals(i) = kappa(4); kappa4_vals(i) = kappa(4);
kappa5_vals(i) = kappa(5); kappa5_vals(i) = kappa(5);
kappa6_vals(i) = kappa(6); kappa6_vals(i) = kappa(6);
end end
%% Evolve to bimodal %% Evolve across a first-order-like transition
% First-order-like distribution evolution with significant bimodality
clear; clc; clear; clc;
N_params = 50; N_params = 50;
N_reps = 500; N_reps = 50;
alpha_values = linspace(0, 45, N_params); alpha_values = linspace(0, 45, N_params);
all_data = cell(1, N_params); all_data = cell(1, N_params);
bimodal_start = 20; % start earlier % Define transition regions
transition_width = 5; % wider window for smoothness 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 for i = 1:N_params
alpha = alpha_values(i); alpha = alpha_values(i);
if alpha < (bimodal_start - transition_width) if alpha < skewed_start
% Pure skewed unimodal, smaller max skewness for subtlety % Region I: Narrow unimodal at high mean
skew_strength = 3 * (alpha / bimodal_start); data = normrnd(mu_high, sigma_initial, [N_reps, 1]);
data = skewnormrnd(1, 0.1, skew_strength, N_reps);
elseif alpha <= (bimodal_start + transition_width) elseif alpha < bimodal_start
% Smooth transition window % Region II: Slightly skewed
t = (alpha - (bimodal_start - transition_width)) / (2 * transition_width); 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);
% Weights transition 1 -> 0.7 elseif alpha < bimodal_end
w1 = 1 - 0.3 * t; % Region III: Bimodal with fixed or slowly drifting peak positions
w2 = 1 - w1; t = (alpha - bimodal_start) / (bimodal_end - bimodal_start); % t in [0, 1]
% Peak separation smaller (max delta = 0.3) % Increased separation between peaks
delta_max = 0.3; drift_amount = 0.3; % larger = more drift toward final mean
delta = delta_max * t; sep_offset = 0.25; % larger = more initial separation between peaks
mu1 = 1 - delta; % Peaks start separated and move toward mu_low
mu2 = 1 + delta; 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 = 0.1; sigma1 = sigma_initial + 0.02 * (1 - abs(0.5 - t) * 2);
sigma2 = 0.1; 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); N1 = round(N_reps * w1);
N2 = N_reps - N1; N2 = N_reps - N1;
mode1_samples = normrnd(mu1, sigma1, [N1, 1]); mode1 = normrnd(mu1, sigma1, [N1, 1]);
mode2_samples = normrnd(mu2, sigma2, [N2, 1]); mode2 = normrnd(mu2, sigma2, [N2, 1]);
data = [mode1_samples; mode2_samples]; data = [mode1; mode2];
data = data(randperm(length(data))); data = data(randperm(length(data)));
else else
% After transition: bimodal, but not strongly balanced % Region IV: Final stable low-mean Gaussian
w1 = 0.5; data = normrnd(mu_low, sigma_initial, [N_reps, 1]);
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)));
end end
% Store data and compute cumulants
all_data{i} = data; all_data{i} = data;
kappa = computeCumulants(data,6); kappa = computeCumulants(data, 6);
mean_vals(i) = kappa(1); mean_vals(i) = kappa(1);
var_vals(i) = kappa(2); var_vals(i) = kappa(2);
skew_vals(i) = kappa(3); skew_vals(i) = kappa(3);
kurt_vals(i) = kappa(4); kappa4_vals(i) = kappa(4);
kappa5_vals(i) = kappa(5); kappa5_vals(i) = kappa(5);
kappa6_vals(i) = kappa(6); kappa6_vals(i) = kappa(6);
end 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) === %% === 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 pdf_matrix = zeros(numel(x_grid), N_params); % Now: rows = y, columns = alpha
for i = 1:N_params for i = 1:N_params
@ -194,9 +186,45 @@ c = colorbar;
ylabel(c, 'PDF', 'FontSize', 14, 'Interpreter', 'latex'); ylabel(c, 'PDF', 'FontSize', 14, 'Interpreter', 'latex');
set(gca, 'FontSize', 14); 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 === %% === Plotting ===
figure(1) figure(1)
set(gcf, 'Color', 'w', 'Position', [100 100 950 750]) 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 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) % 1. Mean with error bars (if you have error data, else just plot)
% If no error, replace errorbar with plot or omit error data % If no error, replace errorbar with plot or omit error data
% For now, no error bars assumed % For now, no error bars assumed
subplot(3,2,1); nexttile;
plot(scan_vals, mean_vals, 'o-', 'LineWidth', 1.5, 'MarkerSize', 6); 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); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_1$', 'Interpreter', 'latex', 'FontSize', label_fontsize); ylabel('$\kappa_1$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize); set(gca, 'FontSize', axis_fontsize);
grid on; grid on;
% 2. Variance % 2. Variance
subplot(3,2,2); nexttile;
plot(scan_vals, var_vals, 's-', 'LineWidth', 1.5, 'MarkerSize', 6); 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); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_2$', 'Interpreter', 'latex', 'FontSize', label_fontsize); ylabel('$\kappa_2$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize); set(gca, 'FontSize', axis_fontsize);
grid on; grid on;
% 3. Skewness % 3. Skewness
subplot(3,2,3); nexttile;
plot(scan_vals, skew_vals, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6); 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); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_3$', 'Interpreter', 'latex', 'FontSize', label_fontsize); ylabel('$\kappa_3$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize); set(gca, 'FontSize', axis_fontsize);
grid on; grid on;
% 4. Kurtosis % 4. Binder Cumulant
subplot(3,2,4); nexttile;
plot(scan_vals, kurt_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6); plot(scan_vals, kappa4_vals, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
title('Fourth-order Cumulant of Distribution', 'FontSize', title_fontsize, 'Interpreter', 'latex'); title('Binder Cumulant', 'FontSize', title_fontsize, 'Interpreter', 'latex');
xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize); xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', 'FontSize', label_fontsize);
ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', label_fontsize); ylabel('$\kappa_4$', 'Interpreter', 'latex', 'FontSize', label_fontsize);
set(gca, 'FontSize', axis_fontsize); set(gca, 'FontSize', axis_fontsize);
grid on; 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) % Super title (you can customize the string)
sgtitle('Cumulants of a simulated evolving distribution', ... sgtitle('Cumulants of a simulated evolving distribution', ...
'FontWeight', 'bold', 'FontSize', 18, 'Interpreter', 'latex'); 'FontWeight', 'bold', 'FontSize', 18, 'Interpreter', 'latex');