From dd4145ccd27a522ce9255f9db61017cd91956a7a Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 17 Oct 2025 18:49:54 +0200 Subject: [PATCH] New and modified fitting and plotting routines (Comparison scripts will be added later). --- ...ssianCurvesToAngularSpectralDistribution.m | 65 +++++---- ...ussianCurvesToRadialSpectralDistribution.m | 4 +- .../runInteractiveASDTwoGaussianFitGUI.m | 10 +- .../runInteractiveRSDTwoGaussianFitGUI.m | 15 ++- .../plotTwoModeGaussianFitsOnRawASD.m | 124 ++++++++++++++++++ ...aw.m => plotTwoModeGaussianFitsOnRawRSD.m} | 4 +- .../runAngularSpectralDistributionAnalysis.m | 23 +++- .../runAngularSpectralDistributionAnalysis.m | 117 ++++++++++++++++- .../runAngularSpectralDistributionAnalysis.m | 23 +++- .../runAngularSpectralDistributionAnalysis.m | 116 +++++++++++++++- 10 files changed, 467 insertions(+), 34 deletions(-) create mode 100644 Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m rename Data-Analyzer/+Plotter/{plotTwoModeGaussianFitsOnRaw.m => plotTwoModeGaussianFitsOnRawRSD.m} (97%) diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m index 890b4ec..baead39 100644 --- a/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToAngularSpectralDistribution.m @@ -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 diff --git a/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m index 875798d..685790a 100644 --- a/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m +++ b/Data-Analyzer/+Analyzer/fitTwoGaussianCurvesToRadialSpectralDistribution.m @@ -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 diff --git a/Data-Analyzer/+Analyzer/runInteractiveASDTwoGaussianFitGUI.m b/Data-Analyzer/+Analyzer/runInteractiveASDTwoGaussianFitGUI.m index 3b079d9..c59b509 100644 --- a/Data-Analyzer/+Analyzer/runInteractiveASDTwoGaussianFitGUI.m +++ b/Data-Analyzer/+Analyzer/runInteractiveASDTwoGaussianFitGUI.m @@ -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 diff --git a/Data-Analyzer/+Analyzer/runInteractiveRSDTwoGaussianFitGUI.m b/Data-Analyzer/+Analyzer/runInteractiveRSDTwoGaussianFitGUI.m index 5b9acb9..4363ec2 100644 --- a/Data-Analyzer/+Analyzer/runInteractiveRSDTwoGaussianFitGUI.m +++ b/Data-Analyzer/+Analyzer/runInteractiveRSDTwoGaussianFitGUI.m @@ -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'); diff --git a/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m b/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m new file mode 100644 index 0000000..8075e87 --- /dev/null +++ b/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawASD.m @@ -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 diff --git a/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRaw.m b/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawRSD.m similarity index 97% rename from Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRaw.m rename to Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawRSD.m index 0f05419..eb357f5 100644 --- a/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRaw.m +++ b/Data-Analyzer/+Plotter/plotTwoModeGaussianFitsOnRawRSD.m @@ -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 diff --git a/Data-Analyzer/+Scripts/BECToDroplets/runAngularSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDroplets/runAngularSpectralDistributionAnalysis.m index 02e7f74..d81a172 100644 --- a/Data-Analyzer/+Scripts/BECToDroplets/runAngularSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToDroplets/runAngularSpectralDistributionAnalysis.m @@ -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); diff --git a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m index 0a7f4d9..f2b3e52 100644 --- a/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToDropletsToStripes/runAngularSpectralDistributionAnalysis.m @@ -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 --- diff --git a/Data-Analyzer/+Scripts/BECToStripes/runAngularSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToStripes/runAngularSpectralDistributionAnalysis.m index 6631091..c86fc11 100644 --- a/Data-Analyzer/+Scripts/BECToStripes/runAngularSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToStripes/runAngularSpectralDistributionAnalysis.m @@ -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); diff --git a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m index 85c59dc..a97531c 100644 --- a/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m +++ b/Data-Analyzer/+Scripts/BECToStripesToDroplets/runAngularSpectralDistributionAnalysis.m @@ -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 ---