Skip to content

Commit dc06a4b

Browse files
RichardWaiteSTFCRichardWaiteSTFC
authored
Merge pull request #210 from SpinW/more_usability_improvements_sw_fitpowder
Fix bug in scale factor estimation Add new function to fit bg and scale to all data (fixing model) Support any ndbase optimizer in fit_background Add methods to set and clear background region Add method to plot background region Initialise background parameters with sensible values when switching background strategy Handle edges/nans properly in avg. of error when taking 1D cuts Support general polynomial background in Q and en Use 2D data Q bins when replace with 1D cuts (previously just used nQ bins between qmin and qmax) Update tutorial Add method to export data into structs (or array of structs for 1D) - for fitbenchmarking Add ability to set background parameters, bounds etc. by string label Add method to print parameters nicely
2 parents ba7a369 + 6ae56d9 commit dc06a4b

File tree

4 files changed

+777
-133
lines changed

4 files changed

+777
-133
lines changed

‎+sw_tests/+unit_tests/unittest_sw_fitpowder.m‎

Lines changed: 237 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,9 @@ function setup_spinw_obj_and_expected_result(testCase)
3333
-Inf Inf;
3434
-Inf Inf;
3535
0 Inf], ...
36-
'ibg', []);
36+
'ibg', [], ...
37+
'npoly_modQ', 1, ...
38+
'npoly_en', 1);
3739
testCase.default_fields = fieldnames(testCase.default_fitpow);
3840
end
3941
end
@@ -76,13 +78,41 @@ function test_init_data_1d_indep_bg(testCase)
7678
testCase.verify_results(out, expected_fitpow);
7779
end
7880

79-
function test_set_background_strategy(testCase)
81+
function test_set_background_strategy_to_planar(testCase)
8082
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
8183
testCase.fit_func, testCase.j1, "independent");
84+
% set some background parameters for 1D cuts equivalent to a planar bg
85+
% with slope_en=1, slope_q=2, intercept = 3
86+
out.set_bg_parameters(1, 1); % en_slope = 1
87+
out.set_bg_parameters(2, 11, 1); % intercept = 4*2 + 3 for cut 1
88+
out.set_bg_parameters(2, 13, 2); % intercept = 5*2 + 3 for cut 2
8289
out.set_background_strategy("planar");
8390
expected_fitpow = testCase.default_fitpow;
91+
expected_fitpow.params(2:end-1) = [2,1,3];
8492
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
85-
testCase.verify_results(out, expected_fitpow);
93+
testCase.verify_results(out, expected_fitpow, ...
94+
testCase.default_fields, ...
95+
'abs_tol', 1e-10);
96+
end
97+
98+
function test_set_background_strategy_to_indep(testCase)
99+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
100+
testCase.fit_func, testCase.j1, "planar");
101+
% set some background parameters (correpsonding to planar bg
102+
% with slope_en=1, slope_q=2, intercept = 3
103+
out.set_bg_parameters(1:3, [2,1,3]); % en_slope = 2
104+
out.set_background_strategy("independent");
105+
expected_fitpow = testCase.default_fitpow;
106+
% add extra background param
107+
expected_fitpow.params = expected_fitpow.params([1:2,2:end],:);
108+
expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:);
109+
expected_fitpow.params(2:2:end-1) = 1; % en slope
110+
expected_fitpow.params(3) = 11; % intercept = 4*2 + 3 for cut 1
111+
expected_fitpow.params(5) = 13; % intercept = 5*2 + 3 for cut 2
112+
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
113+
testCase.verify_results(out, expected_fitpow, ...
114+
testCase.default_fields, ...
115+
'abs_tol', 1e-10);
86116
end
87117

88118
function test_replace_2D_data_with_1D_cuts(testCase)
@@ -91,7 +121,7 @@ function test_replace_2D_data_with_1D_cuts(testCase)
91121
qcens = [4, 5];
92122
out.replace_2D_data_with_1D_cuts(qcens-0.5, qcens+0.5)
93123
expected_fitpow = testCase.default_fitpow;
94-
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
124+
expected_fitpow.modQ_cens = qcens;
95125
testCase.verify_results(out, expected_fitpow);
96126
end
97127

