From a7e78cd6f510ee4507531a65b8cbcc01c1b4b991 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 8 Jan 2025 16:09:21 +0000 Subject: [PATCH 01/33] Fix bug in estimate_scale_factor - background already added to ycalc Didn't break a test as bg is all zeros in unit test - have now added a constant bg so this bug would have been caught. --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 3 ++- swfiles/sw_fitpowder.m | 6 +++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 925ce96d..1bc2ff72 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -406,11 +406,12 @@ function test_calc_uncertainty(testCase) function test_estimate_scale_factor(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... testCase.fit_func, testCase.j1); + out.set_bg_parameters(3, 0.05); % set constant bg out.powspec_args.dE = 0.1; % constant energy resolution out.powspec_args.hermit = true; out.estimate_scale_factor() expected_fitpow = testCase.default_fitpow; - expected_fitpow.params(end) = 17.6; + expected_fitpow.params(end) = 17.36; testCase.verify_results(out, expected_fitpow, ... testCase.default_fields, 'abs_tol', 1e-1); end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index e98efd89..1da91f3c 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -573,7 +573,11 @@ function estimate_scale_factor(obj) params = obj.params; params(end) = 1; [ycalc, bg] = calc_spinwave_spec(obj, params); - scale = max(obj.y - mean(bg(:)), [], "all")/max(ycalc, [], "all"); + if obj.ndim == 1 && obj.background_strategy=="planar" + % integrate nQ |Q| points for each cut + bg = obj.rebin_powspec_to_1D_cuts(bg); + end + scale = max(obj.y - bg, [], "all")/max(ycalc - bg, [], "all"); obj.params(end) = scale; end From f7dd48e8e765989e874fe10a3a375fc235983b30 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 8 Jan 2025 16:36:26 +0000 Subject: [PATCH 02/33] Add new function to fit bg and scale parameters to all data --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 15 +++++++++++++++ swfiles/sw_fitpowder.m | 13 +++++++++++++ 2 files changed, 28 insertions(+) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 1bc2ff72..dffdae25 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -342,6 +342,21 @@ function test_fit_background(testCase, fit_params) 'abs_tol', 1e-4); end + function test_fit_background_and_scale(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.fix_bg_parameters(1:2); % fix slopes of background to 0 + out.set_bg_parameters(3, 1.5); % initial guess + out.fit_background_and_scale(); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(end-1) = 0.0029; + expected_fitpow.params(end) = 15.47; + expected_fitpow.bounds(2:3,:) = 0; % fixed bg slopes + testCase.verify_results(out, expected_fitpow, ... + testCase.default_fields, ... + 'abs_tol', 1e-3); + end + function test_calc_cost_func_of_background_indep(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent", 2); diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 1da91f3c..b819c9b4 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -504,6 +504,19 @@ function clear_cache(obj) varargin{:}); end + function varargout = fit_background_and_scale(obj, varargin) + % fix all model parameters + initial_bounds = obj.bounds; % store so bounds can be reset + obj.fix_model_parameters(1:obj.nparams_model); + varargout = cell(1,nargout(obj.optimizer)); + [varargout{:}] = obj.fit(varargin{:}); + % reset bounds + obj.bounds = initial_bounds; + % overwrite non model parameters + obj.params = varargout{1}(:); % assume first output (as for ndbase) + end + + function [param_errors, varargout] = calc_uncertainty(obj, params, varargin) % Function to estimate parameter errors by evaluating hessian. % By default will only evaluate derivatives for parameters From 127156ea393a273f6f61b0086042b1d997d8e084 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 8 Jan 2025 16:51:51 +0000 Subject: [PATCH 03/33] Support any ndbase optimizer in fit_background --- swfiles/sw_fitpowder.m | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index b819c9b4..b3562ad7 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -549,7 +549,7 @@ function clear_cache(obj) varargout = {cov}; end - function fit_background(obj, varargin) + function varargout = fit_background(obj, varargin) if isempty(obj.ibg) obj.find_indices_and_mean_of_bg_bins(); end @@ -567,18 +567,13 @@ function fit_background(obj, varargin) else fobj = @obj.calc_cost_func_of_background; % minimise scalar end - - [bg_params, ~, stat] = ndbase.simplex([], fobj, ... - obj.params(ibg_par)', ... - 'lb', obj.bounds(ibg_par,1)', ... - 'ub', obj.bounds(ibg_par,2)', ... - varargin{:}); - if stat.exitFlag == 1 - obj.params(ibg_par) = bg_params; - else - warning('spinw:sw_fitpowder:fit_background', ... - 'Fit failed - parameters not updated.'); - end + varargout = cell(1,nargout(obj.optimizer)); + [varargout{:}] = obj.optimizer([], fobj, ... + obj.params(ibg_par)', ... + 'lb', obj.bounds(ibg_par,1)', ... + 'ub', obj.bounds(ibg_par,2)', ... + varargin{:}); + obj.params(ibg_par) = varargout{1}(:); end function estimate_scale_factor(obj) From d04fe5ccd3804a3081ea64d14c539a36a346e1dc Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Jan 2025 09:38:27 +0000 Subject: [PATCH 04/33] Add method to plot background region on figures --- swfiles/sw_fitpowder.m | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index b3562ad7..906f1a65 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -649,7 +649,7 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) ax = subplot(1,1,1); box on; hold on; h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, y); - h.AlphaData = obj.y > 0; % make empty bins transparent + h.AlphaData = double(obj.y > 0); % make empty bins transparent cbar = colorbar(ax); cbar.Label.String = "Intensity"; xlabel(ax, "$\left|Q\right| (\AA^{-1})$", 'interpreter','latex'); @@ -658,6 +658,16 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) xlim(ax, [obj.modQ_cens(1), obj.modQ_cens(end)]); end + function plot_background_region(obj) + if ~isempty(obj.ibg) + if obj.ndim == 1 + obj.plot_background_region_1d() + else + obj.plot_background_region_2d() + end + end + end + function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) % optionally plot ycalc provided, otherwise will plot % fitpow.ycalc if not empty @@ -809,6 +819,29 @@ function fix_parameters(obj, iparams) obj.ibg = iseed(isort(ipt:end)); end end + + function plot_background_region_2d(obj) + im = findobj(gca, 'type', 'Image'); + is_valid = ~isempty(im) && all(size(im.CData) == size(obj.y)); + assert(is_valid, 'sw_fitpowder:invalidinput', ... + ['This function requires active axes to have an' ... + ' image/colorfill plot of the data.']); + im.AlphaData(obj.ibg) = 0.5; + end + + function plot_background_region_1d(obj) + fig = gcf(); + axs = flip(fig.Children); + assert(obj.ncuts==numel(axs), 'sw_fitpowder:invalidinput', ... + ['This function requires active figure to have same' ... + ' number of axes as 1D cuts']); + [ien, icuts] = ind2sub(size(obj.y), obj.ibg); + for icut = 1:obj.ncuts + ibg_cut = ien(icuts == icut); + plot(axs(icut), obj.ebin_cens(ibg_cut), ... + obj.y(ibg_cut,icut), 'xb') + end + end end methods (Static=true, Hidden=true, Access = private) function data_struct = convert_horace_to_struct(data) From 9aa01e429ff8652c0cc533e5084c8c7262c11dca Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Jan 2025 12:10:06 +0000 Subject: [PATCH 05/33] Add function (and tests) to set background region --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 22 +++++++++++++++ swfiles/sw_fitpowder.m | 27 ++++++++++++++++++- 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index dffdae25..0e611c4b 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -466,6 +466,28 @@ function test_add_1Dcuts_after_2D_data(testCase) expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; % integrtate over nQ pts testCase.verify_results(out, expected_fitpow); end + + function test_set_bg_region_data_2d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.set_bg_region(0,1.5); % for all Q + out.set_bg_region(2.5,3.5,4.5,inf); % for highest Q + expected_fitpow = testCase.default_fitpow; + expected_fitpow.ibg = [1;4;6]; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_region_data_1d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + out.set_bg_region(0,1.5); % for all cuts + out.set_bg_region(2.5,3.5,2); % for last cut + expected_fitpow.ibg = [1;4;6]; + testCase.verify_results(out, expected_fitpow); + end + end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 906f1a65..2866f26f 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -415,6 +415,31 @@ function clear_cache(obj) obj.reset_errors_of_bg_bins() end + function set_bg_region(obj, en_lo, en_hi, varargin) + % 2D data: obj.set_bg_region(en_lo, en_hi, q_lo, q_hi) + % 1D data: obj.set_bg_region(en_lo, en_hi, icuts) + % both: + % obj.set_bg_region(en_lo, en_hi) % apply to all cuts/Q + if isempty(varargin) + iq = 1:size(obj.y, 2); + elseif obj.ndim == 1 && numel(varargin)==1 + iq = varargin{1}; + elseif obj.ndim==2 && numel(varargin)==2 + [q_lo, q_hi] = varargin{:}; + iq = obj.modQ_cens < q_hi & obj.modQ_cens > q_lo; + else + error('spinw:sw_fitpowder:set_bg_region', ... + ['Wrong number of additional aruguments: 2 required' ... + ' for 2D problem (qlo and qhi)' ... + 'and 1 required for 1D problem (icut indices)']); + end + ien = obj.ebin_cens < en_hi & obj.ebin_cens > en_lo; + mask = false(size(obj.y)); + mask(ien, iq) = true; + mask = mask & isfinite(obj.y); + obj.ibg = [obj.ibg; find(mask(:))]; + end + function [ycalc, bg] = calc_spinwave_spec(obj, params) model_params = reshape(params(1:obj.nparams_model), 1, []); if obj.do_cache && ~isempty(obj.model_params_cached) && all(abs(model_params - obj.model_params_cached) < 1e-10) @@ -819,7 +844,7 @@ function fix_parameters(obj, iparams) obj.ibg = iseed(isort(ipt:end)); end end - + function plot_background_region_2d(obj) im = findobj(gca, 'type', 'Image'); is_valid = ~isempty(im) && all(size(im.CData) == size(obj.y)); From 4c25403e739e8a3d5c78ea5cf45d6c503e406f18 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Jan 2025 16:05:29 +0000 Subject: [PATCH 06/33] Set sensible background when changing strategy --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 32 +++++++++++- swfiles/sw_fitpowder.m | 52 ++++++++++++++----- 2 files changed, 68 insertions(+), 16 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 0e611c4b..0ffceed5 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -76,13 +76,41 @@ function test_init_data_1d_indep_bg(testCase) testCase.verify_results(out, expected_fitpow); end - function test_set_background_strategy(testCase) + function test_set_background_strategy_to_planar(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent"); + % set some background parameters (correpsonding to planar bg + % with slope_en=slope_q=intercept =2 + out.set_bg_parameters(1, 2); % en_slope = 2 + out.set_bg_parameters(2, 10, 1); % intercept = 10 for cut 1 + out.set_bg_parameters(2, 12, 2); % intercept = 12 for cut 2 out.set_background_strategy("planar"); expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(2:end-1) = 2; expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; - testCase.verify_results(out, expected_fitpow); + testCase.verify_results(out, expected_fitpow, ... + testCase.default_fields, ... + 'abs_tol', 1e-10); + end + + function test_set_background_strategy_to_indep(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "planar"); + % set some background parameters (correpsonding to planar bg + % with slope_en=slope_q=intercept =2 + out.set_bg_parameters(1:3, [2,2,2]); % en_slope = 2 + out.set_background_strategy("independent"); + expected_fitpow = testCase.default_fitpow; + % add extra background param + expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + expected_fitpow.params(2:2:end-1) = 2; % en slope + expected_fitpow.params(3) = 10; % intercept = 10 for cut 1 + expected_fitpow.params(5) = 12; % intercept = 12 for cut 2 + expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + testCase.verify_results(out, expected_fitpow, ... + testCase.default_fields, ... + 'abs_tol', 1e-10); end function test_replace_2D_data_with_1D_cuts(testCase) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 2866f26f..e65935df 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -215,22 +215,44 @@ function initialise_background_parameters_and_bounds(obj) end function set_background_strategy(obj, strategy) - if (strategy == "independent" && obj.ndim==1) || strategy == "planar" - obj.background_strategy = strategy; - else - error('sw_fitpowder:set_background_strategy', ... - 'Parameter indices supplied must be within number of model parameters'); - end - if obj.background_strategy == "planar" + if strategy == "planar" obj.fbg = obj.fbg_planar; - else + elseif strategy == "independent" + if obj.ndim == 2 + error('sw_fitpowder:invalidinput', ... + ['Can only have independent background strategy' ... + ' for 1D datasets.']) + end obj.fbg = obj.fbg_indep; end + % store previous bg parameters before reset size of param array + ibg_par = (obj.nparams_model+1):(numel(obj.params)-1); + bg_pars = obj.params(ibg_par); + % reset param array d bounds etc. + obj.background_strategy = strategy; obj.nparams_bg = obj.get_nparams_in_background_func(); if ~isempty(obj.params) - warning('sw_fitpowder:set_background_strategy', ... - 'Overwriting background parmaweters with 0.'); obj.initialise_background_parameters_and_bounds(); + if any(abs(bg_pars) > obj.zero_abs_tol) && obj.ndim == 1 + % set good guess for background pars + modQs = obj.get_modQ_cens_of_cuts(); + if obj.background_strategy == "independent" + % keep same energy slope (scaled by nQ) + obj.set_bg_parameters(1, bg_pars(1)); + % set intercept (zero energy) for reach cut separately + bg_pars = num2cell(bg_pars); + intercepts = obj.fbg_planar(0, modQs, bg_pars{:}); + for icut = 1:obj.ncuts + obj.set_bg_parameters(2, intercepts(icut), icut); + end + elseif obj.background_strategy == "planar" + % average the energy slope + slope_en = mean(bg_pars(1:2:end)); + % fit intercepts to linear Q + pval = polyfit(modQs, bg_pars(2:2:end), 1); + obj.set_bg_parameters(1:3, [slope_en, pval]); + end + end end end @@ -356,9 +378,6 @@ function replace_2D_data_with_1D_cuts(obj, qmins, qmaxs, background_strategy) obj.add_data(cuts); if nargin > 3 && background_strategy ~= obj.background_strategy obj.set_background_strategy(background_strategy) - warning('spinw:sw_fitpowder:replace_2D_data_with_1D_cuts', ... - ['Background strategy changed - background ' ... - 'parameters and bounds will be cleared.']); end end @@ -648,7 +667,7 @@ function plot_result(obj, params, varargin) end function plot_1d_cuts_on_data(obj, ycalc, varargin) - modQs = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); + modQs = obj.get_modQ_cens_of_cuts(); for icut = 1:obj.ncuts ax = subplot(1, obj.ncuts, icut); hold on; box on; @@ -739,6 +758,11 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) resid = resid(isfinite(resid) & e > obj.zero_abs_tol); end + + function modQ_cens = get_modQ_cens_of_cuts(obj) + modQ_cens = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); + end + function cut = cut_2d_data(obj, qmin, qmax) assert(obj.ndim ==2, ... 'sw_fitpowder:invalidinput', ... From 2869cff75d201485b44129d505ee88f0aa5d0c41 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 9 Jan 2025 16:14:47 +0000 Subject: [PATCH 07/33] Handle edges/nans properly in avg. of error when taking 1D cuts --- swfiles/sw_fitpowder.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index e65935df..84852140 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -769,8 +769,9 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) 'This function is only valid for 2D data'); cut = struct('x', obj.ebin_cens, 'qmin', qmin, 'qmax', qmax); ikeep = obj.modQ_cens > qmin & obj.modQ_cens <= qmax; + ifinite = isfinite(obj.y(:, ikeep)); cut.y = mean(obj.y(:, ikeep), 2, 'omitnan'); - cut.e = sqrt(sum(obj.e(:, ikeep).^2, 2))/sum(ikeep); + cut.e = sqrt(sum(obj.e(:, ikeep).^2, 2))./sum(ifinite, 2); end function ycalc = rebin_powspec_to_1D_cuts(obj, ycalc) % sum up successive nQ points along |Q| axis (dim=2) From cc82c7fa487f50fb9526558e3f6a52478704049b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 7 Feb 2025 14:42:46 +0000 Subject: [PATCH 08/33] Add method to clear background region --- swfiles/sw_fitpowder.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 84852140..05e00860 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -430,8 +430,12 @@ function set_caching(obj, do_cache) function clear_cache(obj) obj.ycalc_cached = []; obj.model_params_cached = []; + obj.clear_background_region(); + end + + function clear_background_region(obj) obj.ibg = []; - obj.reset_errors_of_bg_bins() + obj.reset_errors_of_bg_bins(); end function set_bg_region(obj, en_lo, en_hi, varargin) From 91e673b55e6e305c214f63ef63e66ae0ea9e215a Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 17 Mar 2025 16:12:36 +0000 Subject: [PATCH 09/33] Support general polynomial background in Q and en --- swfiles/sw_fitpowder.m | 73 ++++++++++++++++++++++++++++++------------ 1 file changed, 53 insertions(+), 20 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 05e00860..4250cd11 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -161,15 +161,18 @@ end properties (SetAccess = private) - background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) - fbg_planar = @(en, modQ, slope_en, slope_modQ, intercept) slope_en*en(:) + slope_modQ*modQ(:)' + intercept; - fbg_indep = @(en, slope_en, intercept) slope_en*en(:) + intercept; + background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p) + fbg_planar = @(en, modQ, p, npoly_en) polyval([reshape(p(1:npoly_en),1,[]), 0],en(:)) + ... + polyval(p(npoly_en+1:end), modQ(:)'); % en^n,...en^1, q^n,..q^1, const + fbg_indep = @(en, p) polyval(p, en(:)); do_cache = true ycalc_cached = [] model_params_cached = [] liveplot_counter = 0 ibg = [] bg_errors = [] + npoly_modQ = 1; + npoly_en = 1; end properties (Constant) @@ -214,9 +217,38 @@ function initialise_background_parameters_and_bounds(obj) obj.bounds = [obj.bounds(1:obj.nparams_model,:); bg_bounds; obj.bounds(end, :)]; end + function set_bg_npoly_modQ(obj, npoly_modQ) + if obj.background_strategy == "independent" + error('sw_fitpowder:invalidinput', ... + 'set_bg_npoly_modQ only applies to planar background'); + end + if npoly_modQ < 0 + error('sw_fitpowder:invalidinput', ... + 'npoly_modQ must be positive int'); + end + if obj.ndim==1 && npoly_modQ > obj.ncuts - 1 + error('sw_fitpowder:invalidinput', ... + 'Not enough 1D cuts to fit polynomial of order npoly_modQ.'); + end + % zero-init bg parmas and nparams_bg + obj.npoly_modQ = int32(npoly_modQ); + obj.nparams_bg = obj.get_nparams_in_background_func(); + obj.initialise_background_parameters_and_bounds(); + end + + function set_bg_npoly_en(obj, npoly_en) + if npoly_en < 0 + error('sw_fitpowder:invalidinput', ... + 'npoly_en must be positive int'); + end + obj.npoly_en = int32(npoly_en); + obj.nparams_bg = obj.get_nparams_in_background_func(); + obj.initialise_background_parameters_and_bounds(); + end + function set_background_strategy(obj, strategy) if strategy == "planar" - obj.fbg = obj.fbg_planar; + obj.fbg = @(en, modQ, p) obj.fbg_planar(en, modQ, p, obj.npoly_en); elseif strategy == "independent" if obj.ndim == 2 error('sw_fitpowder:invalidinput', ... @@ -237,20 +269,22 @@ function set_background_strategy(obj, strategy) % set good guess for background pars modQs = obj.get_modQ_cens_of_cuts(); if obj.background_strategy == "independent" - % keep same energy slope (scaled by nQ) - obj.set_bg_parameters(1, bg_pars(1)); + % keep same energy bg + obj.set_bg_parameters(1:obj.npoly_en, bg_pars(1:obj.npoly_en)); % set intercept (zero energy) for reach cut separately - bg_pars = num2cell(bg_pars); - intercepts = obj.fbg_planar(0, modQs, bg_pars{:}); + intercepts = obj.fbg_planar(0, modQs, bg_pars, obj.npoly_en); for icut = 1:obj.ncuts - obj.set_bg_parameters(2, intercepts(icut), icut); + obj.set_bg_parameters(obj.npoly_en+1, intercepts(icut), icut); end elseif obj.background_strategy == "planar" - % average the energy slope - slope_en = mean(bg_pars(1:2:end)); - % fit intercepts to linear Q - pval = polyfit(modQs, bg_pars(2:2:end), 1); - obj.set_bg_parameters(1:3, [slope_en, pval]); + % reshape each col is the params for a cut + bg_pars = reshape(bg_pars, [], obj.ncuts); + % fit constant/intercepts to poly in Q + pval = polyfit(modQs, bg_pars(obj.npoly_en+1,:), obj.npoly_modQ); + obj.set_bg_parameters(obj.npoly_en+1:obj.nparams_bg, pval); + % average energy dep. params (only very rough) + bg_pars_en = mean(bg_pars(1:obj.npoly_en,:), 2); + obj.set_bg_parameters(1:obj.npoly_en, bg_pars_en); end end end @@ -405,16 +439,15 @@ function crop_q_range(obj, qmin, qmax) function bg = calc_background(obj, bg_params) % add background - bg_params = num2cell(bg_params); if obj.background_strategy == "planar" - bg = obj.fbg(obj.ebin_cens, obj.modQ_cens, bg_params{:}); + bg = obj.fbg(obj.ebin_cens, obj.modQ_cens, bg_params); elseif obj.ndim == 1 % add energy dependent background to each cut bg = zeros(size(obj.y)); istart = 1; for icut = 1:obj.ncuts iend = istart + obj.nparams_bg - 1; - bg(:,icut) = obj.fbg(obj.ebin_cens(:), bg_params{istart:iend}); + bg(:,icut) = obj.fbg(obj.ebin_cens(:), bg_params(istart:iend)); istart = iend + 1; end end @@ -694,7 +727,7 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) end function ax = plot_2d_data(obj, y) - ax = subplot(1,1,1); + ax = gca; box on; hold on; h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, y); h.AlphaData = double(obj.y > 0); % make empty bins transparent @@ -786,9 +819,9 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) function nbg_pars = get_nparams_in_background_func(obj) % get background parameters if obj.background_strategy == "planar" - nbg_pars = nargin(obj.fbg) - 2; % dependent variables energy, modQ + nbg_pars = obj.npoly_en + obj.npoly_modQ + 1; % dependent variables energy, modQ else - nbg_pars = nargin(obj.fbg) - 1; % dependent variable energy + nbg_pars = obj.npoly_en + 1; % dependent variable energy end end From 4464c7fe31176327fc954da94dd1efed2972fb57 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 17 Mar 2025 17:07:26 +0000 Subject: [PATCH 10/33] Add unit tests for setting polynomial order --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 73 ++++++++++++++++++- swfiles/sw_fitpowder.m | 5 +- 2 files changed, 73 insertions(+), 5 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 0ffceed5..28cc777a 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -33,7 +33,9 @@ function setup_spinw_obj_and_expected_result(testCase) -Inf Inf; -Inf Inf; 0 Inf], ... - 'ibg', []); + 'ibg', [], ... + 'npoly_modQ', 1, ... + 'npoly_en', 1); testCase.default_fields = fieldnames(testCase.default_fitpow); end end @@ -359,7 +361,7 @@ function test_fit_background(testCase, fit_params) out.y(1) = 10; % higher so other bins are background out.fix_bg_parameters(1:2); % fix slopes of background to 0 out.set_bg_parameters(3, 1.5); % initial guess - out.fit_background(fit_params{:}) + out.fit_background(fit_params{:}); expected_fitpow = testCase.default_fitpow; expected_fitpow.y(1) = 10; expected_fitpow.ibg = [3;6;2;5;4]; @@ -516,6 +518,73 @@ function test_set_bg_region_data_1d(testCase) testCase.verify_results(out, expected_fitpow); end + function test_set_npoly_modQ(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.set_bg_npoly_modQ(2); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.npoly_modQ = 2; + expected_fitpow.params = [expected_fitpow.params(1:end-1); + 0; + expected_fitpow.params(end)]; + expected_fitpow.bounds = [expected_fitpow.bounds(1:end-1,:); + -Inf, Inf; + expected_fitpow.bounds(end,:)]; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_npoly_en(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.set_bg_npoly_en(2); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.npoly_en = 2; + % add extra row in params and bounds + expected_fitpow.params = expected_fitpow.params([1:2,2:end]); + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + testCase.verify_results(out, expected_fitpow); + end + + function test_set_npoly_modQ_1d_data_indep_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent"); + testCase.verifyError(... + @() out.set_bg_npoly_modQ(2), ... + 'sw_fitpowder:invalidinput'); + end + + function test_set_npoly_en_1d_data_indep_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent"); + out.set_bg_npoly_en(2) + expected_fitpow = testCase.default_fitpow; + expected_fitpow.npoly_en = 2; + expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + % add 3 extra rows: + % 3 params x 2 cuts vs 3 planar parmas order=1 + expected_fitpow.params = expected_fitpow.params([1:4,2:end]); + expected_fitpow.bounds = expected_fitpow.bounds([1:4,2:end],:); + testCase.verify_results(out, expected_fitpow); + end + + function test_set_npoly_modQ_1d_data_planar_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "planar"); + % try adding order larger than numebr cuts + testCase.verifyError(... + @() out.set_bg_npoly_modQ(2), ... + 'sw_fitpowder:invalidinput'); + % set it to constant in modQ + out.set_bg_npoly_modQ(0) + expected_fitpow = testCase.default_fitpow; + expected_fitpow.npoly_modQ = 0; + expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + % remove a row + expected_fitpow.params = expected_fitpow.params([1,3:end]); + expected_fitpow.bounds = expected_fitpow.bounds([1,3:end],:); + testCase.verify_results(out, expected_fitpow); + end + end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 4250cd11..69645a70 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -230,8 +230,7 @@ function set_bg_npoly_modQ(obj, npoly_modQ) error('sw_fitpowder:invalidinput', ... 'Not enough 1D cuts to fit polynomial of order npoly_modQ.'); end - % zero-init bg parmas and nparams_bg - obj.npoly_modQ = int32(npoly_modQ); + obj.npoly_modQ = round(npoly_modQ); obj.nparams_bg = obj.get_nparams_in_background_func(); obj.initialise_background_parameters_and_bounds(); end @@ -241,7 +240,7 @@ function set_bg_npoly_en(obj, npoly_en) error('sw_fitpowder:invalidinput', ... 'npoly_en must be positive int'); end - obj.npoly_en = int32(npoly_en); + obj.npoly_en = round(npoly_en); obj.nparams_bg = obj.get_nparams_in_background_func(); obj.initialise_background_parameters_and_bounds(); end From 0ad70f084a3eda247a77b4aa6fc11f58c57c0818 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 28 Apr 2025 18:10:26 +0100 Subject: [PATCH 11/33] Use 2D data Q bins when replace 1D cuts Also use same method to cut data in plotting 1D cuts of 2D data as when replacing data with 1D cuts. Updated tests --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 6 +- swfiles/sw_fitpowder.m | 57 +++++++++++++------ 2 files changed, 44 insertions(+), 19 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 28cc777a..ce6a9e3d 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -121,7 +121,7 @@ function test_replace_2D_data_with_1D_cuts(testCase) qcens = [4, 5]; out.replace_2D_data_with_1D_cuts(qcens-0.5, qcens+0.5) expected_fitpow = testCase.default_fitpow; - expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + expected_fitpow.modQ_cens = qcens; testCase.verify_results(out, expected_fitpow); end @@ -132,7 +132,7 @@ function test_replace_2D_data_with_1D_cuts_specify_bg(testCase) out.replace_2D_data_with_1D_cuts(qcens-0.5, qcens+0.5,... "independent") expected_fitpow = testCase.default_fitpow; - expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; + expected_fitpow.modQ_cens = qcens; % not nQ values as using bins in 2D data % add extra background param expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); @@ -406,7 +406,7 @@ function test_calc_cost_func_of_background_planar(testCase) bg_pars(end) = 1; out.estimate_constant_background(); % so ibg is set cost = out.calc_cost_func_of_background(bg_pars); - testCase.verify_val(cost, 10); + testCase.verify_val(cost, 5); end function test_set_errors_of_bg_bins(testCase) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 69645a70..35ad27c4 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -138,6 +138,7 @@ ebin_cens modQ_cens nQ = 10 % number |Q| bins to calc. in integration limits of 1D cut + modQ_icuts ndim ncuts = 0; % spinw funtion inputs @@ -395,7 +396,13 @@ function add_data(obj, data) obj.y = [obj.y cut.y(:)]; obj.e = [obj.e cut.e(:)]; dQ = (cut.qmax - cut.qmin)/obj.nQ; - obj.modQ_cens = [obj.modQ_cens, (cut.qmin+dQ/2):dQ:cut.qmax]; + if isfield(cut, 'qs') + obj.modQ_cens = [obj.modQ_cens, cut.qs(:)']; + obj.modQ_icuts = [obj.modQ_icuts, ones(1, numel(cut.qs))]; + else + obj.modQ_cens = [obj.modQ_cens, ((cut.qmin+dQ/2):dQ:cut.qmax)]; + obj.modQ_icuts = [obj.modQ_icuts, ones(1, obj.nQ)]; + end end end obj.powspec_args.Evect = obj.ebin_cens(:)'; % row vect @@ -715,7 +722,7 @@ function plot_1d_cuts_on_data(obj, ycalc, varargin) ylim(ax, [ymin, ymax]); xlabel(ax, 'Energy (meV)') ylabel(ax, 'Intensity'); - title(ax, num2str(modQs(icut), 2) + " $\AA^{-1}$", 'interpreter','latex') + title(ax, num2str(modQs(icut), "%.2f") + " $\AA^{-1}$", 'interpreter','latex') end end @@ -762,22 +769,23 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) figure("color","white"); ncuts = numel(qmins); for icut = 1:ncuts - ikeep = obj.modQ_cens > qmins(icut) & obj.modQ_cens <= qmaxs(icut); ax = subplot(1, ncuts, icut); hold on; box on; - ycut = sum(obj.y(:,ikeep), 2); - plot(ax, obj.ebin_cens, ycut, 'ok'); if ~isempty(ycalc) - plot(ax, obj.ebin_cens, sum(ycalc(:,ikeep), 2), '-r'); + cut = obj.cut_2d_data(qmins(icut), qmaxs(icut), ycalc); + plot(ax, cut.x, cut.ycalc, '-r'); + else + cut = obj.cut_2d_data(qmins(icut), qmaxs(icut)); end + plot(ax, cut.x, cut.y, 'ok'); % calc xlims - ifinite = isfinite(ycut); + ifinite = isfinite(cut.y); istart = find(ifinite, 1, 'first'); iend = find(ifinite, 1, 'last'); - xlim(ax, [obj.ebin_cens(istart), obj.ebin_cens(iend)]); + xlim(ax, [cut.x(istart), cut.x(iend)]); xlabel(ax, 'Energy (meV)') ylabel(ax, 'Intensity'); - title(ax, num2str(0.5*(qmaxs(icut)+qmins(icut)), 2) + " $\AA^{-1}$", 'interpreter','latex') + title(ax, num2str(0.5*(cut.qmin +cut.qmax), "%.2f") + " $\AA^{-1}$", 'interpreter','latex') end end @@ -799,20 +807,37 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) modQ_cens = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); end - function cut = cut_2d_data(obj, qmin, qmax) + function cut = cut_2d_data(obj, qmin, qmax, ycalc) assert(obj.ndim ==2, ... 'sw_fitpowder:invalidinput', ... 'This function is only valid for 2D data'); - cut = struct('x', obj.ebin_cens, 'qmin', qmin, 'qmax', qmax); ikeep = obj.modQ_cens > qmin & obj.modQ_cens <= qmax; + cut = struct('x', obj.ebin_cens, 'qmin', qmin, 'qmax', qmax, ... + 'qs', obj.modQ_cens(ikeep)); ifinite = isfinite(obj.y(:, ikeep)); cut.y = mean(obj.y(:, ikeep), 2, 'omitnan'); cut.e = sqrt(sum(obj.e(:, ikeep).^2, 2))./sum(ifinite, 2); - end - function ycalc = rebin_powspec_to_1D_cuts(obj, ycalc) - % sum up successive nQ points along |Q| axis (dim=2) - ycalc = reshape(ycalc, size(ycalc,1), obj.nQ, []); - ycalc = squeeze(mean(ycalc, 2, 'omitnan')); + if nargin == 4 && ~isempty(ycalc) + % also cut ycalc + cut.ycalc = mean(ycalc(:, ikeep), 2, 'omitnan'); + end + end + function ycalc_1d = rebin_powspec_to_1D_cuts(obj, ycalc) + % check if regular nQ points for each cut + counts = groupcounts(obj.modQ_icuts); + unique_counts = unique(counts); + if numel(unique_counts) == 1 && unique_counts == obj.nQ + % avg successive nQ points along |Q| axis (dim=2) + ycalc_1d = reshape(ycalc, size(ycalc,1), obj.nQ, []); + ycalc_1d = squeeze(mean(ycalc, 2, 'omitnan')); + else + % different num Q points per cut, loop over cuts + ycalc_1d = zeros(size(ycalc)); + for icut = 1:obj.ncuts + ikeep = obj.modQ_icuts == icut; + ycalc_1d(:, icut) = mean(ycalc(:, ikeep), 2, 'omitnan'); + end + end end function nbg_pars = get_nparams_in_background_func(obj) From df70190fe870a2d23b39f300a9fa4d68ba779631 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 29 Apr 2025 17:49:48 +0100 Subject: [PATCH 12/33] Set SigP to 0 for fixed parameters in lm4 Now pfit and SigP have same shape --- swfiles/+ndbase/lm4.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/swfiles/+ndbase/lm4.m b/swfiles/+ndbase/lm4.m index ace81059..c31cf2de 100644 --- a/swfiles/+ndbase/lm4.m +++ b/swfiles/+ndbase/lm4.m @@ -223,14 +223,14 @@ end % collect output pOpt = cost_func_wrap.get_bound_parameters(p); + perr = zeros(size(pOpt)); fVal = cost_val / ndof; if exit_flag > 0 % converged on solution - calculate errors cov = pinv(hess) * 2.0 * fVal; - perr = sqrt(diag(cov)); + perr(cost_func_wrap.ifree) = sqrt(diag(cov)); else message = "Failed to converge in MaxIter"; - perr = zeros(size(pOpt)); end end From 40eca62c82a6222987083c46d50b5ddeda6a327b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 29 Apr 2025 17:51:03 +0100 Subject: [PATCH 13/33] Fix bug for 1D cuts now using 2D Q bins --- swfiles/sw_fitpowder.m | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 35ad27c4..f864c88c 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -389,7 +389,8 @@ function add_data(obj, data) obj.ndim = 1; obj.ebin_cens = data(1).x; obj.ncuts = numel(data); - for cut = data + for icut = 1:obj.ncuts + cut = data(icut); assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... 'sw_fitpowder:invalidinput', ... 'Input cell does not have correct fields'); @@ -398,10 +399,10 @@ function add_data(obj, data) dQ = (cut.qmax - cut.qmin)/obj.nQ; if isfield(cut, 'qs') obj.modQ_cens = [obj.modQ_cens, cut.qs(:)']; - obj.modQ_icuts = [obj.modQ_icuts, ones(1, numel(cut.qs))]; + obj.modQ_icuts = [obj.modQ_icuts, icut*ones(1, numel(cut.qs))]; else obj.modQ_cens = [obj.modQ_cens, ((cut.qmin+dQ/2):dQ:cut.qmax)]; - obj.modQ_icuts = [obj.modQ_icuts, ones(1, obj.nQ)]; + obj.modQ_icuts = [obj.modQ_icuts, icut*ones(1, obj.nQ)]; end end end @@ -824,7 +825,7 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) end function ycalc_1d = rebin_powspec_to_1D_cuts(obj, ycalc) % check if regular nQ points for each cut - counts = groupcounts(obj.modQ_icuts); + counts = groupcounts(obj.modQ_icuts(:)); unique_counts = unique(counts); if numel(unique_counts) == 1 && unique_counts == obj.nQ % avg successive nQ points along |Q| axis (dim=2) @@ -832,7 +833,7 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) ycalc_1d = squeeze(mean(ycalc, 2, 'omitnan')); else % different num Q points per cut, loop over cuts - ycalc_1d = zeros(size(ycalc)); + ycalc_1d = zeros(size(ycalc, 1), obj.ncuts); for icut = 1:obj.ncuts ikeep = obj.modQ_icuts == icut; ycalc_1d(:, icut) = mean(ycalc(:, ikeep), 2, 'omitnan'); From 7a6a70e25edc345d67d590596795c23ac1f8146c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 29 Apr 2025 18:41:37 +0100 Subject: [PATCH 14/33] Update tutorial to include setting bg strategy and setting bg region --- tutorials/publish/tutorial39.m | 190 ++++++++++++++++++++++++--------- 1 file changed, 140 insertions(+), 50 deletions(-) diff --git a/tutorials/publish/tutorial39.m b/tutorials/publish/tutorial39.m index b52e514b..99c0b272 100644 --- a/tutorials/publish/tutorial39.m +++ b/tutorials/publish/tutorial39.m @@ -40,7 +40,7 @@ Q = linspace(0.1,3,60); spec = tri.powspec(Q,'Evect',E,'nRand',2e2); % scale data and add background -spec.swConv= scale*(spec.swConv + 0.1*mean(spec.swConv, 'all')); +spec.swConv= scale*(spec.swConv + 0.1*mean(spec.swConv, "all", 'omitnan')); % crop to kinematic range of instrument and apply some arbitrary resolution spec = sw_instrument(spec, 'dQ', dQ, 'dE', eres, 'ThetaMin', 3.5, 'Ei', Ei); % add noise @@ -67,11 +67,6 @@ % The class will call other spinw algorithms such as powspec and % sw_instrument. The parameters used in these calls can be set by the user % by modifying structs as shown below. -% -% Data can be cropped in |Q| and energy transfer - for example it is -% advisible not to include the low-energy transfer region close to the -% elastic line (relative to the resolution). - fit_func = @(obj, p) matparser(obj, 'param', p, 'mat', {'J_1'}, 'init', true); J1 = 0.85; @@ -83,26 +78,91 @@ fitpow.powspec_args.dE = eres; % emergy resolution fitpow.powspec_args.fastmode = true; fitpow.powspec_args.neutron_output = true; -fitpow.powspec_args.nRand = 2e2; % low for speed (typically want > 1e3) +fitpow.powspec_args.nRand = 1e3; % low for speed (typically want > 1e3) fitpow.powspec_args.hermit = true; % set parameters passsed to sw_instrument fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei); -% crop data -fitpow.crop_energy_range(1.5*eres(0), inf); % typically want to avoid elastic line -fitpow.crop_q_range(0.25, 3); -fitpow.powspec_args.dE = eres; -fitpow.powspec_args.fastmode = true; -fitpow.powspec_args.neutron_output = true; -% estimate background and scale factor +%% Estimate the background automatically +% Before fitting the data it is necessary to first get a good estimate for +% the background and the scale factor. +% There is a means or trying to automatically determine the background +% region - it is worth trying this first + fitpow.estimate_constant_background(); +fitpow.estimate_scale_factor() + +% view the background region on 2D colorfill plot of data +figure("color","white") +fitpow.plot_2d_data(fitpow.y); +fitpow.plot_background_region(); + +% view 1D cuts of data and the initial guess +qcens = 0.75:0.75:2.25; +dq = 0.1; +fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, fitpow.params); + +%% User supplied background region +% If the automatic background estimation doesn't look reasonable then it is +% possible to define regions of |Q| and energy transfer that are background +% and use these for fitting. + +% set background regions +%----------------------- + +fitpow.clear_background_region(); % in case previously set + +% set background region with en > 3.75 meV +fitpow.set_bg_region(3.75, inf); +% set background region with +% 0.65 < |Q| < 0.9 Ang^-1 +% 0.8 < en < 1.2 meV +fitpow.set_bg_region(0.8, 1.2, 0.65, 0.9); + +% view the background region on 2D colorfill plot of data +figure("color","white") +fitpow.plot_2d_data(fitpow.y); +fitpow.plot_background_region(); + +% fit background (to voxels in background region) +%----------------------------------------------------- + +% fit only a constant/flat background +% default background is planar with parameters: +% [slope_energy, slope_modQ, intercept] +fitpow.fix_bg_parameters(1:2); % fix slopes to initial value (0) +fitpow.set_bg_parameter_bounds(3, 0.0, []) % Set lb of constant bg = 0 + +[pfit, ~, stat] = fitpow.fit_background(); fitpow.estimate_scale_factor(); +% view 1D cuts of data and the initial guess +fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, fitpow.params); + + +%% Cropping data +% Data can be cropped in |Q| and energy transfer - for example it is +% advisible not to include the low-energy transfer region close to the +% elastic line (relative to the resolution). + +% crop data +fitpow.crop_energy_range(1.5*eres(0), inf); % typically want to avoid elastic line +fitpow.crop_q_range(0.25, 3); -% plot data to inspect - any arguments after parameters are passed to -% contour plot (which shows the result the calculated intensity) -fitpow.plot_result(fitpow.params, 'EdgeAlpha', 0.9, 'LineWidth', 1 ,'EdgeColor', 'k'); +% plot the cropped data and the inital simulated spectrum +fig = figure("color","white"); hold all; +fig.Position(3) = 1.5*fig.Position(3); +subplot(1,2,1) +ax_dat = fitpow.plot_2d_data(fitpow.y); +subplot(1,2,2) +[ycalc, ~] = fitpow.calc_spinwave_spec(fitpow.params); +ax_calc = fitpow.plot_2d_data(ycalc); +ax_calc.CLim = ax_dat.CLim; +for ax = [ax_dat, ax_calc] + ax.YLim(1) = 1.5*eres(0); +end +sgtitle(fig, "Initial guess") %% Fit and inspect the result @@ -114,20 +174,14 @@ % The cost function minimised in the fit can be chi-squared ("chisq") or % the unweighted sum of the squared residuals ("Rsq"). Also the minimiser % can be set to any callable as long as it shares the API of the ndbase -% minimisers. Typically users would be recommended ndbase.lm or +% minimisers. Typically users would be recommended to use ndbase.lm4 or % ndbase.simplex. % set cost function and minimser fitpow.cost_function = "Rsq"; % "Rsq" or "chisq" fitpow.optimizer = @ndbase.simplex; % or e.g. ndbase.simplex -% Fix background parameters (planar background of form -% default background is planar with parameters: -% [slope_energy, slope_modQ, intercept] -fitpow.fix_bg_parameters(1:2); % fix slopes to initial value (0) -fitpow.set_bg_parameter_bounds(3, 0.0, []) % Set lb of constant bg = 0 -% fit background -fitpow.fit_background(); +% Fix the constant background at the result of the previous fit fitpow.fix_bg_parameters(3); % fix fitted constant bg % set bounds on first (and only in this case) model parameter @@ -139,8 +193,6 @@ fitpow.plot_result(pfit, 10, 'EdgeAlpha', 0.9, 'LineWidth', 1 ,'EdgeColor', 'k') % make 1D plots -qcens = [0.8:0.4:1.6]; -dq = 0.05; fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, pfit) @@ -150,35 +202,73 @@ % Alternatively, having loaded in 2D data the sw_fitpowder class has a % convenient method to take 1D cuts and replace the data. % -% To evaluate the cost function for 1D cuts, the class averages the model -% over nQ points. Depending on the size of the 2D data and the value of nQ -% it can be much quicker to fit 1D cuts. -% +% If 2D data are replaced with 1D cuts then the class evaluates the +% cost function using the average of the model over the |Q| bins in the +% 2D data within the integration limits. +% If the class is initialised with 1D cuts then the class averages the model +% over nQ points. + +% set inital parameters to those of 2D fit +fitpow.params = pfit(:); + +% set data to 1D cuts +fitpow.replace_2D_data_with_1D_cuts(qcens-dq, qcens+dq); + +% fit data with lm4 minimiser (which calculates parameter uncertainties) +% this requires a larger nRand as derivatives in LM algorithm sensitive to +% noise introduced in powder averaging. +% if uncertainties come out imaginary try increasing diff_step +fitpow.optimizer = @ndbase.lm4; +fitpow.powspec_args.nRand = 2e4; +fit_kwargs = {'resid_handle', true, 'diff_step', 1e-3}; + +[pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); + +% display result +par_names = ["J1", "SlopeEn", "SlopeQ", "Const", "Scale"]; +["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] + +% "Name" "Value" "Error" +% "J1" "1.004828" "0.00432656" +% "SlopeEn" "0" "0" +% "SlopeQ" "0" "0" +% "Const" "0.93551168" "0" +% "Scale" "1018.0549" "12.6781" + +fitpow.plot_result(pfit) + +%% Change background strategy % The background for 1D cuts can be handled in two ways: % % - A planar background as in the 2D case % % - A linear background in energy transfer that's independent for each cut. % -% The background parameters and bounds are reset if the background -% strategy is changed to independent. - -% set nQ before replacing data (reduced from default, 10, for speed) -fitpow.nQ = 5; -% change background strategyto independent (planar is default) -fitpow.replace_2D_data_with_1D_cuts(qcens-dq, qcens+dq, "independent"); - -% independent linear background for all cuts with parameters: -% [slope_energy, intercept] -fitpow.fix_bg_parameters(1); % fix slopes to initial value (0) for all cuts -% e.g. fitpow.fix_bg_parameters(1, 2) would fix the bg slope of the second cut -fitpow.set_bg_parameter_bounds(2, 0.0, []) % Set lb of constant bg = 0 -% fit background -fitpow.estimate_constant_background(); -fitpow.fit_background(); -fitpow.fix_bg_parameters(2); % fix fitted constant bg - -[pfit, cost_val, stat] = fitpow.fit(); +% The default is planar - but this can be changed after the class has been +% initialised. +% The class will attempt to get reasonable guesses for the +% new background parameters, but note the background parameter bounds will +% be reset. + +fitpow.set_background_strategy("independent"); + +[pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); + +% display result +bg_names = repmat(["SlopeEn_cut", "Const_cut"]', fitpow.ncuts, 1); +bg_names = bg_names + repelem([1:fitpow.ncuts]', fitpow.nparams_bg,1); +par_names = ["J1", bg_names', "Scale"]; +["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] + +% "Name" "Value" "Error" +% "J1" "1.0072907" "0.00407115" +% "SlopeEn_cut1" "-1.1362877" "0.400225" +% "Const_cut1" "2.056183" "0.949017" +% "SlopeEn_cut2" "1.7884139" "0.460254" +% "Const_cut2" "-6.8161423" "1.58929" +% "SlopeEn_cut3" "0.53730869" "0.372718" +% "Const_cut3" "-3.3143013" "0.994353" +% "Scale" "1193.6718" "32.4533" fitpow.plot_result(pfit) From 7b470aa3f7f35a0ccfd76dde57dcd2bc11018f85 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 16 May 2025 14:17:41 +0100 Subject: [PATCH 15/33] Add method to export data in struct form --- swfiles/sw_fitpowder.m | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index f864c88c..8f1ce4a7 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -361,6 +361,25 @@ function set_scale_bounds(obj, lb, ub) obj.set_bounds(size(obj.bounds,1), lb, ub); end + function data = export_data(obj) + if obj.ndim == 2 + data.y = obj.y'; + data.e = obj.e'; + data.x = {obj.ebin_cens, obj.modQ_cens}; + else + % 1D - array of structs returned + data = repelem(struct('x', obj.ebin_cens), obj.ncuts); + for icut = 1:obj.ncuts + data(icut).y = obj.y(:, icut); + data(icut).e = obj.e(:, icut); + % find q limits + data(icut).qs = obj.modQ_cens(obj.modQ_icuts == icut); + data(icut).qmin = data(icut).qs(1); + data(icut).qmax = data(icut).qs(end); + end + end + end + function add_data(obj, data) if ~isa(data, "struct") data = arrayfun(@obj.convert_horace_to_struct, data); From f8e4567f03ca7d8d44a310d046d4f5773641e25b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 19 May 2025 15:58:03 +0100 Subject: [PATCH 16/33] Fix rebin_powspec_to_1D_cuts bug causing test failure --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 2 +- swfiles/sw_fitpowder.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index ce6a9e3d..ccba3a27 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -406,7 +406,7 @@ function test_calc_cost_func_of_background_planar(testCase) bg_pars(end) = 1; out.estimate_constant_background(); % so ibg is set cost = out.calc_cost_func_of_background(bg_pars); - testCase.verify_val(cost, 5); + testCase.verify_val(cost, 10); end function test_set_errors_of_bg_bins(testCase) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 8f1ce4a7..db72e3f8 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -849,7 +849,7 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) if numel(unique_counts) == 1 && unique_counts == obj.nQ % avg successive nQ points along |Q| axis (dim=2) ycalc_1d = reshape(ycalc, size(ycalc,1), obj.nQ, []); - ycalc_1d = squeeze(mean(ycalc, 2, 'omitnan')); + ycalc_1d = squeeze(mean(ycalc_1d, 2, 'omitnan')); else % different num Q points per cut, loop over cuts ycalc_1d = zeros(size(ycalc, 1), obj.ncuts); From 17e2b1c5c0b0185d7d25ead56174dac85039bf55 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 19 May 2025 16:09:48 +0100 Subject: [PATCH 17/33] Add unit tests for export_data method --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 15 +++++++++++++++ swfiles/sw_fitpowder.m | 12 ++++++------ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index ccba3a27..68fae591 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -585,6 +585,21 @@ function test_set_npoly_modQ_1d_data_planar_bg(testCase) testCase.verify_results(out, expected_fitpow); end + function test_export_data_2d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + data = out.export_data(); + testCase.verify_val(data, testCase.data_2d); + end + + function test_export_data_1d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1); + data = out.export_data(); + testCase.verify_val(data, testCase.data_1d_cuts); + end + + end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index db72e3f8..4cf34455 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -370,12 +370,12 @@ function set_scale_bounds(obj, lb, ub) % 1D - array of structs returned data = repelem(struct('x', obj.ebin_cens), obj.ncuts); for icut = 1:obj.ncuts - data(icut).y = obj.y(:, icut); - data(icut).e = obj.e(:, icut); - % find q limits - data(icut).qs = obj.modQ_cens(obj.modQ_icuts == icut); - data(icut).qmin = data(icut).qs(1); - data(icut).qmax = data(icut).qs(end); + data(icut).y = obj.y(:, icut)'; + data(icut).e = obj.e(:, icut)'; + % find q limits (bin edges) + qs = obj.modQ_cens(obj.modQ_icuts == icut); + data(icut).qmin = qs(1) - 0.5*(qs(2)-qs(1)); + data(icut).qmax = qs(end) + 0.5*(qs(end)-qs(end-1)); end end end From 484018d37b6e354c30882d02b94e5fc822fcf2fd Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 19 May 2025 16:10:30 +0100 Subject: [PATCH 18/33] Update tutorial with changing polynomial order --- tutorials/publish/tutorial39.m | 64 +++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 29 deletions(-) diff --git a/tutorials/publish/tutorial39.m b/tutorials/publish/tutorial39.m index 99c0b272..e582a8d4 100644 --- a/tutorials/publish/tutorial39.m +++ b/tutorials/publish/tutorial39.m @@ -120,10 +120,6 @@ % 0.8 < en < 1.2 meV fitpow.set_bg_region(0.8, 1.2, 0.65, 0.9); -% view the background region on 2D colorfill plot of data -figure("color","white") -fitpow.plot_2d_data(fitpow.y); -fitpow.plot_background_region(); % fit background (to voxels in background region) %----------------------------------------------------- @@ -195,7 +191,6 @@ % make 1D plots fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, pfit) - %% Fit 1D cuts % As discussed the sw_fitpowder class can be initialised with a cell array % of xye structs correpsonding to 1D cuts at constant |Q|. @@ -214,28 +209,41 @@ % set data to 1D cuts fitpow.replace_2D_data_with_1D_cuts(qcens-dq, qcens+dq); +% change background to quadratic in |Q| +fitpow.set_bg_npoly_modQ(2); % zeros previously optimised background +% Fix linear terms of background +% bg params ordered as en^n,...en^1, q^n,..q^1, const +fitpow.fix_bg_parameters(3); % fix term ~|Q| +fitpow.fix_bg_parameters(1); % fix term ~energy +% fit background and scale (fix model params) to all bins (not just bg +% region -0 only advisable if near good solution) +[pfit, ~, stat] = fitpow.fit_background_and_scale(); +% fitpow.plot_result(fitpow.params) + % fit data with lm4 minimiser (which calculates parameter uncertainties) % this requires a larger nRand as derivatives in LM algorithm sensitive to % noise introduced in powder averaging. % if uncertainties come out imaginary try increasing diff_step fitpow.optimizer = @ndbase.lm4; -fitpow.powspec_args.nRand = 2e4; +fitpow.powspec_args.nRand = 1e4; fit_kwargs = {'resid_handle', true, 'diff_step', 1e-3}; [pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); +fitpow.plot_result(pfit); + % display result -par_names = ["J1", "SlopeEn", "SlopeQ", "Const", "Scale"]; +par_names = ["J1", "Bg_En1", "Bg_Q2", "Bg_Q1", "Bg_Q0", "Scale"]; ["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] -% "Name" "Value" "Error" -% "J1" "1.004828" "0.00432656" -% "SlopeEn" "0" "0" -% "SlopeQ" "0" "0" -% "Const" "0.93551168" "0" -% "Scale" "1018.0549" "12.6781" +% "Name" "Value" "Error" +% "J1" "1.0027473" "0.00360547" +% "Bg_En1" "0" "0" +% "Bg_Q2" "0.02066691" "0.095079" +% "Bg_Q1" "0" "0" +% "Bg_Q0" "-0.036883773" "0.37194" +% "Scale" "1082.117" "19.405" -fitpow.plot_result(pfit) %% Change background strategy % The background for 1D cuts can be handled in two ways: @@ -250,27 +258,25 @@ % new background parameters, but note the background parameter bounds will % be reset. +fitpow.params = pfit(:); fitpow.set_background_strategy("independent"); [pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); +fitpow.plot_result(pfit) + % display result -bg_names = repmat(["SlopeEn_cut", "Const_cut"]', fitpow.ncuts, 1); +bg_names = repmat(["Bg_En1_cut", "Bg_En0_cut"]', fitpow.ncuts, 1); bg_names = bg_names + repelem([1:fitpow.ncuts]', fitpow.nparams_bg,1); par_names = ["J1", bg_names', "Scale"]; ["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] -% "Name" "Value" "Error" -% "J1" "1.0072907" "0.00407115" -% "SlopeEn_cut1" "-1.1362877" "0.400225" -% "Const_cut1" "2.056183" "0.949017" -% "SlopeEn_cut2" "1.7884139" "0.460254" -% "Const_cut2" "-6.8161423" "1.58929" -% "SlopeEn_cut3" "0.53730869" "0.372718" -% "Const_cut3" "-3.3143013" "0.994353" -% "Scale" "1193.6718" "32.4533" - -fitpow.plot_result(pfit) - - - +% "Name" "Value" "Error" +% "J1" "1.0059739" "0.00403765" +% "Bg_En1_cut1" "-0.85631989" "0.376039" +% "Bg_En0_cut1" "1.7546252" "0.892913" +% "Bg_En1_cut2" "0.52590868" "0.433608" +% "Bg_En0_cut2" "-2.3563309" "1.49539" +% "Bg_En1_cut3" "-0.41373845" "0.350704" +% "Bg_En0_cut3" "0.52287368" "0.936986" +% "Scale" "1131.9098" "30.4495" From e934108ffd8a4c96328fe7a4ff3b9aaeecb5470c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 2 Jun 2025 16:00:28 +0100 Subject: [PATCH 19/33] Swap order of bg parameters (Q first) and add bg param labels This is done so that there is always an E0 label regardless of background strategy --- swfiles/sw_fitpowder.m | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 4cf34455..a49cce55 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -163,8 +163,8 @@ properties (SetAccess = private) background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p) - fbg_planar = @(en, modQ, p, npoly_en) polyval([reshape(p(1:npoly_en),1,[]), 0],en(:)) + ... - polyval(p(npoly_en+1:end), modQ(:)'); % en^n,...en^1, q^n,..q^1, const + fbg_planar = @(en, modQ, p, npoly_modQ) polyval(p(npoly_modQ+1:end), en(:)) + ... + polyval([reshape(p(1:npoly_modQ),1,[]), 0], modQ(:)'); % q^n,..q^1, en^n,...en^1, const fbg_indep = @(en, p) polyval(p, en(:)); do_cache = true ycalc_cached = [] @@ -174,6 +174,7 @@ bg_errors = [] npoly_modQ = 1; npoly_en = 1; + bg_param_labels = []; end properties (Constant) @@ -216,6 +217,9 @@ function initialise_background_parameters_and_bounds(obj) % insert in middle of params and bound obj.params = [obj.params(1:obj.nparams_model); bg_params; obj.params(end)]; obj.bounds = [obj.bounds(1:obj.nparams_model,:); bg_bounds; obj.bounds(end, :)]; + % set labels + obj.bg_param_labels = [repmat("Q", obj.npoly_modQ, 1) + [obj.npoly_modQ:-1:1]'; + repmat("E", obj.npoly_en+1, 1) + [obj.npoly_en:-1:0]']; end function set_bg_npoly_modQ(obj, npoly_modQ) @@ -260,19 +264,19 @@ function set_background_strategy(obj, strategy) % store previous bg parameters before reset size of param array ibg_par = (obj.nparams_model+1):(numel(obj.params)-1); bg_pars = obj.params(ibg_par); - % reset param array d bounds etc. + % reset param array bounds etc. obj.background_strategy = strategy; obj.nparams_bg = obj.get_nparams_in_background_func(); if ~isempty(obj.params) obj.initialise_background_parameters_and_bounds(); if any(abs(bg_pars) > obj.zero_abs_tol) && obj.ndim == 1 - % set good guess for background pars + % attempt good guess for background pars modQs = obj.get_modQ_cens_of_cuts(); if obj.background_strategy == "independent" - % keep same energy bg - obj.set_bg_parameters(1:obj.npoly_en, bg_pars(1:obj.npoly_en)); + % keep same energy bg (excl const bg) + obj.set_bg_parameters(1:obj.npoly_en, bg_pars(obj.npoly_modQ+1:end-1)); % set intercept (zero energy) for reach cut separately - intercepts = obj.fbg_planar(0, modQs, bg_pars, obj.npoly_en); + intercepts = obj.fbg_planar(0, modQs, bg_pars, obj.npoly_modQ); for icut = 1:obj.ncuts obj.set_bg_parameters(obj.npoly_en+1, intercepts(icut), icut); end @@ -280,11 +284,12 @@ function set_background_strategy(obj, strategy) % reshape each col is the params for a cut bg_pars = reshape(bg_pars, [], obj.ncuts); % fit constant/intercepts to poly in Q - pval = polyfit(modQs, bg_pars(obj.npoly_en+1,:), obj.npoly_modQ); - obj.set_bg_parameters(obj.npoly_en+1:obj.nparams_bg, pval); + pval = polyfit(modQs, bg_pars(end,:), obj.npoly_modQ); + obj.set_bg_parameters(1:obj.npoly_modQ, pval(1:end-1)); + obj.set_bg_parameters(obj.nparams_bg, pval(end)); % const bg % average energy dep. params (only very rough) bg_pars_en = mean(bg_pars(1:obj.npoly_en,:), 2); - obj.set_bg_parameters(1:obj.npoly_en, bg_pars_en); + obj.set_bg_parameters(obj.npoly_modQ+1:obj.nparams_bg-1, bg_pars_en); end end end From 12657f66e4511daecb93e7ec6885d20764b2d611 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 2 Jun 2025 16:07:18 +0100 Subject: [PATCH 20/33] Update changing bg strategy tests with non-identical bg parameters --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 68fae591..3bfcb76e 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -82,13 +82,13 @@ function test_set_background_strategy_to_planar(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent"); % set some background parameters (correpsonding to planar bg - % with slope_en=slope_q=intercept =2 - out.set_bg_parameters(1, 2); % en_slope = 2 - out.set_bg_parameters(2, 10, 1); % intercept = 10 for cut 1 - out.set_bg_parameters(2, 12, 2); % intercept = 12 for cut 2 + % with slope_en=1, slope_q=2, intercept = 3 + out.set_bg_parameters(1, 1); % en_slope = 1 + out.set_bg_parameters(2, 11, 1); % intercept = 4*2 + 3 for cut 1 + out.set_bg_parameters(2, 13, 2); % intercept = 5*2 + 3 for cut 2 out.set_background_strategy("planar"); expected_fitpow = testCase.default_fitpow; - expected_fitpow.params(2:end-1) = 2; + expected_fitpow.params(2:end-1) = [2,1,3]; expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; testCase.verify_results(out, expected_fitpow, ... testCase.default_fields, ... @@ -99,16 +99,16 @@ function test_set_background_strategy_to_indep(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "planar"); % set some background parameters (correpsonding to planar bg - % with slope_en=slope_q=intercept =2 - out.set_bg_parameters(1:3, [2,2,2]); % en_slope = 2 + % with slope_en=1, slope_q=2, intercept = 3 + out.set_bg_parameters(1:3, [2,1,3]); % en_slope = 2 out.set_background_strategy("independent"); expected_fitpow = testCase.default_fitpow; % add extra background param expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); - expected_fitpow.params(2:2:end-1) = 2; % en slope - expected_fitpow.params(3) = 10; % intercept = 10 for cut 1 - expected_fitpow.params(5) = 12; % intercept = 12 for cut 2 + expected_fitpow.params(2:2:end-1) = 1; % en slope + expected_fitpow.params(3) = 11; % intercept = 10 for cut 1 + expected_fitpow.params(5) = 13; % intercept = 12 for cut 2 expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; testCase.verify_results(out, expected_fitpow, ... testCase.default_fields, ... From 718275e0b3f07523067e4641d3196a17cc94edcc Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 2 Jun 2025 18:24:14 +0100 Subject: [PATCH 21/33] Add ability to set background parameters, bounds etc. using string --- swfiles/sw_fitpowder.m | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index a49cce55..c69c9229 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -207,19 +207,21 @@ function initialise_parameters_and_bounds(obj, model_params) function initialise_background_parameters_and_bounds(obj) % zero intialise background parameters + obj.bg_param_labels = repmat("E", obj.npoly_en+1, 1) + [obj.npoly_en:-1:0]'; if obj.background_strategy == "independent" nparams_bg_total = obj.ncuts * obj.nparams_bg; else nparams_bg_total = obj.nparams_bg; + % add modQ polynopmial label + obj.bg_param_labels = [obj.bg_param_labels; + repmat("Q", obj.npoly_modQ, 1) + [obj.npoly_modQ:-1:1]']; end bg_params = zeros(nparams_bg_total, 1); bg_bounds = [-inf, inf].*ones(nparams_bg_total, 1); % insert in middle of params and bound obj.params = [obj.params(1:obj.nparams_model); bg_params; obj.params(end)]; obj.bounds = [obj.bounds(1:obj.nparams_model,:); bg_bounds; obj.bounds(end, :)]; - % set labels - obj.bg_param_labels = [repmat("Q", obj.npoly_modQ, 1) + [obj.npoly_modQ:-1:1]'; - repmat("E", obj.npoly_en+1, 1) + [obj.npoly_en:-1:0]']; + end function set_bg_npoly_modQ(obj, npoly_modQ) @@ -303,9 +305,6 @@ function fix_model_parameters(obj, iparams) end function fix_bg_parameters(obj, iparams_bg, icuts) - if any(iparams_bg > obj.nparams_bg) - error('sw_fitpowder', 'Parameter indices supplied must be within number of background parameters'); - end if nargin < 3 icuts = 0; end @@ -324,9 +323,6 @@ function set_model_parameters(obj, iparams, values) end function set_bg_parameters(obj, iparams_bg, values, icuts) - if any(iparams_bg > obj.nparams_bg) - error('sw_fitpowder', 'Parameter indices supplied must be within number of background parameters'); - end if nargin < 4 icuts = 0; end @@ -350,9 +346,6 @@ function set_model_parameter_bounds(obj, iparams, lb, ub) end function set_bg_parameter_bounds(obj, iparams_bg, lb, ub, icuts) - if any(iparams_bg > obj.nparams_bg) - error('sw_fitpowder', 'Parameter indices supplied must be within number of background parameters'); - end if nargin < 5 icuts = 0; end @@ -874,7 +867,26 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) end end + function iparams_bg = get_index_of_background_parameter_from_string(obj, bg_param_labels) + nlabels = numel(bg_param_labels); + iparams_bg = zeros(1, nlabels); + for ilabel = 1:nlabels + iparam_bg = find(obj.bg_param_labels == bg_param_labels(ilabel)); + assert(~isempty(iparam_bg), ... + 'sw_fitpowder:invalidinput', ... + ['Background parameter ' char(bg_param_labels(ilabel)) ' not found.']); + iparams_bg(ilabel) = iparam_bg; + end + end + function iparams = get_index_of_background_parameters(obj, iparams_bg, icuts) + % get index if parameter string passed + if isstring(iparams_bg) + iparams_bg = obj.get_index_of_background_parameter_from_string(iparams_bg); + end + if any(iparams_bg > obj.nparams_bg) + error('sw_fitpowder:invalidinput', 'Parameter indices supplied must be within number of background parameters'); + end if any(icuts == 0) if obj.background_strategy == "independent" icuts = 1:obj.ncuts; % apply to all cuts @@ -888,6 +900,7 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) iparams = [iparams, iparams_bg + obj.nparams_model + (icut-1)*obj.nparams_bg]; end end + function plot_1d_or_2d(obj, ycalc, varargin) if obj.liveplot_interval == 0 figure("color","white"); From 4031e9f1c5ca0bfc944947de890d347392a8d54c Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 2 Jun 2025 18:25:02 +0100 Subject: [PATCH 22/33] Add tests for setting bg parameters etc. with string labels --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 46 ++++++++++++++++++- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 3bfcb76e..0fc7299b 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -107,8 +107,8 @@ function test_set_background_strategy_to_indep(testCase) expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); expected_fitpow.params(2:2:end-1) = 1; % en slope - expected_fitpow.params(3) = 11; % intercept = 10 for cut 1 - expected_fitpow.params(5) = 13; % intercept = 12 for cut 2 + expected_fitpow.params(3) = 11; % intercept = 4*2 + 3 for cut 1 + expected_fitpow.params(5) = 13; % intercept = 5*2 + 3 for cut 2 expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; testCase.verify_results(out, expected_fitpow, ... testCase.default_fields, ... @@ -599,7 +599,49 @@ function test_export_data_1d(testCase) testCase.verify_val(data, testCase.data_1d_cuts); end + function test_set_bg_param_with_name_planar_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.set_bg_parameters("Q1", -5); + out.set_bg_parameters("E1", 5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(2:4) = [5, 0, -5]; % E1, E0, Q1 + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_param_with_invalid_name(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + testCase.verifyError(... + @() out.set_bg_parameters("Q3", -5), ... + 'sw_fitpowder:invalidinput'); + end + + function test_set_bg_param_with_name_indep_bg_all_cuts(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + bg_pars = [-5, 5]; + out.set_bg_parameters("E1", bg_pars(1)); + out.set_bg_parameters("E0", bg_pars(2)); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = [expected_fitpow.params(1); + bg_pars(:); bg_pars(:); 1]; + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + testCase.verify_results(out, expected_fitpow); + end + function test_fix_bg_param_with_array_of_names(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + bg_pars = 1:3; + bg_labels = ["E1","E0","Q1"]; + out.set_bg_parameters(bg_labels, bg_pars); + out.fix_bg_parameters(bg_labels); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(2:4) = bg_pars; + expected_fitpow.bounds(2:4, :) = repmat(bg_pars(:), 1,2); + testCase.verify_results(out, expected_fitpow); + end end end From 4dfae058c2d026c21400231967b5728bcf52abfe Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 5 Jun 2025 15:44:43 +0100 Subject: [PATCH 23/33] Fix bug in background param label order --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 4 ++-- swfiles/sw_fitpowder.m | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 0fc7299b..3721608b 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -605,7 +605,7 @@ function test_set_bg_param_with_name_planar_bg(testCase) out.set_bg_parameters("Q1", -5); out.set_bg_parameters("E1", 5); expected_fitpow = testCase.default_fitpow; - expected_fitpow.params(2:4) = [5, 0, -5]; % E1, E0, Q1 + expected_fitpow.params(2:4) = [-5, 5,0]; % Q1, E1, E0 testCase.verify_results(out, expected_fitpow); end @@ -634,7 +634,7 @@ function test_fix_bg_param_with_array_of_names(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... testCase.fit_func, testCase.j1); bg_pars = 1:3; - bg_labels = ["E1","E0","Q1"]; + bg_labels = ["Q1","E1","E0"]; out.set_bg_parameters(bg_labels, bg_pars); out.fix_bg_parameters(bg_labels); expected_fitpow = testCase.default_fitpow; diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index c69c9229..629a539a 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -213,8 +213,8 @@ function initialise_background_parameters_and_bounds(obj) else nparams_bg_total = obj.nparams_bg; % add modQ polynopmial label - obj.bg_param_labels = [obj.bg_param_labels; - repmat("Q", obj.npoly_modQ, 1) + [obj.npoly_modQ:-1:1]']; + obj.bg_param_labels = [repmat("Q", obj.npoly_modQ, 1) + [obj.npoly_modQ:-1:1]'; + obj.bg_param_labels;]; end bg_params = zeros(nparams_bg_total, 1); bg_bounds = [-inf, inf].*ones(nparams_bg_total, 1); @@ -701,7 +701,7 @@ function estimate_constant_background(obj) bg = mean(obj.y(obj.ibg)); end % set constant background assuming last bg parameter - obj.set_bg_parameters(obj.nparams_bg, bg); % set constant background + obj.set_bg_parameters("E0", bg); % set constant background end function set_errors_of_bg_bins(obj, val) From a7b4b0eca57ebf2cebcfa147824ba362d72533fa Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 5 Jun 2025 15:46:34 +0100 Subject: [PATCH 24/33] Add method to nicely print/display parameters --- swfiles/sw_fitpowder.m | 54 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 629a539a..c1acee62 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -807,6 +807,60 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) end end + function disp_params(obj,pfit, perr) + col_headers = ["Index", "Initial"]; + pars = obj.params(:); + par_fmt = "%8g\t"; + fmt_str = "%8.0f\t" + par_fmt; + if nargin > 1 + assert(numel(pfit) == numel(obj.params), ... + 'sw_fitpowder:invalidinput', ... + 'Fit params vector has invlaid length.'); + col_headers = [col_headers, "Best Fit"]; + pars = [pars, pfit(:)]; + fmt_str = fmt_str + par_fmt; + end + if nargin == 3 + assert(numel(perr) == numel(obj.params), ... + 'sw_fitpowder:invalidinput', ... + 'Fit error vector has invlaid length.'); + col_headers = [col_headers, "Error"]; + pars = [pars, perr(:)]; + fmt_str = fmt_str + par_fmt; + end + fmt_str = fmt_str + "\n"; + % print model params + fprintf("---\nMODEL\n---\n") + fprintf([repmat('%s\t\t', 1,numel(col_headers)), '\n'], col_headers); + imodel = [1:obj.nparams_model]'; + fprintf(fmt_str, [imodel, pars(imodel,:)]'); + % print scale + scales = sprintf(par_fmt, pars(end,:)); + fprintf("%8s\t%s\n", "SCALE", scales); + % print background + fprintf("---\nBACKGROUND\n---\n") + bg_pars = pars(obj.nparams_model+1:size(pars, 1)-1, :); + labels = obj.bg_param_labels; + col_headers = col_headers(2:end); % excl. Index + bg_par_index = [1:obj.nparams_bg]'; + if obj.ndim == 1 && obj.background_strategy == "independent" + % add icut column + col_headers = ["Cut Index", col_headers]; + fmt_str = "%.0f\t" + fmt_str; + cut_index = repelem([1:obj.ncuts]', obj.nparams_bg,1); + bg_pars = [cut_index , bg_pars]; + bg_par_index = repmat(bg_par_index, obj.ncuts, 1); + labels = repmat(labels, obj.ncuts, 1); + end + bg_pars = [bg_par_index, bg_pars]; + col_headers = ["Label", "Index", col_headers]; + fprintf([repmat('%s\t', 1,numel(col_headers)), '\n'], col_headers); + for irow = 1:size(bg_pars, 1) + row = sprintf(fmt_str, bg_pars(irow,:)); + fprintf("%5s\t%s", labels(irow), row); + end + end + end % private From 016b6046320244de85c7633574b225d8ff1028e7 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 5 Jun 2025 15:47:15 +0100 Subject: [PATCH 25/33] Update tutorial with parameter display --- tutorials/publish/tutorial39.m | 136 +++++++++++++++++++++++---------- 1 file changed, 97 insertions(+), 39 deletions(-) diff --git a/tutorials/publish/tutorial39.m b/tutorials/publish/tutorial39.m index e582a8d4..c18dd801 100644 --- a/tutorials/publish/tutorial39.m +++ b/tutorials/publish/tutorial39.m @@ -91,7 +91,7 @@ % region - it is worth trying this first fitpow.estimate_constant_background(); -fitpow.estimate_scale_factor() +fitpow.estimate_scale_factor(); % view the background region on 2D colorfill plot of data figure("color","white") @@ -107,6 +107,25 @@ % If the automatic background estimation doesn't look reasonable then it is % possible to define regions of |Q| and energy transfer that are background % and use these for fitting. +% +% General polynomial backgrounds in both |Q| and energy transfer are +% supported (governed by the attributes npoly_modQ and npoly_en). +% The default background is planar of the form +% +% Q1*|Q| + E1*en + E0 +% +% The background parameters can be referenced by index or label +% (e.g. "E0" for constant background). +% The labels are of the form "QN" or "EN" (for |Q| and energy +% polynomial coeficients repectively) where N represents the order of the +% polynomial. The background parameters are indexed in the following order +% +% [QN, QN-1, ..., Q1, EM, EM-1, ..., E1, E0] +% +% where N = npoly_modQ and M=npoly_en. Note there is no Q0 parameter, this +% would correspond to a constant which would be perfectly correlated with +% E0. + % set background regions %----------------------- @@ -125,10 +144,10 @@ %----------------------------------------------------- % fit only a constant/flat background -% default background is planar with parameters: -% [slope_energy, slope_modQ, intercept] -fitpow.fix_bg_parameters(1:2); % fix slopes to initial value (0) -fitpow.set_bg_parameter_bounds(3, 0.0, []) % Set lb of constant bg = 0 +fitpow.fix_bg_parameters(["Q1", "E1"]); % fix slopes to initial value (0) +fitpow.set_bg_parameter_bounds("E0", 0.0, []) % Set lb of constant bg = 0 + +% general polynomial backgrounds are supported the parameters are ordered [pfit, ~, stat] = fitpow.fit_background(); fitpow.estimate_scale_factor(); @@ -136,11 +155,29 @@ % view 1D cuts of data and the initial guess fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, fitpow.params); +% print inital parameters +fitpow.disp_params(); +% --- +% MODEL +% --- +% Index Initial +% 1 0.85 +% SCALE 1316.58 +% --- +% BACKGROUND +% --- +% Label Index Initial +% Q1 1 0 +% E1 2 0 +% E0 3 0.976825 %% Cropping data % Data can be cropped in |Q| and energy transfer - for example it is % advisible not to include the low-energy transfer region close to the % elastic line (relative to the resolution). +% +% You may want to re-optimise the background and scale factor after +% cropping % crop data fitpow.crop_energy_range(1.5*eres(0), inf); % typically want to avoid elastic line @@ -175,10 +212,10 @@ % set cost function and minimser fitpow.cost_function = "Rsq"; % "Rsq" or "chisq" -fitpow.optimizer = @ndbase.simplex; % or e.g. ndbase.simplex +fitpow.optimizer = @ndbase.simplex; % or e.g. ndbase.lm4 % Fix the constant background at the result of the previous fit -fitpow.fix_bg_parameters(3); % fix fitted constant bg +fitpow.fix_bg_parameters("E0"); % fix fitted constant bg % set bounds on first (and only in this case) model parameter fitpow.set_model_parameter_bounds(1, 0, []) % Force J_1 to be AFM @@ -191,10 +228,26 @@ % make 1D plots fitpow.plot_1d_cuts_of_2d_data(qcens-dq, qcens+dq, pfit) +% print parameters +fitpow.disp_params(pfit); +% --- +% MODEL +% --- +% Index Initial Best Fit +% 1 0.85 1.00185 +% SCALE 1316.58 1024.67 +% --- +% BACKGROUND +% --- +% Label Index Initial Best Fit +% Q1 1 0 0 +% E1 2 0 0 +% E0 3 0.976825 0.976825 + %% Fit 1D cuts % As discussed the sw_fitpowder class can be initialised with a cell array % of xye structs correpsonding to 1D cuts at constant |Q|. -% Alternatively, having loaded in 2D data the sw_fitpowder class has a +% Alternatively, having loaded in 2D data, the sw_fitpowder class has a % convenient method to take 1D cuts and replace the data. % % If 2D data are replaced with 1D cuts then the class evaluates the @@ -211,10 +264,8 @@ % change background to quadratic in |Q| fitpow.set_bg_npoly_modQ(2); % zeros previously optimised background -% Fix linear terms of background -% bg params ordered as en^n,...en^1, q^n,..q^1, const -fitpow.fix_bg_parameters(3); % fix term ~|Q| -fitpow.fix_bg_parameters(1); % fix term ~energy +% Fix linear terms of background to be 0 +fitpow.fix_bg_parameters(["Q1", "E1"]); % fit background and scale (fix model params) to all bins (not just bg % region -0 only advisable if near good solution) [pfit, ~, stat] = fitpow.fit_background_and_scale(); @@ -232,18 +283,22 @@ fitpow.plot_result(pfit); -% display result -par_names = ["J1", "Bg_En1", "Bg_Q2", "Bg_Q1", "Bg_Q0", "Scale"]; -["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] - -% "Name" "Value" "Error" -% "J1" "1.0027473" "0.00360547" -% "Bg_En1" "0" "0" -% "Bg_Q2" "0.02066691" "0.095079" -% "Bg_Q1" "0" "0" -% "Bg_Q0" "-0.036883773" "0.37194" -% "Scale" "1082.117" "19.405" - +% print parameters +fitpow.disp_params(pfit, stat.sigP); +% --- +% MODEL +% --- +% Index Initial Best Fit Error +% 1 1.00185 1.00091 0.00357452 +% SCALE 1106.42 1105.33 19.8056 +% --- +% BACKGROUND +% --- +% Label Index Initial Best Fit Error +% Q2 1 -0.618568 -0.613758 0.294827 +% Q1 2 0 0 0 +% E1 3 0 0 0 +% E0 4 0.444711 0.455551 0.497835 %% Change background strategy % The background for 1D cuts can be handled in two ways: @@ -265,18 +320,21 @@ fitpow.plot_result(pfit) -% display result -bg_names = repmat(["Bg_En1_cut", "Bg_En0_cut"]', fitpow.ncuts, 1); -bg_names = bg_names + repelem([1:fitpow.ncuts]', fitpow.nparams_bg,1); -par_names = ["J1", bg_names', "Scale"]; -["Name", "Value", "Error"; par_names(:), pfit(:), stat.sigP(:)] - -% "Name" "Value" "Error" -% "J1" "1.0059739" "0.00403765" -% "Bg_En1_cut1" "-0.85631989" "0.376039" -% "Bg_En0_cut1" "1.7546252" "0.892913" -% "Bg_En1_cut2" "0.52590868" "0.433608" -% "Bg_En0_cut2" "-2.3563309" "1.49539" -% "Bg_En1_cut3" "-0.41373845" "0.350704" -% "Bg_En0_cut3" "0.52287368" "0.936986" -% "Scale" "1131.9098" "30.4495" +% print parameters +fitpow.disp_params(pfit, stat.sigP); +% --- +% MODEL +% --- +% Index Initial Best Fit Error +% 1 1.00091 1.00181 0.00361181 +% SCALE 1105.33 1213.84 29.0031 +% --- +% BACKGROUND +% --- +% Label Index Cut Index Initial Best Fit Error +% E1 1 1 0 -1.03769 0.359709 +% E0 2 1 0.0977161 1.79735 0.857853 +% E1 1 2 0 1.85999 0.41829 +% E0 2 2 -0.926966 -7.27229 1.43701 +% E1 1 3 0 0.23947 0.337007 +% E0 2 3 -2.61892 -2.66401 0.904453 From f0c6c2a134589d5886cbf74a7adc03bba397b141 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 11:31:44 +0100 Subject: [PATCH 26/33] Stop scaling covariance by variance of residuals This is only valid in some cases --- swfiles/+ndbase/lm4.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/swfiles/+ndbase/lm4.m b/swfiles/+ndbase/lm4.m index c31cf2de..d83f4640 100644 --- a/swfiles/+ndbase/lm4.m +++ b/swfiles/+ndbase/lm4.m @@ -227,7 +227,7 @@ fVal = cost_val / ndof; if exit_flag > 0 % converged on solution - calculate errors - cov = pinv(hess) * 2.0 * fVal; + cov = pinv(hess) * 2.0; perr(cost_func_wrap.ifree) = sqrt(diag(cov)); else message = "Failed to converge in MaxIter"; From 5263316938aaf1e7db5a8b63ea557df23944ddc5 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 11:33:05 +0100 Subject: [PATCH 27/33] Fix bug in getting qcens of cuts when using binning in 2D data --- swfiles/sw_fitpowder.m | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index c1acee62..de379d3a 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -876,7 +876,10 @@ function disp_params(obj,pfit, perr) function modQ_cens = get_modQ_cens_of_cuts(obj) - modQ_cens = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); + modQ_cens = zeros(1, obj.ncuts); + for icut = 1:obj.ncuts + modQ_cens(icut) = mean(obj.modQ_cens(obj.modQ_icuts == icut)); + end end function cut = cut_2d_data(obj, qmin, qmax, ycalc) From 6e687a1424348fb1b674d06513fcfd1a87cf9fe7 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 11:34:00 +0100 Subject: [PATCH 28/33] Improve plot formatting --- swfiles/sw_fitpowder.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index de379d3a..33e3396b 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -754,7 +754,7 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) ax = gca; box on; hold on; h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, y); - h.AlphaData = double(obj.y > 0); % make empty bins transparent + h.AlphaData = double(obj.e > 0); % make empty bins transparent cbar = colorbar(ax); cbar.Label.String = "Intensity"; xlabel(ax, "$\left|Q\right| (\AA^{-1})$", 'interpreter','latex'); @@ -791,11 +791,11 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) hold on; box on; if ~isempty(ycalc) cut = obj.cut_2d_data(qmins(icut), qmaxs(icut), ycalc); - plot(ax, cut.x, cut.ycalc, '-r'); + plot(ax, cut.x, cut.ycalc, '-r', 'DisplayName', 'Fit'); else cut = obj.cut_2d_data(qmins(icut), qmaxs(icut)); end - plot(ax, cut.x, cut.y, 'ok'); + plot(ax, cut.x, cut.y, 'ok', 'DisplayName', 'Data'); % calc xlims ifinite = isfinite(cut.y); istart = find(ifinite, 1, 'first'); From 78177cfabd93e7bc39cb406b692d5ab51f5e3677 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 14:00:34 +0100 Subject: [PATCH 29/33] Support char and cell of char in bg param setter --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 25 ++++++++++++++++--- swfiles/sw_fitpowder.m | 5 +++- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 3721608b..c4e0a001 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -81,7 +81,7 @@ function test_init_data_1d_indep_bg(testCase) function test_set_background_strategy_to_planar(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent"); - % set some background parameters (correpsonding to planar bg + % set some background parameters for 1D cuts equivalent to a planar bg % with slope_en=1, slope_q=2, intercept = 3 out.set_bg_parameters(1, 1); % en_slope = 1 out.set_bg_parameters(2, 11, 1); % intercept = 4*2 + 3 for cut 1 @@ -182,6 +182,25 @@ function test_set_bg_parameters_planar_bg(testCase) testCase.verify_results(out, expected_fitpow); end + function test_set_bg_parameters_with_char(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.set_bg_parameters('E0', -5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(end-1) = -5; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_parameters_with_cell_of_char(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + bg_pars = [-5, 5]; + out.set_bg_parameters({'Q1', 'E1'}, bg_pars); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(2:3) = bg_pars; + testCase.verify_results(out, expected_fitpow); + end + function test_set_bg_parameters_indep_bg_all_cuts(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent", 1); @@ -621,8 +640,7 @@ function test_set_bg_param_with_name_indep_bg_all_cuts(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1, "independent", 1); bg_pars = [-5, 5]; - out.set_bg_parameters("E1", bg_pars(1)); - out.set_bg_parameters("E0", bg_pars(2)); + out.set_bg_parameters(["E1", "E0"], bg_pars); expected_fitpow = testCase.default_fitpow; expected_fitpow.params = [expected_fitpow.params(1); bg_pars(:); bg_pars(:); 1]; @@ -642,6 +660,7 @@ function test_fix_bg_param_with_array_of_names(testCase) expected_fitpow.bounds(2:4, :) = repmat(bg_pars(:), 1,2); testCase.verify_results(out, expected_fitpow); end + end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 33e3396b..b0091969 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -254,7 +254,7 @@ function set_bg_npoly_en(obj, npoly_en) function set_background_strategy(obj, strategy) if strategy == "planar" - obj.fbg = @(en, modQ, p) obj.fbg_planar(en, modQ, p, obj.npoly_en); + obj.fbg = @(en, modQ, p) obj.fbg_planar(en, modQ, p, obj.npoly_modQ); elseif strategy == "independent" if obj.ndim == 2 error('sw_fitpowder:invalidinput', ... @@ -326,6 +326,9 @@ function set_bg_parameters(obj, iparams_bg, values, icuts) if nargin < 4 icuts = 0; end + if ischar(iparams_bg) || iscell(iparams_bg) % Handles character arrays + iparams_bg = string(iparams_bg); + end for ival = 1:numel(values) iparams = obj.get_index_of_background_parameters(iparams_bg(ival), icuts); obj.params(iparams) = values(ival); From 4734e78276b68232e2561d892ddeb84242c1469d Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 14:35:08 +0100 Subject: [PATCH 30/33] Save file when export data --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 11 +++++++++++ swfiles/sw_fitpowder.m | 5 ++++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index c4e0a001..c477425e 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -611,6 +611,17 @@ function test_export_data_2d(testCase) testCase.verify_val(data, testCase.data_2d); end + function test_export_data_2d_with_filename(testCase) + import matlab.unittest.fixtures.TemporaryFolderFixture + fixture = testCase.applyFixture(TemporaryFolderFixture); + tmp_file = fullfile(fixture.Folder,"data.mat"); + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + data = out.export_data(tmp_file); + testCase.verify_val(data, testCase.data_2d); + testCase.assertTrue(isfile(tmp_file)); + end + function test_export_data_1d(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... testCase.fit_func, testCase.j1); diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index b0091969..d484e8f3 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -362,7 +362,7 @@ function set_scale_bounds(obj, lb, ub) obj.set_bounds(size(obj.bounds,1), lb, ub); end - function data = export_data(obj) + function data = export_data(obj, filename) if obj.ndim == 2 data.y = obj.y'; data.e = obj.e'; @@ -379,6 +379,9 @@ function set_scale_bounds(obj, lb, ub) data(icut).qmax = qs(end) + 0.5*(qs(end)-qs(end-1)); end end + if nargin > 1 + save(filename, 'data'); + end end function add_data(obj, data) From 8e0376ccee02bc0bbb6b996fb7760c1e8284ddfd Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 15:19:01 +0100 Subject: [PATCH 31/33] Add check for integer npoly --- swfiles/sw_fitpowder.m | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index d484e8f3..e10231eb 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -229,25 +229,26 @@ function set_bg_npoly_modQ(obj, npoly_modQ) error('sw_fitpowder:invalidinput', ... 'set_bg_npoly_modQ only applies to planar background'); end - if npoly_modQ < 0 - error('sw_fitpowder:invalidinput', ... - 'npoly_modQ must be positive int'); + if abs(mod(npoly_modQ, 1)) > eps || npoly_modQ < 0 + error('sw_fitpowder:invalidinput', ... + 'npoly_modQ must be a positive integer'); end if obj.ndim==1 && npoly_modQ > obj.ncuts - 1 error('sw_fitpowder:invalidinput', ... 'Not enough 1D cuts to fit polynomial of order npoly_modQ.'); end - obj.npoly_modQ = round(npoly_modQ); + + obj.npoly_modQ = npoly_modQ; obj.nparams_bg = obj.get_nparams_in_background_func(); obj.initialise_background_parameters_and_bounds(); end function set_bg_npoly_en(obj, npoly_en) - if npoly_en < 0 - error('sw_fitpowder:invalidinput', ... - 'npoly_en must be positive int'); + if abs(mod(npoly_en, 1)) > eps || npoly_en < 0 + error('sw_fitpowder:invalidinput', ... + 'npoly_en must be a positive integer'); end - obj.npoly_en = round(npoly_en); + obj.npoly_en = npoly_en; obj.nparams_bg = obj.get_nparams_in_background_func(); obj.initialise_background_parameters_and_bounds(); end From 92bbcd32f2e7f22d6443278c093a0ccf6e70ded4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 3 Sep 2025 15:46:30 +0100 Subject: [PATCH 32/33] Remove comment form tutorial docs --- swfiles/sw_fitpowder.m | 2 +- tutorials/publish/tutorial39.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index e10231eb..eb4263a1 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -420,11 +420,11 @@ function add_data(obj, data) 'Input cell does not have correct fields'); obj.y = [obj.y cut.y(:)]; obj.e = [obj.e cut.e(:)]; - dQ = (cut.qmax - cut.qmin)/obj.nQ; if isfield(cut, 'qs') obj.modQ_cens = [obj.modQ_cens, cut.qs(:)']; obj.modQ_icuts = [obj.modQ_icuts, icut*ones(1, numel(cut.qs))]; else + dQ = (cut.qmax - cut.qmin)/obj.nQ; obj.modQ_cens = [obj.modQ_cens, ((cut.qmin+dQ/2):dQ:cut.qmax)]; obj.modQ_icuts = [obj.modQ_icuts, icut*ones(1, obj.nQ)]; end diff --git a/tutorials/publish/tutorial39.m b/tutorials/publish/tutorial39.m index c18dd801..f909346a 100644 --- a/tutorials/publish/tutorial39.m +++ b/tutorials/publish/tutorial39.m @@ -78,7 +78,7 @@ fitpow.powspec_args.dE = eres; % emergy resolution fitpow.powspec_args.fastmode = true; fitpow.powspec_args.neutron_output = true; -fitpow.powspec_args.nRand = 1e3; % low for speed (typically want > 1e3) +fitpow.powspec_args.nRand = 1e3; fitpow.powspec_args.hermit = true; % set parameters passsed to sw_instrument fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei); From 6ae56d9b4fdfc740348bafd860fb0dc9d5b6cd36 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 9 Sep 2025 14:45:03 +0100 Subject: [PATCH 33/33] Support specifying 1D cuts with array of q-values and test --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 8 ++++++++ swfiles/sw_fitpowder.m | 4 +++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index c477425e..4a8766de 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -672,6 +672,14 @@ function test_fix_bg_param_with_array_of_names(testCase) testCase.verify_results(out, expected_fitpow); end + function test_add_1Dcuts_withs_qs_specified(testCase) + cuts = arrayfun(@(qs) struct('x', 1:3, 'y', 1:3, 'e', 1:3, 'qs', qs), ... + 4:5); + out = sw_fitpowder(testCase.swobj, cuts, ... + testCase.fit_func, testCase.j1); + testCase.verify_results(out, testCase.default_fitpow); + end + end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index eb4263a1..ddc3ab5e 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -415,7 +415,9 @@ function add_data(obj, data) obj.ncuts = numel(data); for icut = 1:obj.ncuts cut = data(icut); - assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... + has_xye = all(isfield(cut, {'x', 'y', 'e'})); + has_qfields = all(isfield(cut, {'qmin', 'qmax'})) || isfield(cut, 'qs'); + assert(has_xye && has_qfields, ... 'sw_fitpowder:invalidinput', ... 'Input cell does not have correct fields'); obj.y = [obj.y cut.y(:)];