diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 925ce96d7..4a8766de4 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 @@ -76,13 +78,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 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 + 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,1,3]; 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=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) = 1; % en slope + 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, ... + 'abs_tol', 1e-10); end function test_replace_2D_data_with_1D_cuts(testCase) @@ -91,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 @@ -102,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],:); @@ -152,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); @@ -331,7 +380,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]; @@ -342,6 +391,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); @@ -406,11 +470,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 @@ -450,6 +515,171 @@ 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 + + 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 + + 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_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); + data = out.export_data(); + 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, 5,0]; % Q1, E1, E0 + 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", "E0"], bg_pars); + 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 = ["Q1","E1","E0"]; + 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 + + 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/+ndbase/lm4.m b/swfiles/+ndbase/lm4.m index ace810592..d83f46403 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)); + cov = pinv(hess) * 2.0; + perr(cost_func_wrap.ifree) = sqrt(diag(cov)); else message = "Failed to converge in MaxIter"; - perr = zeros(size(pOpt)); end end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index e98efd89c..ddc3ab5e9 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 @@ -161,15 +162,19 @@ 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_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 = [] model_params_cached = [] liveplot_counter = 0 ibg = [] bg_errors = [] + npoly_modQ = 1; + npoly_en = 1; + bg_param_labels = []; end properties (Constant) @@ -202,35 +207,94 @@ 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 = [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); % 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, :)]; + 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'); + 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 obj.background_strategy == "planar" - obj.fbg = obj.fbg_planar; - else + 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 = 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 abs(mod(npoly_en, 1)) > eps || npoly_en < 0 + error('sw_fitpowder:invalidinput', ... + 'npoly_en must be a positive integer'); + end + obj.npoly_en = 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 = @(en, modQ, p) obj.fbg_planar(en, modQ, p, obj.npoly_modQ); + 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 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 + % attempt good guess for background pars + modQs = obj.get_modQ_cens_of_cuts(); + if obj.background_strategy == "independent" + % 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_modQ); + for icut = 1:obj.ncuts + obj.set_bg_parameters(obj.npoly_en+1, intercepts(icut), icut); + end + elseif obj.background_strategy == "planar" + % 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(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(obj.npoly_modQ+1:obj.nparams_bg-1, bg_pars_en); + end + end end end @@ -242,9 +306,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 @@ -263,12 +324,12 @@ 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 + 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); @@ -289,9 +350,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 @@ -305,6 +363,28 @@ function set_scale_bounds(obj, lb, ub) obj.set_bounds(size(obj.bounds,1), lb, ub); end + function data = export_data(obj, filename) + 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 (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 + if nargin > 1 + save(filename, 'data'); + end + end + function add_data(obj, data) if ~isa(data, "struct") data = arrayfun(@obj.convert_horace_to_struct, data); @@ -333,14 +413,23 @@ function add_data(obj, data) obj.ndim = 1; obj.ebin_cens = data(1).x; obj.ncuts = numel(data); - for cut = data - assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... + for icut = 1:obj.ncuts + cut = data(icut); + 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(:)]; 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, 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 end end obj.powspec_args.Evect = obj.ebin_cens(:)'; % row vect @@ -356,9 +445,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 @@ -386,16 +472,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 @@ -411,8 +496,37 @@ 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) + % 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) @@ -504,6 +618,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 @@ -536,7 +663,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 @@ -554,18 +681,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) @@ -573,7 +695,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 @@ -584,7 +710,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) @@ -611,7 +737,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; @@ -623,7 +749,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 @@ -634,10 +760,10 @@ 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 = 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'); @@ -646,6 +772,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 @@ -660,22 +796,77 @@ 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', 'DisplayName', 'Fit'); + else + cut = obj.cut_2d_data(qmins(icut), qmaxs(icut)); end + plot(ax, cut.x, cut.y, 'ok', 'DisplayName', 'Data'); % 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 + + 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 @@ -692,31 +883,76 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) resid = resid(isfinite(resid) & e > obj.zero_abs_tol); end - function cut = cut_2d_data(obj, qmin, qmax) + + function modQ_cens = get_modQ_cens_of_cuts(obj) + 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) 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(ikeep); + cut.e = sqrt(sum(obj.e(:, ikeep).^2, 2))./sum(ifinite, 2); + if nargin == 4 && ~isempty(ycalc) + % also cut ycalc + cut.ycalc = mean(ycalc(:, ikeep), 2, 'omitnan'); + end 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')); + 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_1d, 2, 'omitnan')); + else + % different num Q points per cut, loop over cuts + 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'); + end + end end 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 + + 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 @@ -730,6 +966,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"); @@ -797,6 +1034,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) diff --git a/tutorials/publish/tutorial39.m b/tutorials/publish/tutorial39.m index b52e514b8..f909346a5 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,124 @@ 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; 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. +% +% 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 +%----------------------- + +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); + + +% fit background (to voxels in background region) +%----------------------------------------------------- + +% fit only a constant/flat background +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(); + +% 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 +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,21 +207,15 @@ % 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 +fitpow.optimizer = @ndbase.simplex; % or e.g. ndbase.lm4 -% 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(); -fitpow.fix_bg_parameters(3); % fix fitted constant bg +% Fix the constant background at the result of the previous fit +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 @@ -139,48 +226,115 @@ 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) +% 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. % -% 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); + +% change background to quadratic in |Q| +fitpow.set_bg_npoly_modQ(2); % zeros previously optimised background +% 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(); +% 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 = 1e4; +fit_kwargs = {'resid_handle', true, 'diff_step', 1e-3}; + +[pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); + +fitpow.plot_result(pfit); + +% 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: % % - 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 +% 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. -[pfit, cost_val, stat] = fitpow.fit(); - -fitpow.plot_result(pfit) +fitpow.params = pfit(:); +fitpow.set_background_strategy("independent"); +[pfit, cost_val, stat] = fitpow.fit(fit_kwargs{:}); +fitpow.plot_result(pfit) +% 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