Class interface to fit powder spectra

[RichardWaiteSTFC]

Note includes commits from #167 and #163

To-Do:

• [x] Fix bugs in ndbase.lm2 & 3 and ndbase.pso that give incorrect size error when used as optimizer~~
• [x] Add capability to plot fit result

To test (data here)

% test spinw powder fitting
mnf2dat = load("C:\Users\Downloads\mnf2\mnf2.mat");

% init spinw object
J1 = -0.05;
J2 = 0.3;
D = 0.05;

mnf2 = spinw;
mnf2.genlattice('lat_const', [4.87 4.87 3.31], 'angle', [90 90 90]*pi/180, 'sym', 'P 42/m n m');
mnf2.addatom('r', [0 0 0], 'S', 2.5, 'label', 'MMn2', 'color', 'b')
mnf2.gencoupling('maxDistance', 5)
mnf2.addmatrix('label', 'J1', 'value', J1, 'color', 'red');
mnf2.addmatrix('label', 'J2', 'value', J2, 'color', 'green');
mnf2.addcoupling('mat', 'J1', 'bond', 1)
mnf2.addcoupling('mat', 'J2', 'bond', 2)
mnf2.addmatrix('label', 'D', 'value', diag([0 0 D]), 'color', 'black');
mnf2.addaniso('D')
mnf2.genmagstr('mode', 'direct', 'S', [0 0; 0 0; 1 -1])

% define fit_func
fit_func =  @(obj, p) matparser(obj, 'param', p, 'mat', {'J1', 'J2', 'D(3,3)'}, 'init', true);

% Resolution (from PyChop)
Ei = 20; % 20
eres = @(en) 512.17*sqrt((Ei-en).^3 .* ( 8.26326e-10*(0.169+0.4*(Ei./(Ei-en)).^1.5).^2 + 2.81618e-11*(1.169+0.4*(Ei./(Ei-en)).^1.5).^2));
% Q-resolution (parameters for MARI)
e2k = @(en) sqrt( en .* (2*1.67492728e-27*1.60217653e-22) )./1.05457168e-34./1e10;
L1 = 11.8;  % Moderator to Sample distance in m
L2 = 2.5;   % Sample to detector distance in m
ws = 0.05;  % Width of sample in m
wm = 0.12;  % Width of moderator in m
wd = 0.025; % Width of detector in m
ki = e2k(Ei);
a1 = ws/L1; % Angular width of sample seen from moderator
a2 = wm/L1; % Angular width of moderator seen from sample
a3 = wd/L2; % Angular width of detector seen from sample
a4 = ws/L2; % Angular width of sample seen from detector
dQ = 2.35 * sqrt( (ki*a1)^2/12 + (ki*a2)^2/12 + (ki*a3)^2/12 + (ki*a4)^2/12 );

%% 1D cut example

% make cuts
q0 = [0.7, 1.3, 2.0]; % , 2.6];
clear cts;
for iq = 1:numel(q0)
    cts(iq) = cutpow(mnf2dat, q0(iq)-0.1, q0(iq)+0.1, 1.0);
end

% planar bg (fix slopes so constant)
% -------------------------------------

fitpow = sw_fitpowder(mnf2, cts,  fit_func, [J1 J2 D], 'planar');
fitpow.set_caching(true);
fitpow.powspec_args.dE = eres;
fitpow.powspec_args.fastmode = false;
fitpow.powspec_args.neutron_output = true;
fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei);
fitpow.estimate_constant_background()
fitpow.estimate_scale_factor();
fitpow.set_model_parameter_bounds(1:3, [-1 0 -0.2], [0 1 0.2]) % Force J1 to be ferromagnetic to agree with structure
% fix slope and set initial const background
fitpow.set_bg_parameter_bounds(3, 0.0, []) % Set lb of constant bg = 0
fitpow.fix_bg_parameters(1:2); % fix slopes of background to 0
fitpow.cost_function = "Rsq";  % "Rsq" or "chisq"
fitpow.optimizer = @ndbase.simplex;
tic;
result = fitpow.fit();
toc;

% fitpow.plot_result(fitpow.params)  % initial guess
[pfit, cost_val, stat] = result{:};
fitpow.plot_result(pfit)  % result

% independent bg (fix slopes so constant)
% ---------------------------------------------

fitpow = sw_fitpowder(mnf2, cts,  fit_func, [J1 J2 D], 'independent');
fitpow.set_caching(true);
fitpow.powspec_args.dE = eres;
fitpow.powspec_args.fastmode = false;
fitpow.powspec_args.neutron_output = true;
fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei);
fitpow.estimate_constant_background();
fitpow.estimate_scale_factor();
fitpow.set_model_parameter_bounds(1:3, [-1 0 -0.2], [0 1 0.2]) % Force J1 to be ferromagnetic to agree with structure
% fix slope and set initial const background
fitpow.set_bg_parameter_bounds(2, 0.0, []) % Set lb of constant bg = 0
fitpow.fix_bg_parameters(1); % fix slopes of background to 0
fitpow.cost_function = "Rsq";  % "Rsq" or "chisq"
fitpow.optimizer = @ndbase.simplex;
tic;
result = fitpow.fit();
toc;

% fitpow.plot_result(fitpow.params)  % initial guess
[pfit, cost_val, stat] = result{:};
fitpow.plot_result(pfit)  % result

