New and modified fitting and plotting routines (Comparison scripts will be added later).
This commit is contained in:
parent
36d57e5973
commit
dd4145ccd2
@ -34,7 +34,7 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
|
||||
Ncurves = numel(S_theta_all);
|
||||
fitResults = struct('pFit',[],'thetaFit',[],'xFit',[],'yFit',[],...
|
||||
'thetaFine',[],'isValid',[],'fitMaxTheta',[]);
|
||||
'thetaFine',[],'isValid',[], 'residual',[], 'fitMaxTheta',[]);
|
||||
|
||||
rawCurves = struct('x', cell(1, Ncurves), 'theta', cell(1, Ncurves));
|
||||
|
||||
@ -69,9 +69,8 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
S_theta_all_shifted = S_theta_all;
|
||||
theta_recentered = theta_vals;
|
||||
end
|
||||
|
||||
|
||||
for k = 1:Ncurves
|
||||
|
||||
x = S_theta_all_shifted{k};
|
||||
theta = theta_recentered;
|
||||
|
||||
@ -144,25 +143,27 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
|
||||
try
|
||||
pFit = lsqcurvefit(twoGauss, p0, theta, x, lb, ub, optsLSQ);
|
||||
|
||||
% --- Compute RMS of fit vs raw curve in [0, pi/4] ---
|
||||
maskRMS = theta >= 0 & theta <= pi/4;
|
||||
if any(maskRMS)
|
||||
x_raw_subset = x(maskRMS);
|
||||
y_fit_subset = twoGauss(pFit, theta(maskRMS));
|
||||
residual_rms_val = sqrt(mean((x_raw_subset - y_fit_subset).^2));
|
||||
end
|
||||
|
||||
% --- Handle missing or weak secondary peak, discard poor fits ---
|
||||
% --- Handle missing or weak secondary peak ---
|
||||
if pFit(3) <= opts.AmplitudeThreshold
|
||||
% amplitude too small → discard secondary
|
||||
pFit(3) = 0;
|
||||
pFit(4:5) = NaN;
|
||||
else
|
||||
% --- Compute RMS of fit vs raw curve in [0, pi/4] ---
|
||||
maskRMS = theta >= 0 & theta <= pi/4;
|
||||
if any(maskRMS)
|
||||
x_raw_subset = x(maskRMS);
|
||||
y_fit_subset = twoGauss(pFit, theta(maskRMS));
|
||||
rms_val = sqrt(mean((x_raw_subset - y_fit_subset).^2));
|
||||
if rms_val > opts.ResidualThreshold
|
||||
% residual too large → discard secondary
|
||||
pFit(3) = NaN;
|
||||
pFit(4:5) = NaN;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% --- Discard poor fits ---
|
||||
if residual_rms_val > opts.ResidualThreshold
|
||||
% residual too large → discard secondary
|
||||
pFit(3) = NaN;
|
||||
pFit(4:5) = NaN;
|
||||
end
|
||||
|
||||
% --- Map to canonical format with primary amplitude = 1 ---
|
||||
@ -179,6 +180,7 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
fitResults(k).thetaFine = thetaFine;
|
||||
fitResults(k).yFit = yFit;
|
||||
fitResults(k).isValid = true;
|
||||
fitResults(k).residual = residual_rms_val;
|
||||
fitResults(k).fitMaxTheta = fitMaxTheta;
|
||||
continue; % success → next curve
|
||||
end
|
||||
@ -200,14 +202,31 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
optsLSQfree = optimoptions('lsqcurvefit', 'Display', 'off', ...
|
||||
'MaxFunctionEvaluations', 2e4, 'MaxIterations', 2e4);
|
||||
|
||||
pFree = lsqcurvefit(twoGaussFree, p0, theta(:), xSmooth(:), lb, ub, optsLSQfree);
|
||||
pFree = lsqcurvefit(twoGaussFree, p0, theta(:), x(:), lb, ub, optsLSQfree);
|
||||
|
||||
% --- Handle missing secondary peak ---
|
||||
if pFree(4) <= 1e-3
|
||||
% --- Compute RMS of fit vs raw curve in [0, pi/4] ---
|
||||
maskRMS = theta >= 0 & theta <= pi/4;
|
||||
if any(maskRMS)
|
||||
x_raw_subset = x(maskRMS);
|
||||
y_fit_subset = twoGaussFree(pFree, theta(maskRMS));
|
||||
residual_rms_val = sqrt(mean((x_raw_subset - y_fit_subset).^2));
|
||||
end
|
||||
|
||||
%{
|
||||
% --- Handle missing or weak secondary peak ---
|
||||
if pFree(4) <= opts.AmplitudeThreshold
|
||||
pFree(4) = 0; % secondary amplitude = 0
|
||||
pFree(5:6) = NaN; % secondary position & width = NaN
|
||||
end
|
||||
|
||||
|
||||
% --- Discard poor fits ---
|
||||
if rms_val > opts.ResidualThreshold
|
||||
% residual too large → discard secondary
|
||||
pFree(4) = NaN;
|
||||
pFree(5:6) = NaN;
|
||||
end
|
||||
%}
|
||||
|
||||
fitMaxTheta = max(pFree(5), pFree(2));
|
||||
thetaFine = linspace(0, opts.MaxTheta, 500);
|
||||
yFit = twoGaussFree(pFree, thetaFine);
|
||||
@ -219,6 +238,7 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToAngularSpectralDistribu
|
||||
fitResults(k).thetaFine = thetaFine;
|
||||
fitResults(k).yFit = yFit;
|
||||
fitResults(k).isValid = true;
|
||||
fitResults(k).residual = residual_rms_val;
|
||||
fitResults(k).fitMaxTheta = fitMaxTheta;
|
||||
end
|
||||
catch
|
||||
@ -242,11 +262,11 @@ function s = fillNaNStruct()
|
||||
'yFit', nan, ...
|
||||
'thetaFine', nan, ...
|
||||
'isValid', false, ...
|
||||
'residual', nan, ...
|
||||
'fitMaxTheta', nan ...
|
||||
);
|
||||
end
|
||||
|
||||
|
||||
% --- Helper: return struct with zeros (no secondary peak found) ---
|
||||
function s = fillZeroStruct(theta, fitMaxTheta)
|
||||
s = struct( ...
|
||||
@ -256,6 +276,7 @@ function s = fillZeroStruct(theta, fitMaxTheta)
|
||||
'yFit', zeros(size(theta)), ...
|
||||
'thetaFine', zeros(1,500), ...
|
||||
'isValid', false, ...
|
||||
'residual', zeros(size(theta)), ...
|
||||
'fitMaxTheta', fitMaxTheta ...
|
||||
);
|
||||
end
|
||||
|
||||
@ -33,7 +33,7 @@ function [fitResults, rawCurves] = fitTwoGaussianCurvesToRadialSpectralDistribut
|
||||
|
||||
Ncurves = numel(S_k_all);
|
||||
fitResults = struct('pFit',[],'kFit',[],'xFit',[],'yFit',[],...
|
||||
'kFine',[],'isValid',[],'fitMaxKRho',[]);
|
||||
'kFine',[],'isValid',[],'residual',[],'fitMaxKRho',[]);
|
||||
rawCurves = struct('x', cell(1, Ncurves), 'k', cell(1, Ncurves));
|
||||
|
||||
for k = 1:Ncurves
|
||||
@ -130,6 +130,7 @@ function s = fillNaNStructRadial()
|
||||
'yFit', nan, ...
|
||||
'kFine', nan, ...
|
||||
'isValid', false, ...
|
||||
'residual', nan, ...
|
||||
'fitMaxKRho', nan ...
|
||||
);
|
||||
end
|
||||
@ -143,6 +144,7 @@ function s = fillZeroStructRadial(kVals, MaxKRho)
|
||||
'yFit', zeros(size(kVals)), ...
|
||||
'kFine', zeros(1,500), ...
|
||||
'isValid', false, ...
|
||||
'residual', zeros(size(kVals)), ...
|
||||
'fitMaxKRho', MaxKRho ...
|
||||
);
|
||||
end
|
||||
|
||||
@ -112,7 +112,7 @@ function runInteractiveASDTwoGaussianFitGUI(spectral_analysis_results)
|
||||
'Callback',@(~,~) refitAll());
|
||||
|
||||
%% --- Fit info text boxes ---
|
||||
infoLabels = {'Fit validity:', 'Peak positions (rad):', 'Amplitudes:'};
|
||||
infoLabels = {'Fit validity:','Peak positions:','Amplitudes:','Residual (RMS):'};
|
||||
hInfoText = gobjects(numel(infoLabels),1);
|
||||
for i = 1:numel(infoLabels)
|
||||
yText = 0.22 - 0.07*(i-1);
|
||||
@ -195,11 +195,19 @@ function runInteractiveASDTwoGaussianFitGUI(spectral_analysis_results)
|
||||
hInfoText(2).String = '-';
|
||||
hInfoText(3).String = '-';
|
||||
end
|
||||
|
||||
if isfield(fit,'residual') && ~isempty(fit.residual)
|
||||
hInfoText(4).String = sprintf('%.4f', fit.residual);
|
||||
else
|
||||
hInfoText(4).String = '-';
|
||||
end
|
||||
|
||||
else
|
||||
hInfoText(1).String = 'Invalid';
|
||||
hInfoText(1).ForegroundColor = [0.7 0 0];
|
||||
hInfoText(2).String = '-';
|
||||
hInfoText(3).String = '-';
|
||||
hInfoText(4).String = '-';
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@ -113,7 +113,7 @@ function runInteractiveRSDTwoGaussianFitGUI(spectral_analysis_results)
|
||||
'Callback',@(~,~) refitAll());
|
||||
|
||||
%% --- Fit info text boxes ---
|
||||
infoLabels = {'Fit validity:','Peak positions:','Amplitudes:'};
|
||||
infoLabels = {'Fit validity:','Peak positions:','Amplitudes:','Residual (RMS):'};
|
||||
hInfoText = gobjects(numel(infoLabels),1);
|
||||
for i = 1:numel(infoLabels)
|
||||
yText = 0.22 - 0.07*(i-1);
|
||||
@ -129,7 +129,7 @@ function runInteractiveRSDTwoGaussianFitGUI(spectral_analysis_results)
|
||||
'HorizontalAlignment','left',...
|
||||
'BackgroundColor','w');
|
||||
end
|
||||
|
||||
|
||||
%% --- Compute initial fits ---
|
||||
[fitResults, rawData] = Analyzer.fitTwoGaussianCurvesToRadialSpectralDistribution(...
|
||||
spectral_analysis_results.S_k_all,...
|
||||
@ -195,11 +195,18 @@ function runInteractiveRSDTwoGaussianFitGUI(spectral_analysis_results)
|
||||
hInfoText(1).String='Valid'; hInfoText(1).ForegroundColor=[0 0.5 0];
|
||||
hInfoText(2).String='-'; hInfoText(3).String='-';
|
||||
end
|
||||
|
||||
if isfield(fit,'residual') && ~isempty(fit.residual)
|
||||
hInfoText(4).String = sprintf('%.4f', fit.residual);
|
||||
else
|
||||
hInfoText(4).String = '-';
|
||||
end
|
||||
|
||||
else
|
||||
hInfoText(1).String='Invalid'; hInfoText(1).ForegroundColor=[0.7 0 0];
|
||||
hInfoText(2).String='-'; hInfoText(3).String='-';
|
||||
hInfoText(2).String='-'; hInfoText(3).String='-'; hInfoText(4).String='-';
|
||||
end
|
||||
|
||||
|
||||
set(gca, 'YScale', 'log');
|
||||
|
||||
xlabel(hAx,'k_\rho'); ylabel(hAx,'Normalized amplitude');
|
||||
|
||||
124
Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m
Normal file
124
Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m
Normal file
@ -0,0 +1,124 @@
|
||||
function plotTwoModeGaussianFitsOnRawASD(fitResults, rawCurves, nRows, nCols, varargin)
|
||||
%% plotTwoModeGaussianFitsOnRawASD
|
||||
% Author: Karthik
|
||||
% Date: 2025-10-16
|
||||
% Version: 1.0
|
||||
%
|
||||
% Description:
|
||||
% Plots raw radial spectral curves with their two-Gaussian fits in a
|
||||
% compact, paginated layout. Each page shows up to nRows × nCols subplots.
|
||||
%
|
||||
% Inputs:
|
||||
% fitResults - struct array from fitTwoGaussianCurvesToRadialSpectralDistribution
|
||||
% rawCurves - struct array with fields:
|
||||
% .x - raw normalized amplitudes
|
||||
% .k - corresponding k_rho values
|
||||
% nRows - number of subplot rows per page (default: 8)
|
||||
% nCols - number of subplot columns per page (default: 12)
|
||||
% varargin - name-value pairs:
|
||||
% 'FontName', 'FontSize', 'SkipLivePlot', 'SkipSaveFigures',
|
||||
% 'SaveDirectory', 'RepsPerPage'
|
||||
|
||||
% --- Parse optional name-value pairs ---
|
||||
p = inputParser;
|
||||
addParameter(p, 'FontName', 'Arial', @ischar);
|
||||
addParameter(p, 'FontSize', 12, @isnumeric);
|
||||
addParameter(p, 'SkipLivePlot', false, @islogical);
|
||||
addParameter(p, 'SkipSaveFigures', true, @islogical);
|
||||
addParameter(p, 'SaveDirectory', pwd, @ischar);
|
||||
addParameter(p, 'RepsPerPage', [], @isnumeric); % optional override of nRows*nCols
|
||||
parse(p, varargin{:});
|
||||
opts = p.Results;
|
||||
|
||||
% --- Default subplot grid if not provided ---
|
||||
if nargin < 3 || isempty(nRows), nRows = 8; end
|
||||
if nargin < 4 || isempty(nCols), nCols = 12; end
|
||||
|
||||
plotsPerPage = nRows * nCols;
|
||||
if ~isempty(opts.RepsPerPage)
|
||||
plotsPerPage = opts.RepsPerPage;
|
||||
end
|
||||
|
||||
Ncurves = numel(rawCurves);
|
||||
numPages = ceil(Ncurves / plotsPerPage);
|
||||
|
||||
% --- Setup save directory if needed ---
|
||||
if ~opts.SkipSaveFigures
|
||||
saveFolder = fullfile(opts.SaveDirectory, 'Results', 'SavedFigures', 'TwoGaussianFits');
|
||||
if ~exist(saveFolder, 'dir')
|
||||
mkdir(saveFolder);
|
||||
end
|
||||
end
|
||||
|
||||
for pageIdx = 1:numPages
|
||||
repStart = (pageIdx-1)*plotsPerPage + 1;
|
||||
repEnd = min(pageIdx*plotsPerPage, Ncurves);
|
||||
repSubset = repStart:repEnd;
|
||||
N_rows = numel(repSubset);
|
||||
|
||||
% --- Create figure ---
|
||||
if ~opts.SkipLivePlot
|
||||
fig = figure('Color', 'w', 'Position', [50 50 1600 900]);
|
||||
else
|
||||
fig = figure('Color', 'w', 'Position', [50 50 1600 900], 'Visible', 'off');
|
||||
end
|
||||
|
||||
t = tiledlayout(fig, ceil(N_rows/nCols), nCols, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
title(t, sprintf('Two-Gaussian fits | Page %d/%d', pageIdx, numPages), ...
|
||||
'FontSize', opts.FontSize + 2, 'FontWeight', 'bold', 'FontName', opts.FontName);
|
||||
|
||||
for r = 1:N_rows
|
||||
k = repSubset(r);
|
||||
nexttile;
|
||||
hold on; grid on; box on;
|
||||
|
||||
% --- Plot raw curve ---
|
||||
if isfield(rawCurves, 'x') && isfield(rawCurves, 'theta') && ...
|
||||
~isempty(rawCurves(k).x) && all(isfinite(rawCurves(k).x))
|
||||
plot(rawCurves(k).theta, rawCurves(k).x, 'k.-', ...
|
||||
'LineWidth', 1, 'DisplayName', 'Raw data');
|
||||
end
|
||||
|
||||
% --- Overlay fit if valid ---
|
||||
if k <= numel(fitResults)
|
||||
fit = fitResults(k);
|
||||
if isfield(fit, 'isValid') && ~isempty(fit.isValid) && fit.isValid && ...
|
||||
isfield(fit, 'thetaFine') && isfield(fit, 'yFit') && ...
|
||||
~isempty(fit.thetaFine) && all(isfinite(fit.yFit))
|
||||
plot(fit.thetaFine, fit.yFit, 'r-', ...
|
||||
'LineWidth', 1.2, 'DisplayName', 'Two-Gaussian fit');
|
||||
end
|
||||
end
|
||||
|
||||
% --- Labels and title ---
|
||||
xlabel('\theta (rad)', 'FontSize', opts.FontSize);
|
||||
ylabel('N. A.', 'FontSize', opts.FontSize);
|
||||
title(sprintf('Curve %d', k), 'FontSize', opts.FontSize - 2, 'Interpreter', 'none');
|
||||
|
||||
hold off;
|
||||
end
|
||||
|
||||
% --- Render live plot ---
|
||||
if ~opts.SkipLivePlot
|
||||
drawnow;
|
||||
end
|
||||
|
||||
% --- Save figure ---
|
||||
if ~opts.SkipSaveFigures
|
||||
saveFileName = sprintf('TwoGaussianFits_page_%02d.png', pageIdx);
|
||||
if exist('Plotter', 'class')
|
||||
Plotter.saveFigure(fig, ...
|
||||
'SaveFileName', saveFileName, ...
|
||||
'SaveDirectory', saveFolder, ...
|
||||
'SkipSaveFigures', opts.SkipSaveFigures);
|
||||
else
|
||||
saveas(fig, fullfile(saveFolder, saveFileName));
|
||||
end
|
||||
end
|
||||
|
||||
% --- Close invisible figure to free memory ---
|
||||
if opts.SkipLivePlot
|
||||
close(fig);
|
||||
end
|
||||
end
|
||||
end
|
||||
@ -1,5 +1,5 @@
|
||||
function plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, nRows, nCols, varargin)
|
||||
%% plotTwoGaussianFitsWithPagination
|
||||
function plotTwoModeGaussianFitsOnRawRSD(fitResults, rawCurves, nRows, nCols, varargin)
|
||||
%% plotTwoModeGaussianFitsOnRawRSD
|
||||
% Author: Karthik
|
||||
% Date: 2025-10-16
|
||||
% Version: 1.0
|
||||
@ -167,10 +167,31 @@ Plotter.plotSpectralDistributionCumulants(results, ...
|
||||
'PositionThreshold', pi/15, ...
|
||||
'AmplitudeThreshold', 0.15);
|
||||
|
||||
% --- Post-fit diagnostics ---
|
||||
pFit_all = vertcat(fitResults.pFit); % Ncurves x Nparams
|
||||
isValid_all = [fitResults.isValid]; % logical vector
|
||||
|
||||
% Extract last three parameters (usually secondary peak amplitude, position, width)
|
||||
pTail = pFit_all(:, end-2:end);
|
||||
|
||||
% Count curves where *any* of these last three values are NaN or zero
|
||||
numWithNaN = sum(any(isnan(pTail), 2));
|
||||
numWithZero = sum(any(pTail == 0, 2));
|
||||
|
||||
% Count of valid fits
|
||||
numValid = sum(isValid_all);
|
||||
|
||||
% Display summary
|
||||
fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
|
||||
fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
|
||||
fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
|
||||
fprintf('Curves marked as isValid = true: %d\n', numValid);
|
||||
fprintf('----------------------------------------\n\n');
|
||||
|
||||
%% ------------------ Plot Fits on Raw ------------------
|
||||
options.skipLivePlot = true;
|
||||
options.skipSaveFigures = false;
|
||||
Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
|
||||
Plotter.plotTwoModeGaussianFitsOnRawASD(fitResults, rawCurves, 4, 6, ...
|
||||
'SkipLivePlot', options.skipLivePlot, ...
|
||||
'SkipSaveFigures', options.skipSaveFigures, ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
|
||||
@ -167,10 +167,31 @@ Plotter.plotSpectralDistributionCumulants(results, ...
|
||||
'PositionThreshold', pi/15, ...
|
||||
'AmplitudeThreshold', 0.15);
|
||||
|
||||
% --- Post-fit diagnostics ---
|
||||
pFit_all = vertcat(fitResults.pFit); % Ncurves x Nparams
|
||||
isValid_all = [fitResults.isValid]; % logical vector
|
||||
|
||||
% Extract last three parameters (usually secondary peak amplitude, position, width)
|
||||
pTail = pFit_all(:, end-2:end);
|
||||
|
||||
% Count curves where *any* of these last three values are NaN or zero
|
||||
numWithNaN = sum(any(isnan(pTail), 2));
|
||||
numWithZero = sum(any(pTail == 0, 2));
|
||||
|
||||
% Count of valid fits
|
||||
numValid = sum(isValid_all);
|
||||
|
||||
% Display summary
|
||||
fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
|
||||
fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
|
||||
fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
|
||||
fprintf('Curves marked as isValid = true: %d\n', numValid);
|
||||
fprintf('----------------------------------------\n\n');
|
||||
|
||||
%% ------------------ Plot Fits on Raw ------------------
|
||||
options.skipLivePlot = true;
|
||||
options.skipSaveFigures = false;
|
||||
Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
|
||||
Plotter.plotTwoModeGaussianFitsOnRawASD(fitResults, rawCurves, 4, 6, ...
|
||||
'SkipLivePlot', options.skipLivePlot, ...
|
||||
'SkipSaveFigures', options.skipSaveFigures, ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
@ -318,6 +339,100 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
|
||||
'SaveFileName', 'SecondaryPeakWidthMeanWithSEM.fig', ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
|
||||
%% ------------------ 9. Compare mean of individual fits vs fit of mean ------------------
|
||||
|
||||
paramMap = struct('A1',1,'mu1',2,'sigma1',3,'A2',4,'mu2',5,'sigma2',6);
|
||||
paramNames = {'A2', 'mu2', 'sigma2'};
|
||||
paramLabels = {'Secondary peak amplitude', 'Secondary peak position (\theta, rad)', 'Secondary peak width (\sigma, rad)'};
|
||||
|
||||
meanOfIndividualFits = struct();
|
||||
fitOfMeanCurve = struct();
|
||||
|
||||
for p = 1:numel(paramNames)
|
||||
paramName = paramNames{p};
|
||||
paramLabel = paramLabels{p};
|
||||
paramIdx = paramMap.(paramName);
|
||||
|
||||
% --- Mean of individual fits per scan_reference_values ---
|
||||
meanVals = nan(size(scan_reference_values)); % already unique
|
||||
for k = 1:numel(scan_reference_values)
|
||||
thisAlpha = scan_reference_values(k);
|
||||
mask = scan_parameter_values == thisAlpha; % pick all repetitions
|
||||
|
||||
pvals = arrayfun(@(f) f.pFit(paramIdx), fitResults(mask), 'UniformOutput', true);
|
||||
meanVals(k) = mean(pvals(~isnan(pvals) & isfinite(pvals)), 'omitnan');
|
||||
end
|
||||
meanOfIndividualFits.(paramName) = meanVals;
|
||||
|
||||
% --- Fit of mean curve ---
|
||||
switch paramName
|
||||
case 'A2'
|
||||
fitOfMeanCurve.(paramName) = amp2_vals;
|
||||
case 'mu2'
|
||||
fitOfMeanCurve.(paramName) = mu2_vals;
|
||||
case 'sigma2'
|
||||
fitOfMeanCurve.(paramName) = sigma2_vals;
|
||||
end
|
||||
|
||||
% --- Plot comparison ---
|
||||
fig = figure('Name', 'Comparison', 'NumberTitle', 'off');
|
||||
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
|
||||
|
||||
hold on;
|
||||
plot(scan_reference_values, meanOfIndividualFits.(paramName), 'o-', 'LineWidth', 1.8, ...
|
||||
'DisplayName', 'Mean of Individual Fits');
|
||||
plot(scan_reference_values, fitOfMeanCurve.(paramName), 's--', 'LineWidth', 1.8, ...
|
||||
'DisplayName', 'Fit of Mean Curve');
|
||||
hold off;
|
||||
ylim([0, 1.1])
|
||||
xlabel('\alpha (degrees)', 'FontName', options.font, 'FontSize', 16);
|
||||
ylabel(paramLabel, 'FontName', options.font, 'FontSize', 16);
|
||||
title(options.titleString, 'FontName', options.font, 'FontSize', 16);
|
||||
legend('Location', 'northeast');
|
||||
grid on;
|
||||
set(gca, 'FontSize', 14);
|
||||
|
||||
if ~options.skipSaveFigures
|
||||
saveas(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.fig']));
|
||||
exportgraphics(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.png']), 'Resolution', 300);
|
||||
end
|
||||
end
|
||||
|
||||
%% ------------------ 5. Save Fit Results ------------------
|
||||
outputDir = fullfile(options.saveDirectory, 'DToS_FitData');
|
||||
if ~exist(outputDir, 'dir')
|
||||
mkdir(outputDir);
|
||||
end
|
||||
|
||||
% --- Define filenames for each variable ---
|
||||
file_scanParamVals = fullfile(outputDir, 'scan_parameter_values');
|
||||
file_scanRefVals = fullfile(outputDir, 'scan_reference_values');
|
||||
file_fitResults = fullfile(outputDir, 'fitResults.mat');
|
||||
file_rawCurves = fullfile(outputDir, 'rawCurves.mat');
|
||||
file_spectralResults = fullfile(outputDir, 'spectral_analysis_results.mat');
|
||||
file_meanFits = fullfile(outputDir, 'fitOfMeanCurve.mat');
|
||||
file_indivFits = fullfile(outputDir, 'meanOfIndividualFits.mat');
|
||||
|
||||
% --- Save each variable separately ---
|
||||
save(file_scanParamVals, 'scan_parameter_values', '-v7.3');
|
||||
save(file_scanRefVals, 'scan_reference_values', '-v7.3');
|
||||
save(file_fitResults, 'fitResults', '-v7.3');
|
||||
save(file_rawCurves, 'rawCurves', '-v7.3');
|
||||
save(file_spectralResults, 'spectral_analysis_results', '-v7.3');
|
||||
save(file_meanFits, 'fitOfMeanCurve', '-v7.3');
|
||||
save(file_indivFits, 'meanOfIndividualFits', '-v7.3');
|
||||
|
||||
fprintf('\n--- Fit data saved successfully ---\n');
|
||||
fprintf('scanParamVals: %s\n', file_scanParamVals);
|
||||
fprintf('scanRefVals: %s\n', file_scanRefVals);
|
||||
fprintf('fitResults: %s\n', file_fitResults);
|
||||
fprintf('rawCurves: %s\n', file_rawCurves);
|
||||
fprintf('spectral_analysis_results:%s\n', file_spectralResults);
|
||||
fprintf('fitOfMeanCurve: %s\n', file_meanFits);
|
||||
fprintf('meanOfIndividualFits: %s\n', file_indivFits);
|
||||
fprintf('----------------------------------\n\n');
|
||||
|
||||
|
||||
%% Inspect individual realizations of a single parameter
|
||||
|
||||
% --- Recenter curves first ---
|
||||
|
||||
@ -167,10 +167,31 @@ Plotter.plotSpectralDistributionCumulants(results, ...
|
||||
'PositionThreshold', pi/15, ...
|
||||
'AmplitudeThreshold', 0.15);
|
||||
|
||||
% --- Post-fit diagnostics ---
|
||||
pFit_all = vertcat(fitResults.pFit); % Ncurves x Nparams
|
||||
isValid_all = [fitResults.isValid]; % logical vector
|
||||
|
||||
% Extract last three parameters (usually secondary peak amplitude, position, width)
|
||||
pTail = pFit_all(:, end-2:end);
|
||||
|
||||
% Count curves where *any* of these last three values are NaN or zero
|
||||
numWithNaN = sum(any(isnan(pTail), 2));
|
||||
numWithZero = sum(any(pTail == 0, 2));
|
||||
|
||||
% Count of valid fits
|
||||
numValid = sum(isValid_all);
|
||||
|
||||
% Display summary
|
||||
fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
|
||||
fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
|
||||
fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
|
||||
fprintf('Curves marked as isValid = true: %d\n', numValid);
|
||||
fprintf('----------------------------------------\n\n');
|
||||
|
||||
%% ------------------ Plot Fits on Raw ------------------
|
||||
options.skipLivePlot = true;
|
||||
options.skipSaveFigures = false;
|
||||
Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
|
||||
Plotter.plotTwoModeGaussianFitsOnRawASD(fitResults, rawCurves, 4, 6, ...
|
||||
'SkipLivePlot', options.skipLivePlot, ...
|
||||
'SkipSaveFigures', options.skipSaveFigures, ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
|
||||
@ -167,10 +167,31 @@ Plotter.plotSpectralDistributionCumulants(results, ...
|
||||
'PositionThreshold', pi/15, ...
|
||||
'AmplitudeThreshold', 0.15);
|
||||
|
||||
% --- Post-fit diagnostics ---
|
||||
pFit_all = vertcat(fitResults.pFit); % Ncurves x Nparams
|
||||
isValid_all = [fitResults.isValid]; % logical vector
|
||||
|
||||
% Extract last three parameters (usually secondary peak amplitude, position, width)
|
||||
pTail = pFit_all(:, end-2:end);
|
||||
|
||||
% Count curves where *any* of these last three values are NaN or zero
|
||||
numWithNaN = sum(any(isnan(pTail), 2));
|
||||
numWithZero = sum(any(pTail == 0, 2));
|
||||
|
||||
% Count of valid fits
|
||||
numValid = sum(isValid_all);
|
||||
|
||||
% Display summary
|
||||
fprintf('\n--- Fit diagnostics (last 3 params) ---\n');
|
||||
fprintf('Curves with ≥1 NaN in last three pFit elements: %d\n', numWithNaN);
|
||||
fprintf('Curves with ≥1 zero in last three pFit elements: %d\n', numWithZero);
|
||||
fprintf('Curves marked as isValid = true: %d\n', numValid);
|
||||
fprintf('----------------------------------------\n\n');
|
||||
|
||||
%% ------------------ Plot Fits on Raw ------------------
|
||||
options.skipLivePlot = true;
|
||||
options.skipSaveFigures = false;
|
||||
Plotter.plotTwoModeGaussianFitsOnRaw(fitResults, rawCurves, 4, 6, ...
|
||||
Plotter.plotTwoModeGaussianFitsOnRawASD(fitResults, rawCurves, 4, 6, ...
|
||||
'SkipLivePlot', options.skipLivePlot, ...
|
||||
'SkipSaveFigures', options.skipSaveFigures, ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
@ -318,6 +339,99 @@ Plotter.plotMeanWithSE(scan_reference_values, sigma2_vals, ...
|
||||
'SaveFileName', 'SecondaryPeakWidthMeanWithSEM.fig', ...
|
||||
'SaveDirectory', options.saveDirectory);
|
||||
|
||||
%% ------------------ 9. Compare mean of individual fits vs fit of mean ------------------
|
||||
|
||||
paramMap = struct('A1',1,'mu1',2,'sigma1',3,'A2',4,'mu2',5,'sigma2',6);
|
||||
paramNames = {'A2', 'mu2', 'sigma2'};
|
||||
paramLabels = {'Secondary peak amplitude', 'Secondary peak position (\theta, rad)', 'Secondary peak width (\sigma, rad)'};
|
||||
|
||||
meanOfIndividualFits = struct();
|
||||
fitOfMeanCurve = struct();
|
||||
|
||||
for p = 1:numel(paramNames)
|
||||
paramName = paramNames{p};
|
||||
paramLabel = paramLabels{p};
|
||||
paramIdx = paramMap.(paramName);
|
||||
|
||||
% --- Mean of individual fits per scan_reference_values ---
|
||||
meanVals = nan(size(scan_reference_values)); % already unique
|
||||
for k = 1:numel(scan_reference_values)
|
||||
thisAlpha = scan_reference_values(k);
|
||||
mask = scan_parameter_values == thisAlpha; % pick all repetitions
|
||||
|
||||
pvals = arrayfun(@(f) f.pFit(paramIdx), fitResults(mask), 'UniformOutput', true);
|
||||
meanVals(k) = mean(pvals(~isnan(pvals) & isfinite(pvals)), 'omitnan');
|
||||
end
|
||||
meanOfIndividualFits.(paramName) = meanVals;
|
||||
|
||||
% --- Fit of mean curve ---
|
||||
switch paramName
|
||||
case 'A2'
|
||||
fitOfMeanCurve.(paramName) = amp2_vals;
|
||||
case 'mu2'
|
||||
fitOfMeanCurve.(paramName) = mu2_vals;
|
||||
case 'sigma2'
|
||||
fitOfMeanCurve.(paramName) = sigma2_vals;
|
||||
end
|
||||
|
||||
% --- Plot comparison ---
|
||||
fig = figure('Name', 'Comparison', 'NumberTitle', 'off');
|
||||
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);
|
||||
|
||||
hold on;
|
||||
plot(scan_reference_values, meanOfIndividualFits.(paramName), 'o-', 'LineWidth', 1.8, ...
|
||||
'DisplayName', 'Mean of Individual Fits');
|
||||
plot(scan_reference_values, fitOfMeanCurve.(paramName), 's--', 'LineWidth', 1.8, ...
|
||||
'DisplayName', 'Fit of Mean Curve');
|
||||
hold off;
|
||||
ylim([0, 1.2])
|
||||
xlabel('\alpha (degrees)', 'FontName', options.font, 'FontSize', 16);
|
||||
ylabel(paramLabel, 'FontName', options.font, 'FontSize', 16);
|
||||
title(options.titleString, 'FontName', options.font, 'FontSize', 16);
|
||||
legend('Location', 'northeast');
|
||||
grid on;
|
||||
set(gca, 'FontSize', 14);
|
||||
|
||||
if ~options.skipSaveFigures
|
||||
saveas(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.fig']));
|
||||
exportgraphics(fig, fullfile(options.saveDirectory, ['Compare_' paramName '.png']), 'Resolution', 300);
|
||||
end
|
||||
end
|
||||
|
||||
%% ------------------ 5. Save Fit Results ------------------
|
||||
outputDir = fullfile(options.saveDirectory, 'SToD_FitData');
|
||||
if ~exist(outputDir, 'dir')
|
||||
mkdir(outputDir);
|
||||
end
|
||||
|
||||
% --- Define filenames for each variable ---
|
||||
file_scanParamVals = fullfile(outputDir, 'scan_parameter_values');
|
||||
file_scanRefVals = fullfile(outputDir, 'scan_reference_values');
|
||||
file_fitResults = fullfile(outputDir, 'fitResults.mat');
|
||||
file_rawCurves = fullfile(outputDir, 'rawCurves.mat');
|
||||
file_spectralResults = fullfile(outputDir, 'spectral_analysis_results.mat');
|
||||
file_meanFits = fullfile(outputDir, 'fitOfMeanCurve.mat');
|
||||
file_indivFits = fullfile(outputDir, 'meanOfIndividualFits.mat');
|
||||
|
||||
% --- Save each variable separately ---
|
||||
save(file_scanParamVals, 'scan_parameter_values', '-v7.3');
|
||||
save(file_scanRefVals, 'scan_reference_values', '-v7.3');
|
||||
save(file_fitResults, 'fitResults', '-v7.3');
|
||||
save(file_rawCurves, 'rawCurves', '-v7.3');
|
||||
save(file_spectralResults, 'spectral_analysis_results', '-v7.3');
|
||||
save(file_meanFits, 'fitOfMeanCurve', '-v7.3');
|
||||
save(file_indivFits, 'meanOfIndividualFits', '-v7.3');
|
||||
|
||||
fprintf('\n--- Fit data saved successfully ---\n');
|
||||
fprintf('scanParamVals: %s\n', file_scanParamVals);
|
||||
fprintf('scanRefVals: %s\n', file_scanRefVals);
|
||||
fprintf('fitResults: %s\n', file_fitResults);
|
||||
fprintf('rawCurves: %s\n', file_rawCurves);
|
||||
fprintf('spectral_analysis_results:%s\n', file_spectralResults);
|
||||
fprintf('fitOfMeanCurve: %s\n', file_meanFits);
|
||||
fprintf('meanOfIndividualFits: %s\n', file_indivFits);
|
||||
fprintf('----------------------------------\n\n');
|
||||
|
||||
%% Inspect individual realizations of a single parameter
|
||||
|
||||
% --- Recenter curves first ---
|
||||
|
||||
Loading…
Reference in New Issue
Block a user