Class interface to fit powder spectra

[RichardWaiteSTFC]

Note includes commits from #167 and #163

To-Do:

• Fix bugs in ndbase.lm2 & 3 and ndbase.pso that give incorrect size error when used as optimizer~~
• 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]

Test Results

    4 files  ±  0    112 suites  +4   4m 40s ⏱️ - 4m 41s
  711 tests + 25    693 ✅ + 25  18 💤 ±0  0 ❌ ±0 
2 016 runs  +100  1 980 ✅ +100  36 💤 ±0  0 ❌ ±0 

Results for commit d72280c. ± Comparison against base commit 3e8b492.

♻️ This comment has been updated with latest results.

[RichardWaiteSTFC]

2D plot of initial guess

1D plot of initial guess

[mducle]

Looks good, just a couple of minor points. Also, you can delete the fitpow.m file...

But it needs more documentation and some basic tests... (maybe just make sure that the example syntax can run with 1 or 2 iteration so that it doesn't significantly increase the time to run tests...)

[mducle]

Sorry another doc request and also a request to add a "feature" from fitpow - of caching a previous powspec calculation (from a previous iteration). This is only useful for gradient methods like L-M but it saves evaluating the expensive powder spinwave calculations when the algorithm is varying the background or scaling parameters - it effectively saves around (n_bgpar/n_totalpars) fraction of iterations.

[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]

⚠️ 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 ⚠️
swfiles/@spinw/powspec.m	76.92%	9 Missing ⚠️
swfiles/sw_instrument.m	0.00%	9 Missing ⚠️
swfiles/+ndbase/lm3.m	0.00%	5 Missing ⚠️
swfiles/sw_mex.m	0.00%	3 Missing ⚠️

Additional details and impacted files

@@            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     

☔ View full report in Codecov by Sentry.
📢 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!