%% 2D cut example

fitpow = sw_fitpowder(mnf2, mnf2dat,  fit_func, [J1 J2 D]);
% fitpow.liveplot_interval = 2;
fitpow.crop_energy_range(2.0, 8.0);
fitpow.crop_q_range(0.4, 3.5);
fitpow.powspec_args.dE = eres;
fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei);
fitpow.estimate_constant_background()
fitpow.estimate_scale_factor();
fitpow.set_model_parameter_bounds(1:3, [-1 0 -0.2], [0 1 0.2]) % Force J1 to be ferromagnetic to agree with structure
fitpow.set_bg_parameter_bounds(3, 0.0, []) % Set lb of constant bg = 0
fitpow.fix_bg_parameters(1:2); % fix slopes of background to 0
fitpow.cost_function = "Rsq";  % "Rsq" or "chisq"
fitpow.optimizer = @ndbase.simplex;
tic;
result = fitpow.fit('MaxIter', 1);
toc;

% plot
[pfit,cost_val,stat] = result{:};
fitpow.plot_result(pfit, 16, 'EdgeAlpha', 0.9, 'LineWidth', 1) % ,'EdgeColor', 'k')

[github-actions[bot]]

Test Results

    4 files  ±  0    112 suites  +4   4m 40s :stopwatch: - 4m 41s   711 tests + 25    693 :white_check_mark: + 25  18 :zzz: ±0  0 :x: ±0  2 016 runs  +100  1 980 :white_check_mark: +100  36 :zzz: ±0  0 :x: ±0 

Results for commit d72280cd. ± Comparison against base commit 3e8b492d.

:recycle: This comment has been updated with latest results.

[RichardWaiteSTFC]

2D plot of initial guess

1D plot of initial guess

[RichardWaiteSTFC]

Thanks @mducle for the suggestion of caching the spinwave calculation, it does provide considerable speedup when using e.g. the lm3 minimiser. Here are timings for fitting the 3 x1D cuts in the PR description

Background    With/Without Cache
planar             142/223 s
indep.             104/247 s

As expected the speedup is greater for the independent backgrounds as there are more background parameters.

[RichardWaiteSTFC]

Thanks for the detailed review! FYI having added fastmode and neutron_output to powspec I ran 1 iteration of the simplex minimiser on the 2D dataset - here are the rough timings:

Time (s)   fastmode    neutron_output
 45.6        1            1
 60.1        0            1
 69.3        0            0

[codecov-commenter]

:warning: Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 65.13158% with 106 lines in your changes missing coverage. Please review.

Project coverage is 42.11%. Comparing base (735d30a) to head (d72280c). Report is 31 commits behind head on master.

Files	Patch %	Lines
swfiles/sw_fitpowder.m	66.66%	80 Missing :warning:
swfiles/@spinw/powspec.m	76.92%	9 Missing :warning:
swfiles/sw_instrument.m	0.00%	9 Missing :warning:
swfiles/+ndbase/lm3.m	0.00%	5 Missing :warning:
swfiles/sw_mex.m	0.00%	3 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #174 +/- ## ========================================== + Coverage 40.51% 42.11% +1.60% ========================================== Files 240 240 Lines 15981 16079 +98 ========================================== + Hits 6474 6771 +297 + Misses 9507 9308 -199 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[mducle]

@RichardWaiteSTFC Looks good but could you add the changes in the following diff, please:

spinwfitcacheinval.txt

The changes to lm3.m is a quick hack to make sure that "fixed" parameters stay fixed... we'll need to revisit this and rewrite lm3 properly but I need this for the current fitting stuff with the user.

The two lines to invalidate the cache are needed because add_data and apply_energy_mask changes the shape of the y and e members such that calc_background will calculate a bg for the new shape but if the cache is not invalidated it will return the old y values (shape) and give a shape mismatch error.

Edit: This is mainly a problem if you call add_data or apply_energy_mask after a call to estimate_scale_factor as that does a spinwave calculation and saves it to the cache.

[RichardWaiteSTFC]

@RichardWaiteSTFC Looks good but could you add the changes in the following diff, please:

spinwfitcacheinval.txt

The changes to lm3.m is a quick hack to make sure that "fixed" parameters stay fixed... we'll need to revisit this and rewrite lm3 properly but I need this for the current fitting stuff with the user.

The two lines to invalidate the cache are needed because add_data and apply_energy_mask changes the shape of the y and e members such that calc_background will calculate a bg for the new shape but if the cache is not invalidated it will return the old y values (shape) and give a shape mismatch error.

Edit: This is mainly a problem if you call add_data or apply_energy_mask after a call to estimate_scale_factor as that does a spinwave calculation and saves it to the cache.

Thanks Duc - there is already a method called clear_cache - would it be OK to call that rather than touch the obj.model_params_cached (which is currently private)? It deletes the ycalc array so there's a bit of overhead re-assigning the memory but maybe that would have to happen anyway is the shape changes? I can change the behaviour of that anyway if you prefer?

[mducle]

Doh! Sorry I should have read the code more carefully. Yes, using clear_cache is much better!

[RichardWaiteSTFC]

Doh! Sorry I should have read the code more carefully. Yes, using clear_cache is much better!

No worries, it was a good spot about cropping the data and the cache!