@@ -102,7 +132,7 @@ function test_replace_2D_data_with_1D_cuts_specify_bg(testCase)
102132
out.replace_2D_data_with_1D_cuts(qcens-0.5, qcens+0.5,...
103133
"independent")
104134
expected_fitpow = testCase.default_fitpow;
105-
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
135+
expected_fitpow.modQ_cens = qcens; % not nQ values as using bins in 2D data
106136
% add extra background param
107137
expected_fitpow.params = expected_fitpow.params([1:2,2:end],:);
108138
expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:);
@@ -152,6 +182,25 @@ function test_set_bg_parameters_planar_bg(testCase)
152182
testCase.verify_results(out, expected_fitpow);
153183
end
154184

185+
function test_set_bg_parameters_with_char(testCase)
186+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
187+
testCase.fit_func, testCase.j1);
188+
out.set_bg_parameters('E0', -5);
189+
expected_fitpow = testCase.default_fitpow;
190+
expected_fitpow.params(end-1) = -5;
191+
testCase.verify_results(out, expected_fitpow);
192+
end
193+
194+
function test_set_bg_parameters_with_cell_of_char(testCase)
195+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
196+
testCase.fit_func, testCase.j1);
197+
bg_pars = [-5, 5];
198+
out.set_bg_parameters({'Q1', 'E1'}, bg_pars);
199+
expected_fitpow = testCase.default_fitpow;
200+
expected_fitpow.params(2:3) = bg_pars;
201+
testCase.verify_results(out, expected_fitpow);
202+
end
203+
155204
function test_set_bg_parameters_indep_bg_all_cuts(testCase)
156205
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
157206
testCase.fit_func, testCase.j1, "independent", 1);
@@ -331,7 +380,7 @@ function test_fit_background(testCase, fit_params)
331380
out.y(1) = 10; % higher so other bins are background
332381
out.fix_bg_parameters(1:2); % fix slopes of background to 0
333382
out.set_bg_parameters(3, 1.5); % initial guess
334-
out.fit_background(fit_params{:})
383+
out.fit_background(fit_params{:});
335384
expected_fitpow = testCase.default_fitpow;
336385
expected_fitpow.y(1) = 10;
337386
expected_fitpow.ibg = [3;6;2;5;4];
@@ -342,6 +391,21 @@ function test_fit_background(testCase, fit_params)
342391
'abs_tol', 1e-4);
343392
end
344393

