Skip to content

More usability improvements sw fitpowder #210

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
RichardWaiteSTFC merged 33 commits into from
Sep 10, 2025
Merged

More usability improvements sw fitpowder #210

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
a7e78cd
Fix bug in estimate_scale_factor - background already added to ycalc
RichardWaiteSTFC Jan 8, 2025
f7dd48e
Add new function to fit bg and scale parameters to all data
RichardWaiteSTFC Jan 8, 2025
127156e
Support any ndbase optimizer in fit_background
RichardWaiteSTFC Jan 8, 2025
d04fe5c
Add method to plot background region on figures
RichardWaiteSTFC Jan 9, 2025
9aa01e4
Add function (and tests) to set background region
RichardWaiteSTFC Jan 9, 2025
4c25403
Set sensible background when changing strategy
RichardWaiteSTFC Jan 9, 2025
2869cff
Handle edges/nans properly in avg. of error when taking 1D cuts
RichardWaiteSTFC Jan 9, 2025
cc82c7f
Add method to clear background region
RichardWaiteSTFC Feb 7, 2025
91e673b
Support general polynomial background in Q and en
RichardWaiteSTFC Mar 17, 2025
4464c7f
Add unit tests for setting polynomial order
RichardWaiteSTFC Mar 17, 2025
0ad70f0
Use 2D data Q bins when replace 1D cuts
RichardWaiteSTFC Apr 28, 2025
df70190
Set SigP to 0 for fixed parameters in lm4
RichardWaiteSTFC Apr 29, 2025
40eca62
Fix bug for 1D cuts now using 2D Q bins
RichardWaiteSTFC Apr 29, 2025
7a6a70e
Update tutorial to include setting bg strategy and setting bg region
RichardWaiteSTFC Apr 29, 2025
7b470aa
Add method to export data in struct form
RichardWaiteSTFC May 16, 2025
f8e4567
Fix rebin_powspec_to_1D_cuts bug causing test failure
RichardWaiteSTFC May 19, 2025
17e2b1c
Add unit tests for export_data method
RichardWaiteSTFC May 19, 2025
484018d
Update tutorial with changing polynomial order
RichardWaiteSTFC May 19, 2025
e934108
Swap order of bg parameters (Q first) and add bg param labels
RichardWaiteSTFC Jun 2, 2025
12657f6
Update changing bg strategy tests with non-identical bg parameters
RichardWaiteSTFC Jun 2, 2025
718275e
Add ability to set background parameters, bounds etc. using string
RichardWaiteSTFC Jun 2, 2025
4031e9f
Add tests for setting bg parameters etc. with string labels
RichardWaiteSTFC Jun 2, 2025
4dfae05
Fix bug in background param label order
RichardWaiteSTFC Jun 5, 2025
a7b4b0e
Add method to nicely print/display parameters
RichardWaiteSTFC Jun 5, 2025
016b604
Update tutorial with parameter display
RichardWaiteSTFC Jun 5, 2025
f0c6c2a
Stop scaling covariance by variance of residuals
RichardWaiteSTFC Sep 3, 2025
5263316
Fix bug in getting qcens of cuts when using binning in 2D data
RichardWaiteSTFC Sep 3, 2025
6e687a1
Improve plot formatting
RichardWaiteSTFC Sep 3, 2025
78177cf
Support char and cell of char in bg param setter
RichardWaiteSTFC Sep 3, 2025
4734e78
Save file when export data
RichardWaiteSTFC Sep 3, 2025
8e0376c
Add check for integer npoly
RichardWaiteSTFC Sep 3, 2025
92bbcd3
Remove comment form tutorial docs
RichardWaiteSTFC Sep 3, 2025
6ae56d9
Support specifying 1D cuts with array of q-values and test
RichardWaiteSTFC Sep 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
244 changes: 237 additions & 7 deletions +sw_tests/+unit_tests/unittest_sw_fitpowder.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Comment on lines +87 to +88
Copy link
Member

@mducle mducle Jun 25, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to the comment in line 85, the intercept index is 3 - but here you're changing parameter index 2 (which should slope_q, right?)

Copy link
Collaborator Author

@RichardWaiteSTFC RichardWaiteSTFC Sep 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, just checked this and the comment on line 85 refers to the equivalent planar background (which are checked on line 91). The parameters being set here are for the individual 1D cuts - i.e. parameter 2 corresponds to E0 - this is what disp_params gives

---
BACKGROUND
---
Label	Index	Cut Index	Initial	
   E1	1	       1	       1	
   E0	2	       1	      11	
   E1	1	       2	       1	
   E0	2	       2	      13

I will update the comment on line 85 to be clearer

Copy link
Member

@mducle mducle Sep 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@RichardWaiteSTFC bump...

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)
Expand All @@ -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

Expand All @@ -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],:);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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];
Expand All @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
6 changes: 3 additions & 3 deletions swfiles/+ndbase/lm4.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading
Loading