394+
function test_fit_background_and_scale(testCase)
395+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
396+
testCase.fit_func, testCase.j1);
397+
out.fix_bg_parameters(1:2); % fix slopes of background to 0
398+
out.set_bg_parameters(3, 1.5); % initial guess
399+
out.fit_background_and_scale();
400+
expected_fitpow = testCase.default_fitpow;
401+
expected_fitpow.params(end-1) = 0.0029;
402+
expected_fitpow.params(end) = 15.47;
403+
expected_fitpow.bounds(2:3,:) = 0; % fixed bg slopes
404+
testCase.verify_results(out, expected_fitpow, ...
405+
testCase.default_fields, ...
406+
'abs_tol', 1e-3);
407+
end
408+
345409
function test_calc_cost_func_of_background_indep(testCase)
346410
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
347411
testCase.fit_func, testCase.j1, "independent", 2);
@@ -406,11 +470,12 @@ function test_calc_uncertainty(testCase)
406470
function test_estimate_scale_factor(testCase)
407471
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
408472
testCase.fit_func, testCase.j1);
473+
out.set_bg_parameters(3, 0.05); % set constant bg
409474
out.powspec_args.dE = 0.1; % constant energy resolution
410475
out.powspec_args.hermit = true;
411476
out.estimate_scale_factor()
412477
expected_fitpow = testCase.default_fitpow;
413-
expected_fitpow.params(end) = 17.6;
478+
expected_fitpow.params(end) = 17.36;
414479
testCase.verify_results(out, expected_fitpow, ...
415480
testCase.default_fields, 'abs_tol', 1e-1);
416481
end
@@ -450,6 +515,171 @@ function test_add_1Dcuts_after_2D_data(testCase)
450515
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d; % integrtate over nQ pts
451516
testCase.verify_results(out, expected_fitpow);
452517
end
518+
519+
function test_set_bg_region_data_2d(testCase)
520+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
521+
testCase.fit_func, testCase.j1);
522+
out.set_bg_region(0,1.5); % for all Q
523+
out.set_bg_region(2.5,3.5,4.5,inf); % for highest Q
524+
expected_fitpow = testCase.default_fitpow;
525+
expected_fitpow.ibg = [1;4;6];
526+
testCase.verify_results(out, expected_fitpow);
527+
end
528+
529+
function test_set_bg_region_data_1d(testCase)
530+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
531+
testCase.fit_func, testCase.j1);
532+
expected_fitpow = testCase.default_fitpow;
533+
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
534+
out.set_bg_region(0,1.5); % for all cuts
535+
out.set_bg_region(2.5,3.5,2); % for last cut
536+
expected_fitpow.ibg = [1;4;6];
537+
testCase.verify_results(out, expected_fitpow);
538+
end
539+
540+
function test_set_npoly_modQ(testCase)
541+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
542+
testCase.fit_func, testCase.j1);
543+
out.set_bg_npoly_modQ(2);
544+
expected_fitpow = testCase.default_fitpow;
545+
expected_fitpow.npoly_modQ = 2;
546+
expected_fitpow.params = [expected_fitpow.params(1:end-1);
547+
0;
548+
expected_fitpow.params(end)];
549+
expected_fitpow.bounds = [expected_fitpow.bounds(1:end-1,:);
550+
-Inf, Inf;
551+
expected_fitpow.bounds(end,:)];
552+
testCase.verify_results(out, expected_fitpow);
553+
end
554+
555+
function test_set_npoly_en(testCase)
556+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
557+
testCase.fit_func, testCase.j1);
558+
out.set_bg_npoly_en(2);
559+
expected_fitpow = testCase.default_fitpow;
560+
expected_fitpow.npoly_en = 2;
561+
% add extra row in params and bounds
562+
expected_fitpow.params = expected_fitpow.params([1:2,2:end]);
563+
expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:);
564+
testCase.verify_results(out, expected_fitpow);
565+
end
566+
567+
function test_set_npoly_modQ_1d_data_indep_bg(testCase)
568+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
569+
testCase.fit_func, testCase.j1, "independent");
570+
testCase.verifyError(...
571+
@() out.set_bg_npoly_modQ(2), ...
572+
'sw_fitpowder:invalidinput');
573+
end
574+
575+
function test_set_npoly_en_1d_data_indep_bg(testCase)
576+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
577+
testCase.fit_func, testCase.j1, "independent");
578+
out.set_bg_npoly_en(2)
579+
expected_fitpow = testCase.default_fitpow;
580+
expected_fitpow.npoly_en = 2;
581+
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
582+
% add 3 extra rows:
583+
% 3 params x 2 cuts vs 3 planar parmas order=1
584+
expected_fitpow.params = expected_fitpow.params([1:4,2:end]);
585+
expected_fitpow.bounds = expected_fitpow.bounds([1:4,2:end],:);
586+
testCase.verify_results(out, expected_fitpow);
587+
end
588+
589+
function test_set_npoly_modQ_1d_data_planar_bg(testCase)
590+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
591+
testCase.fit_func, testCase.j1, "planar");
592+
% try adding order larger than numebr cuts
593+
testCase.verifyError(...
594+
@() out.set_bg_npoly_modQ(2), ...
595+
'sw_fitpowder:invalidinput');
596+
% set it to constant in modQ
597+
out.set_bg_npoly_modQ(0)
598+
expected_fitpow = testCase.default_fitpow;
599+
expected_fitpow.npoly_modQ = 0;
600+
expected_fitpow.modQ_cens = testCase.default_modQ_cens_1d;
601+
% remove a row
602+
expected_fitpow.params = expected_fitpow.params([1,3:end]);
603+
expected_fitpow.bounds = expected_fitpow.bounds([1,3:end],:);
604+
testCase.verify_results(out, expected_fitpow);
605+
end
606+
607+
function test_export_data_2d(testCase)
608+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
609+
testCase.fit_func, testCase.j1);
610+
data = out.export_data();
611+
testCase.verify_val(data, testCase.data_2d);
612+
end
613+
614+
function test_export_data_2d_with_filename(testCase)
615+
import matlab.unittest.fixtures.TemporaryFolderFixture
616+
fixture = testCase.applyFixture(TemporaryFolderFixture);
617+
tmp_file = fullfile(fixture.Folder,"data.mat");
618+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
619+
testCase.fit_func, testCase.j1);
620+
data = out.export_data(tmp_file);
621+
testCase.verify_val(data, testCase.data_2d);
622+
testCase.assertTrue(isfile(tmp_file));
623+
end
624+
625+
function test_export_data_1d(testCase)
626+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
627+
testCase.fit_func, testCase.j1);
628+
data = out.export_data();
629+
testCase.verify_val(data, testCase.data_1d_cuts);
630+
end
631+
632+
function test_set_bg_param_with_name_planar_bg(testCase)
633+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
634+
testCase.fit_func, testCase.j1);
635+
out.set_bg_parameters("Q1", -5);
636+
out.set_bg_parameters("E1", 5);
637+
expected_fitpow = testCase.default_fitpow;
638+
expected_fitpow.params(2:4) = [-5, 5,0]; % Q1, E1, E0
639+
testCase.verify_results(out, expected_fitpow);
640+
end
641+
642+
function test_set_bg_param_with_invalid_name(testCase)
643+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
644+
testCase.fit_func, testCase.j1);
645+
testCase.verifyError(...
646+
@() out.set_bg_parameters("Q3", -5), ...
647+
'sw_fitpowder:invalidinput');
648+
end
649+
650+
function test_set_bg_param_with_name_indep_bg_all_cuts(testCase)
651+
out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ...
652+
testCase.fit_func, testCase.j1, "independent", 1);
653+
bg_pars = [-5, 5];
654+
out.set_bg_parameters(["E1", "E0"], bg_pars);
655+
expected_fitpow = testCase.default_fitpow;
656+
expected_fitpow.params = [expected_fitpow.params(1);
657+
bg_pars(:); bg_pars(:); 1];
658+
expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:);
659+
testCase.verify_results(out, expected_fitpow);
660+
end
661+
662+
function test_fix_bg_param_with_array_of_names(testCase)
663+
out = sw_fitpowder(testCase.swobj, testCase.data_2d, ...
664+
testCase.fit_func, testCase.j1);
665+
bg_pars = 1:3;
666+
bg_labels = ["Q1","E1","E0"];
667+
out.set_bg_parameters(bg_labels, bg_pars);
668+
out.fix_bg_parameters(bg_labels);
669+
expected_fitpow = testCase.default_fitpow;
670+
expected_fitpow.params(2:4) = bg_pars;
671+
expected_fitpow.bounds(2:4, :) = repmat(bg_pars(:), 1,2);
672+
testCase.verify_results(out, expected_fitpow);
673+
end
674+
675+
function test_add_1Dcuts_withs_qs_specified(testCase)
676+
cuts = arrayfun(@(qs) struct('x', 1:3, 'y', 1:3, 'e', 1:3, 'qs', qs), ...
677+
4:5);
678+
out = sw_fitpowder(testCase.swobj, cuts, ...
679+
testCase.fit_func, testCase.j1);
680+
testCase.verify_results(out, testCase.default_fitpow);
681+
end
682+
453683
end
454684

455685
end

‎swfiles/+ndbase/lm4.m‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -223,14 +223,14 @@
223223
end
224224
% collect output
225225
pOpt = cost_func_wrap.get_bound_parameters(p);
226+
perr = zeros(size(pOpt));
226227
fVal = cost_val / ndof;
227228
if exit_flag > 0
228229
% converged on solution - calculate errors
229-
cov = pinv(hess) * 2.0 * fVal;
230-
perr = sqrt(diag(cov));
230+
cov = pinv(hess) * 2.0;
231+
perr(cost_func_wrap.ifree) = sqrt(diag(cov));
231232
else
232233
message = "Failed to converge in MaxIter";
233-
perr = zeros(size(pOpt));
234234
end
235235
end
236236

0 commit comments

Comments
 (0)