From 1cdd36b71fc42a17e3a5c5608168f82e2e22f6b0 Mon Sep 17 00:00:00 2001 From: Duc Le Date: Fri, 19 Jan 2024 01:26:25 +0000 Subject: [PATCH 01/43] Add fitpow method for fitting powder data Change powspec to call spinwave() once to take advantage of chunking or mex threading Refactor sw_freemem() to poll memory max only every 10s Change sw_egrid() to not calculate Bose factor for T=0 (default) The above changes save ~30% time per iteration in a powder fit (now spinwave() eval is ~50% of time per iteration, was ~30%) --- swfiles/@spinw/fitpow.m | 472 +++++++++++++++++++++++++++++++++++++++ swfiles/@spinw/powspec.m | 67 +++--- swfiles/sw_freemem.m | 14 ++ 3 files changed, 519 insertions(+), 34 deletions(-) create mode 100644 swfiles/@spinw/fitpow.m diff --git a/swfiles/@spinw/fitpow.m b/swfiles/@spinw/fitpow.m new file mode 100644 index 000000000..fba9585f2 --- /dev/null +++ b/swfiles/@spinw/fitpow.m @@ -0,0 +1,472 @@ +function fitsp = fitpow(obj, data, varargin) +% fits experimental powder spin wave data +% +% ### Syntax +% +% `fitsp = fitpow(obj, data, Name, Value)` +% +% ### Description +% `fitsp = fitpow(obj, data, Name, Value)` takes a set of energy cuts of powder +% spin wave spectrum and fits this. +% +% ### Input Arguments +% +% `obj` +% : [spinw] object. +% +% `data` +% : input data cuts which may be an array of any of: +% * a Horace 1D sqw or d1d object +% * a struct with fields 'x', 'y', 'e', 'qmin', 'qmax' +% * a struct with fields 'IX', 'qmin', 'qmax' +% * a cell array {x y e qmin qmax} +% * a cell array {IX_dataset_1d qmin qmax} +% where x is a vector of energy transfer, +% y is a vector of intensity, +% e is a vector of uncertainties (errorbar) +% qmin is a scalar denoting the lower Q value of the cut +% qmax is a scalar denoting the upper Q value of the cut +% IX / IX_dataset_1d is a Horace IX_dataset_1D object +% Instead of the pair of scalars qmin, qmax, you specify a +% 2-vector: [qmin qmax] under the key 'qrange', e.g. +% * a struct with fields 'x', 'y', 'e', 'qrange' +% * a cell array {x y e [qmin qmax]} +% +% ### Name-Value Pair Arguments +% +% `'func'` +% : Function to change the Hamiltonian in `obj`, it needs to have the +% following header: +% ``` +% obj = @func(obj, x) +% ``` +% +% `'xmin'` +% : Lower limit of the optimisation parameters, optional. +% +% `'xmax'` +% : Upper limit of the optimisation parameters, optional. +% +% `'x0'` +% : Starting value of the optimisation parameters. If empty +% or undefined, random values are used within the given limits. +% +% `'optimizer'` +% : String that determines the type of optimizer to use, possible values: +% * `'pso'` Particle-swarm optimizer, see [ndbase.pso], +% default. +% * `'simplex'` Matlab built-in simplex optimizer, see [fminsearch](www.mathworks.ch/help/matlab/ref/fminsearch.html). +% +% `'nRun'` +% : Number of consecutive fitting runs, each result is saved in the +% output `fitsp.x` and `fitsp.R` arrays. If the Hamiltonian given by the +% random `x` parameters is incompatible with the ground state, +% those `x` values will be omitted and new random `x` values will be +% generated instead. Default value is 1. +% +% `'nMax'` +% : Maximum number of runs, including the ones that produce error +% (due to incompatible ground state). Default value is 1000. +% +% `'hermit'` +% : Method for matrix diagonalization, for details see [spinw.spinwave]. +% +% `'epsilon'` +% : Small number that controls wether the magnetic structure is +% incommensurate or commensurate, default value is $10^{-5}$. +% +% `'imagChk'` +% : Checks that the imaginary part of the spin wave dispersion is +% smaller than the energy bin size. +% If false, will not check +% If 'penalize' will apply a penalty to iterations that yield imaginary modes +% If true, will stop the fit if an iteration gives imaginary modes +% Default is `penalize`. +% +% Parameters for visualizing the fit results: +% +% `'plot'` +% : If `true`, the measured dispersion is plotted together with the +% fit. Default is `true`. +% +% `'iFact'` +% : Factor of the plotted simulated spin wave intensity (red +% ellipsoids). +% +% `'lShift'` +% : Vertical shift of the `Q` point labels on the plot. +% +% Optimizer options: +% +% `'TolX'` +% : Minimum change of` x` when convergence reached, default +% value is $10^{-4}$. +% +% `'TolFun'` +% : Minimum change of the R value when convergence reached, +% default value is $10^{-5}$. +% +% `'MaxFunEvals'` +% : Maximum number of function evaluations, default value is +% $10^7$. +% +% `'MaxIter'` +% : Maximum number of iterations for the [ndbse.pso] optimizer. +% Default value is 20. +% +% ### Output Arguments +% +% Output `fitsp` is struct type with the following fields: +% * `obj` Copy of the input `obj`, with the best fitted +% Hamiltonian parameters. +% * `x` Final values of the fitted parameters, dimensions are +% $[n_{run}\times n_{par}]$. The rows of `x` are sorted according +% to increasing R values. +% * `redX2` Reduced $\chi^2_\eta$ value, goodness of the fit stored in a column +% vector with $n_{run}$ number of elements, sorted in increasing +% order. $\chi^2_\eta$ is defined as: +% +% $\begin{align} +% \chi^2_\eta &= \frac{\chi^2}{\eta},\\ +% \eta &= n-m+1, +% \end{align}$ +% where \\eta is the degree of freedom, $n$ number of +% observations and $m$ is the number of fitted parameters. +% +% * `exitflag` Exit flag of the `fminsearch` command. +% * `output` Output of the `fminsearch` command. +% +% {{note As a rule of thumb when the variance of the measurement error is +% known a priori, \\chi$^2_\eta$\\gg 1 indicates a poor model fit. A +% \\chi$^2_\eta$\\gg 1 indicates that the fit has not fully captured the +% data (or that the error variance has been underestimated). In principle, +% a value of \\chi$^2_\eta$= 1 indicates that the extent of the match +% between observations and estimates is in accord with the error variance. +% A \\chi$^2_\eta$ < 1 indicates that the model is 'over-fitting' the data: +% either the model is improperly fitting noise, or the error variance has +% been overestimated.}} +% +% Any other option used by [spinw.spinwave] function are also accepted. +% +% ### See Also +% +% [spinw.spinwave] \| [spinw.matparser] \| [sw_egrid] \| [sw_neutron] +% +pref = swpref; +tid0 = pref.tid; + +if nargin < 2 + error('spinw:fitpow:MissingData','No input data cuts given') +end + +inpForm.fname = {'epsilon' 'xmin' 'xmax' 'x0' 'func' 'fibo' 'nRand'}; +inpForm.defval = {1e-5 [] [] [] [] true 100}; +inpForm.size = {[1 1] [1 -3] [1 -4] [1 -5] [1 1] [1 1] [1 1]}; +inpForm.soft = {true true true true false false false}; + +inpForm.fname = [inpForm.fname {'formfact' 'formfactfun' 'conv_func' 'nQ' 'intensity_scale'}]; +inpForm.defval = [inpForm.defval {false @sw_mff @(s)s 10 0}]; +inpForm.size = [inpForm.size {[1 -2] [1 1] [1 1] [1 1] [1 -9]}]; +inpForm.soft = [inpForm.soft {false false false false true}]; + +inpForm.fname = [inpForm.fname {'tolx' 'tolfun' 'maxfunevals' 'nRun' 'plot'}]; +inpForm.defval = [inpForm.defval {1e-4 1e-5 1e3 1 true}]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 1]}]; +inpForm.soft = [inpForm.soft {false false false false false}]; + +inpForm.fname = [inpForm.fname {'nMax' 'hermit' 'iFact' 'lShift' 'optimizer'}]; +inpForm.defval = [inpForm.defval {1e3 true 1/30 0 'pso' }]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 -7] }]; +inpForm.soft = [inpForm.soft {false false false false false }]; + +inpForm.fname = [inpForm.fname {'maxiter' 'sw' 'optmem' 'tid' 'imagChk'}]; +inpForm.defval = [inpForm.defval {20 1 0 tid0 'penalize'}]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 -8] }]; +inpForm.soft = [inpForm.soft {false false false false false }]; + +param = sw_readparam(inpForm, varargin{:}); + +% number of parameters (length of x) +nPar = max(max(length(param.xmin),length(param.xmax)),length(param.x0)); +nRun = param.nRun; + +% Parses the input data +[lindat, data] = parse_cuts(data); + +% setup parameters for spinw.spinwave() +param0 = param; +%param0.plot = false; +param0.tid = 0; + +% Handle scale factor +param0.fit_separate = false; +param0.fit_together = false; +param0.scale_factor = 1.0; +if ischar(param.intensity_scale) || isstring(param.intensity_scale) + if strcmp(param.intensity_scale, 'fit_separate') + param0.fit_separate = true; + elseif strcmp(param.intensity_scale, 'fit_together') + param0.fit_together = true; + end +elseif isfloat(param.intensity_scale) + param0.scale_factor = param.intensity_scale; +end + +% Parses the imagChk parameter - we need to set it to true/false for spinw.spinwave() +if strncmp(param.imagChk, 'penalize', 1) + param0.imagChk = true; + param0.penalize_imag = true; +else + param0.imagChk = logical(param.imagChk(1)); + param0.penalize_imag = false; +end + +optfunname = sprintf('ndbase.%s', param.optimizer); +assert(~isempty(which(optfunname)), 'spinw:fitpow:invalidoptimizer', 'Optimizer %s not recognised', param.optimizer); +optimizer = str2func(optfunname); + +x = zeros(nRun,nPar); +redX2 = zeros(nRun,1); + +sw_timeit(0, 1, param.tid,'Fitting powder spin wave spectra'); + +idx = 1; + +while idx <= nRun + if ~isempty(param.x0) + x0 = param.x0; + elseif isempty(param.xmax) && isempty(param.xmin) + x0 = rand(1,nPar); + elseif isempty(param.xmax) + x0 = param.xmin + rand(1,nPar); + elseif isempty(param.xmin) + x0 = param.xmax - rand(1,nPar); + else + x0 = rand(1,nPar).*(param.xmax-param.xmin)+param.xmin; + end + + [x(idx,:),~, output(idx)] = optimizer(lindat, @(x,p)pow_fitfun(obj, lindat, data, param.func, p, param0), x0, ... + 'lb', param.xmin, 'ub', param.xmax, ... + 'TolX', param.tolx, 'TolFun', param.tolfun, 'MaxIter', param.maxiter); + redX2(idx) = output.redX2; + + idx = idx + 1; + sw_timeit(idx/nRun*100,0,param.tid); +end + +sw_timeit(100,2,param.tid); + +% Sort results +[redX2, sortIdx] = sort(redX2); +x = x(sortIdx,:); + +% set the best fit to the spinw object +param.func(obj,x(1,:)); + +% Store all output in a struct variable. +fitsp.obj = copy(obj); +fitsp.x = x; +fitsp.redX2 = redX2; +%fitsp.exitflag = exitflag; +fitsp.output = output; + +end + +function out = horace1d_tostruct(dat) + assert(isa(dat, 'd1d') || isa(dat, 'sqw'), 'spinw:fitpow:invalidinput', 'Input must be a Horace 1D object, cell or struct'); + if isa(dat, 'sqw') + assert(dat.dimensions == 1, 'spinw:fitpow:invalidinput', 'Input SQW object must be 1D'); + dd = dat.data; + else + dd = dat; + end + assert(strcmp(dat.ulabel{1}, '|Q|'), 'spinw:fitpow:invalidinput', 'Input Horace object is not a powder cut'); + out = struct('x', dd.p{1}, 'y', dd.s, 'e', dd.e, 'qmin', dd.iint(1,1), 'qmax', dd.iint(2,1)); +end + +function out = verify_cut_struct(dat) + if all(isfield(dat, {'x', 'y', 'e', 'qmin', 'qmax'})) + out = dat; + elseif all(isfield(dat, {'x', 'y', 'e', 'qrange'})) + out = struct('x', dat.x, 'y', dat.y, 'e', dat.e, 'qmin', dat.qrange(1), 'qmax', dat.qrange(2)); + elseif all(isfield(dat, {'IX', 'qmin', 'qmax'})) + out = struct('x', dat.IX.x, 'y', dat.IX.signal, 'e', dat.IX.error, 'qmin', dat.qmin, 'qmax', dat.qmax); + elseif all(isfield(dat, {'IX', 'qrange'})) + out = struct('x', dat.IX.x, 'y', dat.IX.signal, 'e', dat.IX.error, 'qmin', dat.qrange(1), 'qmax', dat.qrange(2)); + else + error('spinw:fitpow:invalidinput', 'Input structure must have fields (x,y,e) or (IX) and (qmin,qmax) or (qrange)'); + end +end + +function out = xyecell_tostruct(dat) + errmsg = 'Cell input should consist of IX_dataset objects or numeric arrays'; + if numel(dat) > 3 + assert(all(cellfun(@isnumeric, dat)), 'spinw:fitpow:invalidinput', errmsg); + xye = dat(1:3); + remd = dat(4:end); + elseif numel(dat) > 1 + assert(isa(dat{1}, 'IX_dataset_1d'), 'spinw:fitpow:invalidinput', errmsg); + xye = {dat{1}.x, dat{1}.signal, dat{1}.error}; + remd = dat(2:end); + end + if numel(remainder) == 2 + out = struct('x', xye{1}, 'y', xye{2}, 'e', xye{3}, 'qmin', remd{1}, 'qmax', remd{2}); + else + out = struct('x', xye{1}, 'y', xye{2}, 'e', xye{3}, 'qmin', remd{1}(1), 'qmax', remd{1}(2)); + end +end + +function [lindat, outstruct] = parse_cuts(input) + %outstruct = struct(); + if iscell(input) + for ii = 1:numel(input) + if isstruct(input{ii}) + outstruct(ii) = verify_cut_struct(input{ii}); + elseif iscell(input{ii}) + outstruct(ii) = xyecell_tostruct(input{ii}); + else + outstruct(ii) = horace1d_tostruct(input{ii}); + end + end + elseif isstruct(input) + for ii = 1:numel(input) + outstruct(ii) = verify_cut_struct(input(ii)); + end + else + for ii = 1:numel(input) + outstruct(ii) = horace1d_tostruct(input(ii)); + end + end + for ii = 1:numel(outstruct) + xx = outstruct(ii).x(:).'; + df = diff(xx) ./ 2; + outstruct(ii).evect = [xx(1)-df(1), xx(1:end-1) + df, xx(end) + df(end)]; + end + lindat = linearise_data(outstruct); +end + +function lindat = linearise_data(outstruct) + % Linearise the data + [x, y, e] = deal([], [], []); + is_same_e = true; + for ii = 1:numel(outstruct) + x = [x; ii+([1:numel(outstruct(ii).y)]/1e9)]; + y = [y; outstruct(ii).y]; + e = [e; outstruct(ii).e]; + if ii == 1 + evec = outstruct(ii).x; + else + if ~is_same_e && (numel(evec) ~= numel(outstruct(ii).x) ... + || sum(abs(evec(:) - outstruc(ii).x(:))) > 1e-5) + is_same_e = false; + end + end + end + lindat = struct('x', x, 'y', y, 'e', e, 'is_same_e', is_same_e); +end + + +function yCalc = pow_fitfun(obj, lindat, data, parfunc, x, param) +% calculates the powder spin wave spectrum to directly compare to data +% +% yCalc = SW_FITFUN(obj, data, swfunc, x, param) +% +% swobj spinw object +% data Structure contains the experimental data +% swfunc Function to change the Hamiltonian in obj, of the form: +% obj = swfunc(obj,x); +% x Actual parameter values. +% param Parameters for the spinwave function. + +persistent swp hf lines; +if isempty(swp) + swp = swpref; +end +fid0 = swp.fid; +swp.fid = 0; +cleanupObj = onCleanup(@()swpref.setpref('fid', fid0)); + +parfunc(obj,x); + +if param.plot + if isempty(hf) || ~isvalid(hf) + hf = findobj('Tag', 'fitpow_plot'); + end + if isempty(hf) || ~isvalid(hf) + hf = figure('Tag', 'fitpow_plot'); + for ii = 1:numel(data) + subplot(1, numel(data), ii); hold all; + errorbar(data(ii).x, data(ii).y, data(ii).e, '.k'); + end + lines = []; + end + set(0, 'CurrentFigure', hf); +end + +fitstruc = data; +sf = param.scale_factor; +nD = numel(data); +nY = 0; + +if lindat.is_same_e + nY = nD * numel(data(1).x); + qq = zeros(nD * param.nQ); + for ii = 1:nD + i0 = (ii - 1) * param.nQ + 1; + qq(i0:(i0 + param.nQ - 1)) = linspace(data(ii).qmin, data(ii).qmax, param.nQ); + end + spec = obj.powspec(qq, 'Evect', data(1).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... + 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... + 'imagChk', param.imagChk); + spec = param.conv_func(spec); + for ii = 1:nD + i0 = (ii - 1) * param.nQ + 1; + fitstruc(ii).y = sum(spec.swConv(:, i0:(i0 + param.nQ - 1)), 2); + end +else + for ii = 1:nD + qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); + spec = obj.powspec(qq, 'Evect', data(ii).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... + 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... + 'imagChk', param.imagChk); + spec = param.conv_func(spec); + fitstruc(ii).y = sum(spec.swConv, 2); + nY = nY + numel(fitstruc(ii).y); + end +end + +yCalc = zeros(nY, 1); +i0 = 1; +for ii = 1:numel(data) + if param.fit_separate + sf = 1.0; + ssf = fminsearch(@(x)sum(abs(data(ii).y(:) - x*fitstruc(ii).y(:)), 'omitnan'), max(data(ii).y) / max(fitstruc(ii).y)); + if ssf < 1e5 + fitstruc(ii).y = ssf * fitstruc(ii).y; + end + end + i1 = i0 + numel(fitstruc(ii).y) - 1; + yCalc(i0:i1) = fitstruc(ii).y; + i0 = i1 + 1; +end + +if param.fit_together + sf = fminsearch(@(x)sum(abs(lindat.y(:) - x*yCalc(:)), 'omitnan'), max(l0.y) / max(yCalc)); +end + +if sf ~= 1.0 + yCalc = sf * yCalc; +end + +if param.plot + for ii = 1:numel(data) + subplot(1, numel(data), ii); + try %#ok + delete(lines(ii)); + end + lines(ii) = plot(fitstruc(ii).x, sf * fitstruc(ii).y, '-r'); + end + drawnow; +end + +end diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index 6f6e9fdb4..97b8659d7 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -246,7 +246,6 @@ end nQ = numel(hklA); -powSpec = zeros(max(1,nE),nQ); fprintf0(fid,'Calculating powder spectra...\n'); @@ -260,6 +259,7 @@ sw_timeit(0,1,param.tid,'Powder spectrum calculation'); +nk = param.nRand; if param.fibo % apply the Fibonacci numerical integration on a sphere % according to J. Phys. A: Math. Gen. 37 (2004) 11591 @@ -279,11 +279,13 @@ QF(1,:) = cos(theta).*sin(phi); QF(2,:) = cos(theta).*cos(phi); + nk = F; end % lambda value for SCGA, empty will make integration in first loop specQ.lambda = []; +hkl = zeros(3, nk * nQ); for ii = 1:nQ if param.fibo Q = QF*hklA(ii); @@ -291,40 +293,37 @@ rQ = randn(3,param.nRand); Q = bsxfun(@rdivide,rQ,sqrt(sum(rQ.^2)))*hklA(ii); end - hkl = (Q'*obj.basisvector)'/2/pi; - - switch funIdx - case 0 - % general function call allow arbitrary additional parameters to - % pass to the spectral calculation function - warnState = warning('off','sw_readparam:UnreadInput'); - specQ = param.specfun(obj,hkl,varargin{:}); - warning(warnState); - case 1 - % @spinwave - specQ = spinwave(obj,hkl,struct('fitmode',true,'notwin',true,... - 'Hermit',param.hermit,'formfact',param.formfact,... - 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... - 'optmem',param.optmem,'tid',0,'fid',0),'noCheck'); - - case 2 - % @scga - specQ = scga(obj,hkl,struct('fitmode',true,'formfact',param.formfact,... - 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... - 'fid',0,'lambda',specQ.lambda,'nInt',param.nInt,'T',param.T,... - 'plot',false),'noCheck'); - end - - specQ = sw_neutron(specQ,'pol',false); - specQ.obj = obj; - % use edge grid by default - specQ = sw_egrid(specQ,struct('Evect',param.Evect,'T',param.T,'binType',param.binType,... - 'imagChk',param.imagChk,'component',param.component),'noCheck'); - powSpec(:,ii) = sum(specQ.swConv,2)/param.nRand; - sw_timeit(ii/nQ*100,0,param.tid); + i0 = (ii-1) * nk + 1; + hkl(:, i0:(i0+nk-1)) = (Q'*obj.basisvector)'/2/pi; +end + +switch funIdx + case 0 + % general function call allow arbitrary additional parameters to + % pass to the spectral calculation function + warnState = warning('off','sw_readparam:UnreadInput'); + specQ = param.specfun(obj,hkl,varargin{:}); + warning(warnState); + case 1 + % @spinwave + specQ = spinwave(obj,hkl,struct('fitmode',true,'notwin',true,... + 'Hermit',param.hermit,'formfact',param.formfact,... + 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... + 'optmem',param.optmem,'tid',param.tid,'fid',0),'noCheck'); + case 2 + % @scga + specQ = scga(obj,hkl,struct('fitmode',true,'formfact',param.formfact,... + 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... + 'fid',0,'lambda',specQ.lambda,'nInt',param.nInt,'T',param.T,... + 'plot',false),'noCheck'); end -sw_timeit(100,2,param.tid); +specQ = sw_neutron(specQ,'pol',false); +specQ.obj = obj; +% use edge grid by default +specQ = sw_egrid(specQ,struct('Evect',param.Evect,'T',param.T,'binType',param.binType,... + 'imagChk',param.imagChk,'component',param.component),'noCheck'); +powSpec = squeeze(sum(reshape(specQ.swConv, max(1, nE), nk, nQ), 2) / param.nRand); fprintf0(fid,'Calculation finished.\n'); @@ -377,4 +376,4 @@ F = num(end-1); F1 = num(end-2); end -end \ No newline at end of file +end diff --git a/swfiles/sw_freemem.m b/swfiles/sw_freemem.m index 8e1be2269..d2327a7a5 100644 --- a/swfiles/sw_freemem.m +++ b/swfiles/sw_freemem.m @@ -16,7 +16,21 @@ % `mem` % : Size of free memory in bytes. % +persistent t0 m0; +if isempty(t0) || isempty(m0) + m0 = get_free(); + t0 = tic; +else + if isempty(t0) || toc(t0) > 10 % poll only every 10 seconds + m0 = get_free(); + t0 = tic; + end +end +mem = m0; + +end +function mem = get_free() mem = 0; try %#ok From e7ed1673253ad89a7513285a1ccaeb51de50e9de Mon Sep 17 00:00:00 2001 From: Duc Le Date: Sun, 21 Jan 2024 14:07:46 +0000 Subject: [PATCH 02/43] Add intensity and background to fitpow --- swfiles/@spinw/fitpow.m | 456 ++++++++++++++++++++++++++++------------ 1 file changed, 326 insertions(+), 130 deletions(-) diff --git a/swfiles/@spinw/fitpow.m b/swfiles/@spinw/fitpow.m index fba9585f2..09161dc81 100644 --- a/swfiles/@spinw/fitpow.m +++ b/swfiles/@spinw/fitpow.m @@ -1,19 +1,19 @@ function fitsp = fitpow(obj, data, varargin) % fits experimental powder spin wave data -% +% % ### Syntax -% +% % `fitsp = fitpow(obj, data, Name, Value)` -% +% % ### Description % `fitsp = fitpow(obj, data, Name, Value)` takes a set of energy cuts of powder % spin wave spectrum and fits this. -% +% % ### Input Arguments -% +% % `obj` % : [spinw] object. -% +% % `data` % : input data cuts which may be an array of any of: % * a Horace 1D sqw or d1d object @@ -33,97 +33,118 @@ % * a cell array {x y e [qmin qmax]} % % ### Name-Value Pair Arguments -% +% % `'func'` % : Function to change the Hamiltonian in `obj`, it needs to have the % following header: % ``` % obj = @func(obj, x) % ``` -% +% % `'xmin'` % : Lower limit of the optimisation parameters, optional. -% +% % `'xmax'` % : Upper limit of the optimisation parameters, optional. -% +% % `'x0'` % : Starting value of the optimisation parameters. If empty % or undefined, random values are used within the given limits. -% +% % `'optimizer'` % : String that determines the type of optimizer to use, possible values: +% * `'lm3'` (Default) Levenberg-Marquardt as implemented in Horace, see [ndbase.lm3], % * `'pso'` Particle-swarm optimizer, see [ndbase.pso], -% default. % * `'simplex'` Matlab built-in simplex optimizer, see [fminsearch](www.mathworks.ch/help/matlab/ref/fminsearch.html). -% +% % `'nRun'` % : Number of consecutive fitting runs, each result is saved in the % output `fitsp.x` and `fitsp.R` arrays. If the Hamiltonian given by the % random `x` parameters is incompatible with the ground state, % those `x` values will be omitted and new random `x` values will be % generated instead. Default value is 1. -% +% % `'nMax'` % : Maximum number of runs, including the ones that produce error % (due to incompatible ground state). Default value is 1000. -% +% +% `'nQ'` +% : Number of |Q| points to use for each cut to simulate cut thickness +% Default is 10. +% % `'hermit'` % : Method for matrix diagonalization, for details see [spinw.spinwave]. -% -% `'epsilon'` -% : Small number that controls wether the magnetic structure is -% incommensurate or commensurate, default value is $10^{-5}$. -% +% +% `'formfact'` +% : Whether to include the form factor in the calculation (Default: true) +% +% `'formfactfun'` +% : The function used to calculate the form factor (Default: sw_mff) +% +% `'conv_func'` +% : The convolution function used to broaden the data in |Q| and E. +% Default: no convolution is applied. +% % `'imagChk'` % : Checks that the imaginary part of the spin wave dispersion is -% smaller than the energy bin size. +% smaller than the energy bin size. % If false, will not check % If 'penalize' will apply a penalty to iterations that yield imaginary modes % If true, will stop the fit if an iteration gives imaginary modes % Default is `penalize`. -% -% Parameters for visualizing the fit results: -% +% +% `'intensity_scale'` +% : The scale factor for the intensity. Use a negative number for a fixed scale factor, +% positive for it to be fitted from the given initial value, +% and 0 for it to be fitted with an initial value guessed from the data. +% Default: -1 (fixed, no scale factor applied). +% +% `'background'` +% : The type of background to use. Either 'none', 'flat', 'linear', 'quadratic' +% +% `'background_par'` +% : The initial background parameters as a vector from highest order. +% E.g. if using ('background', 'quadratic'), then ('background_par', [a, b, c]) +% where the background used will be a*x.^2 + b*x + c. +% Default: [] (empty vector - initial parameters will be guessed from data). +% +% Parameters for monitoring the fitting: +% % `'plot'` -% : If `true`, the measured dispersion is plotted together with the -% fit. Default is `true`. -% -% `'iFact'` -% : Factor of the plotted simulated spin wave intensity (red -% ellipsoids). -% -% `'lShift'` -% : Vertical shift of the `Q` point labels on the plot. -% +% : If `true`, the data and fits are plotted as the fitting process progresses. +% Default is `true`. +% +% `'verbose'` +% : Whether to print information as the fit progresses +% Default is `true`. +% % Optimizer options: -% +% % `'TolX'` % : Minimum change of` x` when convergence reached, default % value is $10^{-4}$. -% +% % `'TolFun'` % : Minimum change of the R value when convergence reached, % default value is $10^{-5}$. -% -% `'MaxFunEvals'` -% : Maximum number of function evaluations, default value is -% $10^7$. -% -% `'MaxIter'` +% +% `'maxfunevals'` +% : Maximum number of function evaluations, default value is $10^3$. +% +% `'maxiter'` % : Maximum number of iterations for the [ndbse.pso] optimizer. % Default value is 20. -% +% % ### Output Arguments -% +% % Output `fitsp` is struct type with the following fields: % * `obj` Copy of the input `obj`, with the best fitted % Hamiltonian parameters. % * `x` Final values of the fitted parameters, dimensions are -% $[n_{run}\times n_{par}]$. The rows of `x` are sorted according +% $[n_{run}\times n_{par}]$. The rows of `x` are sorted according % to increasing R values. -% * `redX2` Reduced $\chi^2_\eta$ value, goodness of the fit stored in a column -% vector with $n_{run}$ number of elements, sorted in increasing +% * `redX2` Reduced $\chi^2_\eta$ value, goodness of the fit stored in a column +% vector with $n_{run}$ number of elements, sorted in increasing % order. $\chi^2_\eta$ is defined as: % % $\begin{align} @@ -135,7 +156,7 @@ % % * `exitflag` Exit flag of the `fminsearch` command. % * `output` Output of the `fminsearch` command. -% +% % {{note As a rule of thumb when the variance of the measurement error is % known a priori, \\chi$^2_\eta$\\gg 1 indicates a poor model fit. A % \\chi$^2_\eta$\\gg 1 indicates that the fit has not fully captured the @@ -147,11 +168,12 @@ % been overestimated.}} % % Any other option used by [spinw.spinwave] function are also accepted. -% +% % ### See Also -% +% % [spinw.spinwave] \| [spinw.matparser] \| [sw_egrid] \| [sw_neutron] % + pref = swpref; tid0 = pref.tid; @@ -159,25 +181,30 @@ error('spinw:fitpow:MissingData','No input data cuts given') end -inpForm.fname = {'epsilon' 'xmin' 'xmax' 'x0' 'func' 'fibo' 'nRand'}; -inpForm.defval = {1e-5 [] [] [] [] true 100}; -inpForm.size = {[1 1] [1 -3] [1 -4] [1 -5] [1 1] [1 1] [1 1]}; -inpForm.soft = {true true true true false false false}; +inpForm.fname = {'xmin' 'xmax' 'x0' 'func' 'fibo' 'nRand' 'verbose'}; +inpForm.defval = {[] [] [] [] true 100 true}; +inpForm.size = {[1 -3] [1 -4] [1 -5] [1 1] [1 1] [1 1] [1 1]}; +inpForm.soft = {true true true false false false false}; + +inpForm.fname = [inpForm.fname {'formfact' 'formfactfun' 'conv_func' 'nQ'}]; +inpForm.defval = [inpForm.defval {true @sw_mff @(s)s 10}]; +inpForm.size = [inpForm.size {[1 -2] [1 1] [1 1] [1 1]}]; +inpForm.soft = [inpForm.soft {false false false false}]; -inpForm.fname = [inpForm.fname {'formfact' 'formfactfun' 'conv_func' 'nQ' 'intensity_scale'}]; -inpForm.defval = [inpForm.defval {false @sw_mff @(s)s 10 0}]; -inpForm.size = [inpForm.size {[1 -2] [1 1] [1 1] [1 1] [1 -9]}]; -inpForm.soft = [inpForm.soft {false false false false true}]; +inpForm.fname = [inpForm.fname {'intensity_scale' 'background' 'background_par'}]; +inpForm.defval = [inpForm.defval {-1 'linear' []}]; +inpForm.size = [inpForm.size {[1 1] [1 -9] [1 -10]}]; +inpForm.soft = [inpForm.soft {false false true}]; inpForm.fname = [inpForm.fname {'tolx' 'tolfun' 'maxfunevals' 'nRun' 'plot'}]; inpForm.defval = [inpForm.defval {1e-4 1e-5 1e3 1 true}]; inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 1]}]; inpForm.soft = [inpForm.soft {false false false false false}]; -inpForm.fname = [inpForm.fname {'nMax' 'hermit' 'iFact' 'lShift' 'optimizer'}]; -inpForm.defval = [inpForm.defval {1e3 true 1/30 0 'pso' }]; -inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 -7] }]; -inpForm.soft = [inpForm.soft {false false false false false }]; +inpForm.fname = [inpForm.fname {'nMax' 'hermit' 'optimizer'}]; +inpForm.defval = [inpForm.defval {1e3 true 'lm3' }]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 -7] }]; +inpForm.soft = [inpForm.soft {false false false }]; inpForm.fname = [inpForm.fname {'maxiter' 'sw' 'optmem' 'tid' 'imagChk'}]; inpForm.defval = [inpForm.defval {20 1 0 tid0 'penalize'}]; @@ -198,20 +225,6 @@ %param0.plot = false; param0.tid = 0; -% Handle scale factor -param0.fit_separate = false; -param0.fit_together = false; -param0.scale_factor = 1.0; -if ischar(param.intensity_scale) || isstring(param.intensity_scale) - if strcmp(param.intensity_scale, 'fit_separate') - param0.fit_separate = true; - elseif strcmp(param.intensity_scale, 'fit_together') - param0.fit_together = true; - end -elseif isfloat(param.intensity_scale) - param0.scale_factor = param.intensity_scale; -end - % Parses the imagChk parameter - we need to set it to true/false for spinw.spinwave() if strncmp(param.imagChk, 'penalize', 1) param0.imagChk = true; @@ -225,7 +238,7 @@ assert(~isempty(which(optfunname)), 'spinw:fitpow:invalidoptimizer', 'Optimizer %s not recognised', param.optimizer); optimizer = str2func(optfunname); -x = zeros(nRun,nPar); +%x = zeros(nRun,nPar); redX2 = zeros(nRun,1); sw_timeit(0, 1, param.tid,'Fitting powder spin wave spectra'); @@ -244,12 +257,77 @@ else x0 = rand(1,nPar).*(param.xmax-param.xmin)+param.xmin; end + xmin = param.xmin; + xmax = param.xmax; + + % Handle scale factor + if param.intensity_scale < 0 + param0.fit_scale = false; + else + if param.intensity_scale == 0 + sf = guess_intensity(obj, data, param0); + else + sf = param.intensity_scale; + end + param0.fit_scale = true; + x0 = [x0 sf]; + if ~isempty(xmin) + xmin = [xmin 0]; + end + if ~isempty(xmax) + xmax = [xmax 10*sf]; + end + end + + % Handle background + param0.bkgfun = []; + if ~strcmp(param.background, 'none') + xb = param.background_par; + if strcmp(param.background, 'flat') + param0.bkgfun = @(x, p) ones(numel(x),1)*p; + if isempty(xb) + xb = min(lindat.y); + end + xnx = 0; + xmx = max(lindat.y) / 2; + elseif strcmp(param.background, 'linear') + param0.bkgfun = @(x, p) x*p(1) + p(2); + if isempty(xb) + xb = guess_bkg(data, 'linear'); + end + xnx = [-abs(xb(1))*2, 0]; + xmx = [abs(xb(1))*2, max(lindat.y)/2]; + elseif strcmp(param.background, 'quadratic') + param0.bkgfun = @(x, p) x.^2*p(1) + x*p(2) + p(3); + if isempty(xb) + xb = guess_bkg(data, 'quadratic'); + end + xnx = [-abs(xb(1:2))*2, 0]; + xmx = [abs(xb(1:2))*2, max(lindat.y)/2]; + else + error('spinw:fitpow:unknownbackground', 'Background type %s not supported', param.background); + end + xb(isnan(xb)) = 0; + x0 = [x0 xb]; %#ok<*AGROW> + param0.nbkgpar = numel(xb); + if ~isempty(xmin) + xmin = [xmin xnx]; + end + if ~isempty(xmax) + xmax = [xmax xmx]; + end + end + % Only cache powder calcs for gradient methods + if strncmp(param.optimizer, 'lm', 2) + pow_fitfun('clear_cache'); + else + pow_fitfun('no_cache'); + end [x(idx,:),~, output(idx)] = optimizer(lindat, @(x,p)pow_fitfun(obj, lindat, data, param.func, p, param0), x0, ... - 'lb', param.xmin, 'ub', param.xmax, ... - 'TolX', param.tolx, 'TolFun', param.tolfun, 'MaxIter', param.maxiter); + 'lb', xmin, 'ub', xmax, 'TolX', param.tolx, 'TolFun', param.tolfun, 'MaxIter', param.maxiter); redX2(idx) = output.redX2; - + idx = idx + 1; sw_timeit(idx/nRun*100,0,param.tid); end @@ -261,7 +339,7 @@ x = x(sortIdx,:); % set the best fit to the spinw object -param.func(obj,x(1,:)); +param.func(obj,x(1,1:nPar)); % Store all output in a struct variable. fitsp.obj = copy(obj); @@ -281,7 +359,7 @@ dd = dat; end assert(strcmp(dat.ulabel{1}, '|Q|'), 'spinw:fitpow:invalidinput', 'Input Horace object is not a powder cut'); - out = struct('x', dd.p{1}, 'y', dd.s, 'e', dd.e, 'qmin', dd.iint(1,1), 'qmax', dd.iint(2,1)); + out = struct('x', (dd.p{1}(1:end-1) + dd.p{1}(2:end))/2, 'y', dd.s, 'e', dd.e, 'qmin', dd.iint(1,1), 'qmax', dd.iint(2,1)); end function out = verify_cut_struct(dat) @@ -350,21 +428,62 @@ [x, y, e] = deal([], [], []); is_same_e = true; for ii = 1:numel(outstruct) - x = [x; ii+([1:numel(outstruct(ii).y)]/1e9)]; - y = [y; outstruct(ii).y]; - e = [e; outstruct(ii).e]; + x = [x; ii+([1:numel(outstruct(ii).y)].'/1e9)]; %#ok + y = [y; outstruct(ii).y(:)]; + e = [e; outstruct(ii).e(:)]; if ii == 1 evec = outstruct(ii).x; else - if ~is_same_e && (numel(evec) ~= numel(outstruct(ii).x) ... + if ~is_same_e && (numel(evec) ~= numel(outstruct(ii).x) ... || sum(abs(evec(:) - outstruc(ii).x(:))) > 1e-5) is_same_e = false; end end end + e(isnan(y)) = 1; + y(isnan(y)) = 0; + e(e==0) = 1; lindat = struct('x', x, 'y', y, 'e', e, 'is_same_e', is_same_e); end +function out = is_same_data(d0, data) + out = numel(d0) == numel(data) && numel(d0(1).x) == numel(data(1).x) && ... + all(d0(1).x(:) == data(1).x(:)); +end + +function out = guess_bkg(data, bktype) + xx = []; yy = []; + for ii = 1:numel(data) + xx = [xx; data(ii).x(:)]; + yy = [yy; data(ii).y(:)]; + end + my = min(yy); + idx = find(yy < (max(yy)-my)/20+my); + [xx, isr] = sort(xx(idx)); + yy = yy(idx(isr)); + p2 = yy(end) / xx(end).^2; + p1 = (yy(end) - yy(1)) / (xx(end) - xx(1)); + p0 = min(yy); + if strcmp(bktype, 'linear') + out = fminsearch(@(p) sum(abs(yy - (p(1)*xx + p(2)))), [p1 p0]); + elseif strcmp(bktype, 'quadratic') + out = fminsearch(@(p) sum(abs(yy - (p(1)*xx.^2 + p(2)*xx + p(3)))), [p2 p1 p0]); + end +end + +function out = guess_intscale(swobj, data, param) + sfs = zeros(numel(data), 1); + for ii = 1:numel(data) + qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); + spec = swobj.powspec(qq, 'Evect', data(ii).evect, 'nRand', 100, 'fibo', param.fibo, ... + 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... + 'imagChk', param.imagChk); + spec = param.conv_func(spec); + yy = sum(spec.swConv, 2); + sfs(ii) = max(data(ii).y) / max(yy); + end + out = mean(sfs); +end function yCalc = pow_fitfun(obj, lindat, data, parfunc, x, param) % calculates the powder spin wave spectrum to directly compare to data @@ -378,7 +497,21 @@ % x Actual parameter values. % param Parameters for the spinwave function. -persistent swp hf lines; +persistent swp hf lines d0 useCache cache nCache iCache nEval; + +if ischar(obj) + if strcmp(obj, 'clear_cache') + cache = []; + nCache = []; + iCache = 1; + useCache = true; + elseif strcmp(obj, 'no_cache') + useCache = false; + end + nEval = 0; + return; +end + if isempty(swp) swp = swpref; end @@ -386,87 +519,150 @@ swp.fid = 0; cleanupObj = onCleanup(@()swpref.setpref('fid', fid0)); -parfunc(obj,x); +% Strip out background and intensity parameters before feeding the SpinW model +x_rem = x; +if param.verbose + fprintf('Evaluation #%i with parameters:', nEval); fprintf('%f ', x_rem); + nEval = nEval + 1; +end +sf = 1.0; +if param.intensity_scale < 0 + sf = abs(param.intensity_scale); +else + if ~isempty(param.bkgfun) + bkgpar = x_rem((end-param.nbkgpar+1):end); + x_rem = x_rem(1:(end-param.nbkgpar)); + end + if param.fit_scale + sf = x_rem(end); + x_rem = x_rem(1:end-1); + end +end + +% Modify the SpinW model with new parameters +parfunc(obj, x_rem); + +% Other setup +fitstruc = data; +nD = numel(data); +nY = 0; if param.plot - if isempty(hf) || ~isvalid(hf) + if isempty(hf) || ~isvalid(hf) && (isempty(d0) || is_same_data(d0, data)) hf = findobj('Tag', 'fitpow_plot'); end if isempty(hf) || ~isvalid(hf) hf = figure('Tag', 'fitpow_plot'); - for ii = 1:numel(data) - subplot(1, numel(data), ii); hold all; + nR = ceil(nD / 4); + for ii = 1:nD + subplot(nR, ceil(nD / nR), ii); hold all; errorbar(data(ii).x, data(ii).y, data(ii).e, '.k'); + title(sprintf('%4.2f < |Q| < %4.2f', data(ii).qmin, data(ii).qmax)); + box on; end lines = []; + d0 = data; end set(0, 'CurrentFigure', hf); end -fitstruc = data; -sf = param.scale_factor; -nD = numel(data); -nY = 0; - -if lindat.is_same_e - nY = nD * numel(data(1).x); - qq = zeros(nD * param.nQ); - for ii = 1:nD - i0 = (ii - 1) * param.nQ + 1; - qq(i0:(i0 + param.nQ - 1)) = linspace(data(ii).qmin, data(ii).qmax, param.nQ); - end - spec = obj.powspec(qq, 'Evect', data(1).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... - 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... - 'imagChk', param.imagChk); - spec = param.conv_func(spec); - for ii = 1:nD - i0 = (ii - 1) * param.nQ + 1; - fitstruc(ii).y = sum(spec.swConv(:, i0:(i0 + param.nQ - 1)), 2); +not_cached = true; +if ~isempty(cache) && useCache + for ii = 1:numel(cache) + if x_rem == cache{ii}{1} + not_cached = false; + fitstruc = cache{ii}{2}; + break + end end -else - for ii = 1:nD - qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); - spec = obj.powspec(qq, 'Evect', data(ii).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... +end +if not_cached + if lindat.is_same_e + nY = nD * numel(data(1).x); + qq = zeros(nD * param.nQ); + for ii = 1:nD + i0 = (ii - 1) * param.nQ + 1; + qq(i0:(i0 + param.nQ - 1)) = linspace(data(ii).qmin, data(ii).qmax, param.nQ); + end + spec = obj.powspec(qq, 'Evect', data(1).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... 'imagChk', param.imagChk); spec = param.conv_func(spec); - fitstruc(ii).y = sum(spec.swConv, 2); - nY = nY + numel(fitstruc(ii).y); + for ii = 1:nD + i0 = (ii - 1) * param.nQ + 1; + fitstruc(ii).y = sum(spec.swConv(:, i0:(i0 + param.nQ - 1)), 2); + end + else + for ii = 1:nD + qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); + spec = obj.powspec(qq, 'Evect', data(ii).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... + 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... + 'imagChk', param.imagChk); + spec = param.conv_func(spec); + fitstruc(ii).y = sum(spec.swConv, 2); + nY = nY + numel(fitstruc(ii).y); + end + end + if useCache + if isempty(nCache) % Determines cache size from free memory up to max 1GB + nCache = min([numel(x_rem)*5, ceil((min([sw_freemem, 1e9]) / 4) / (numel(fitstruc) * numel(fitstruc(1).y) * 3))]); + end + cache{iCache} = {x_rem, fitstruc}; + iCache = iCache + 1; + if iCache > nCache + iCache = 1; + end end end yCalc = zeros(nY, 1); i0 = 1; for ii = 1:numel(data) - if param.fit_separate - sf = 1.0; - ssf = fminsearch(@(x)sum(abs(data(ii).y(:) - x*fitstruc(ii).y(:)), 'omitnan'), max(data(ii).y) / max(fitstruc(ii).y)); - if ssf < 1e5 - fitstruc(ii).y = ssf * fitstruc(ii).y; - end + %if param.fit_separate + % sf = 1.0; + % ssf = fminsearch(@(x)sum(abs(data(ii).y(:) - x*fitstruc(ii).y(:)), 'omitnan'), max(data(ii).y) / max(fitstruc(ii).y)); + % if ssf < 1e5 + % fitstruc(ii).y = ssf * fitstruc(ii).y; + % end + %end + if sf ~= 1.0 + fitstruc(ii).y = fitstruc(ii).y * sf; end i1 = i0 + numel(fitstruc(ii).y) - 1; - yCalc(i0:i1) = fitstruc(ii).y; + if ~isempty(param.bkgfun) + yCalc(i0:i1) = fitstruc(ii).y + param.bkgfun(data(ii).x, bkgpar); + else + yCalc(i0:i1) = fitstruc(ii).y; + end i0 = i1 + 1; end -if param.fit_together - sf = fminsearch(@(x)sum(abs(lindat.y(:) - x*yCalc(:)), 'omitnan'), max(l0.y) / max(yCalc)); -end +%if param.fit_together +% sf = fminsearch(@(x)sum(abs(lindat.y(:) - x*yCalc(:)), 'omitnan'), max(l0.y) / max(yCalc)); +%end -if sf ~= 1.0 - yCalc = sf * yCalc; +yCalc(isnan(yCalc)) = 0; +if ~param.imagChk + yCalc = abs(yCalc); end if param.plot + nR = ceil(nD / 4); + i0 = 1; for ii = 1:numel(data) - subplot(1, numel(data), ii); + subplot(nR, ceil(nD / nR), ii); try %#ok delete(lines(ii)); end - lines(ii) = plot(fitstruc(ii).x, sf * fitstruc(ii).y, '-r'); + i1 = i0 + numel(fitstruc(ii).y) - 1; + lines(ii) = plot(fitstruc(ii).x, yCalc(i0:i1), '-r'); + i0 = i1 + 1; end drawnow; end +if param.verbose + fprintf('. chi2 = %d\n', sum((yCalc - lindat.y).^2 ./ lindat.e.^2) / numel(yCalc)); +end + end From 3189c02edb6c2cfea737876ce8175fc7d6ad706e Mon Sep 17 00:00:00 2001 From: Duc Le Date: Sun, 21 Jan 2024 22:18:25 +0000 Subject: [PATCH 03/43] Add mex file to calculate q-convolution loop --- external/sw_qconv/sw_qconv.cpp | 85 ++++++++++++++++++++++++++++++++++ swfiles/sw_instrument.m | 21 +++++---- swfiles/sw_mex.m | 3 ++ 3 files changed, 100 insertions(+), 9 deletions(-) create mode 100644 external/sw_qconv/sw_qconv.cpp diff --git a/external/sw_qconv/sw_qconv.cpp b/external/sw_qconv/sw_qconv.cpp new file mode 100644 index 000000000..7d2e3e91b --- /dev/null +++ b/external/sw_qconv/sw_qconv.cpp @@ -0,0 +1,85 @@ +#include "mex.h" +#include +#include +#include +#include + +template +void loop(T *swOut, const mwSize *d1, const double stdG, const double *Qconv, const T *swConv, size_t i0, size_t i1) { + for (size_t ii=i0; ii fG(d1[1], 0.0); + for (size_t jj=0; jj +void do_calc(T *swOut, const mwSize *d1, const double stdG, const double *Qconv, const T *swConv, size_t nThread) { + if (d1[1] > 10*nThread) { + size_t nBlock = d1[1] / nThread; + size_t i0 = 0, i1 = nBlock; + std::vector swv(nThread); + std::vector threads; + for (size_t ii=0; ii, std::ref(swv[ii]), std::ref(d1), stdG, std::ref(Qconv), std::ref(swConv), i0, i1) + ); + i0 = i1; + i1 += nBlock; + if (i1 > d1[1] || ii == (nThread - 2)) { + i1 = d1[1]; } + } + for (size_t ii=0; ii(swOut, d1, stdG, Qconv, swConv, 0, d1[1]); + } +} + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +{ + if (nrhs < 3) { + throw std::runtime_error("sw_qconv: Requires 3 arguments"); + } + size_t nThreads = 8; + if (nrhs == 4) { + nThreads = (size_t)(*mxGetDoubles(prhs[3])); + } + if (mxIsComplex(prhs[1])) { throw std::runtime_error("Arg 2 is complex\n"); } + if (mxIsComplex(prhs[2])) { throw std::runtime_error("Arg 3 is complex\n"); } + const mwSize *d1 = mxGetDimensions(prhs[0]); + const mwSize *d2 = mxGetDimensions(prhs[1]); + if (d1[1] != d2[1]) { throw std::runtime_error("Arg 1 and 2 size mismatch\n"); } + if (mxGetNumberOfElements(prhs[2]) > 1) { throw std::runtime_error("Arg 3 should be scalar\n"); } + double *Qconv = mxGetDoubles(prhs[1]); + double stdG = *(mxGetDoubles(prhs[2])); + if (mxIsComplex(prhs[0])) { + std::complex *swConv = reinterpret_cast*>(mxGetComplexDoubles(prhs[0])); + plhs[0] = mxCreateDoubleMatrix(d1[0], d1[1], mxCOMPLEX); + std::complex *swOut = reinterpret_cast*>(mxGetComplexDoubles(plhs[0])); + do_calc(swOut, d1, stdG, Qconv, swConv, nThreads); + } else { + double *swConv = mxGetDoubles(prhs[0]); + plhs[0] = mxCreateDoubleMatrix(d1[0], d1[1], mxREAL); + double *swOut = mxGetDoubles(plhs[0]); + do_calc(swOut, d1, stdG, Qconv, swConv, nThreads); + } +} diff --git a/swfiles/sw_instrument.m b/swfiles/sw_instrument.m index c51054c4e..a8aa872e0 100644 --- a/swfiles/sw_instrument.m +++ b/swfiles/sw_instrument.m @@ -229,16 +229,19 @@ stdG = param.dQ/2.35482; for jj = 1:nPlot - swConv = spectra.swConv{jj}; - swConvTemp = swConv * 0; - for ii = 1:numel(Qconv) - % Gaussian with intensity normalised to 1, centered on E(ii) - fG = exp(-((Qconv-Qconv(ii))/stdG).^2/2); - fG = fG/sum(fG); - swConvTemp = swConvTemp + swConv(:,ii) * fG; - + if pref.usemex + spectra.swConv{jj} = sw_qconv(spectra.swConv{jj}, Qconv, stdG, pref.nthread); + else + swConv = spectra.swConv{jj}; + swConvTemp = swConv * 0; + for ii = 1:numel(Qconv) + % Gaussian with intensity normalised to 1, centered on E(ii) + fG = exp(-((Qconv-Qconv(ii))/stdG).^2/2); + fG = fG/sum(fG); + swConvTemp = swConvTemp + swConv(:,ii) * fG; + end + spectra.swConv{jj} = swConvTemp; end - spectra.swConv{jj} = swConvTemp; end fprintf0(fid,'Finite instrumental momentum resolution of %5.3f A-1 is applied.\n',param.dQ); diff --git a/swfiles/sw_mex.m b/swfiles/sw_mex.m index 172ffbf41..b51f28384 100755 --- a/swfiles/sw_mex.m +++ b/swfiles/sw_mex.m @@ -45,6 +45,7 @@ function sw_mex(varargin) chol_omp_dir = [sw_rootdir filesep 'external' filesep 'chol_omp']; mtimesx_dir = [sw_rootdir filesep 'external' filesep 'mtimesx']; swloop_dir = [sw_rootdir filesep 'external' filesep 'swloop']; + swqconv_dir = [sw_rootdir filesep 'external' filesep 'sw_qconv']; eigen_ver = '3.4.0'; if ~exist([swloop_dir filesep 'eigen-' eigen_ver], 'dir') cd(swloop_dir); @@ -97,6 +98,8 @@ function sw_mex(varargin) cd(swloop_dir); mex('-v','-R2018a',['-Ieigen-' eigen_ver],'swloop.cpp') end + cd(swqconv_dir); + mex('-R2018a','sw_qconv.cpp') % return back to original folder cd(aDir); end From e2f4abe2c66ccd4d416dcdaace9a4018ce15e5fb Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 7 Mar 2024 09:21:40 +0000 Subject: [PATCH 04/43] Fix variables not assigned bug when maxiter=0 in lm3 --- swfiles/+ndbase/lm3.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/swfiles/+ndbase/lm3.m b/swfiles/+ndbase/lm3.m index 9cd3a8abc..43f682acb 100644 --- a/swfiles/+ndbase/lm3.m +++ b/swfiles/+ndbase/lm3.m @@ -256,7 +256,7 @@ % function values at start yBest = yCalc(:); - +converged=false; if Niter > 0 % optimisation @@ -357,6 +357,8 @@ tmp = repmat(1./sqrt(diag(cov)),[1,NpFree]); cor = tmp.*cov.*tmp'; else + sigP = [] + iter=0 chisqr = chi2Best/nnorm; ok = true; warning('WARNING: Convergence not achieved') From 29d9c013543619f528140bbba83d17bb9129b015 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 7 Mar 2024 09:23:14 +0000 Subject: [PATCH 05/43] Add new class to sw_files --- swfiles/sw_fitpowder.m | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 swfiles/sw_fitpowder.m diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m new file mode 100644 index 000000000..6d32e1eb4 --- /dev/null +++ b/swfiles/sw_fitpowder.m @@ -0,0 +1,24 @@ +classdef sw_fitpowder < handle & matlab.mixin.SetGet + % class to fit powders + + % + % [spinw.copy], [spinw.struct], [Comparing handle and value classes](https://www.mathworks.com/help/matlab/matlab_oop/comparing-handle-and-value-classes.html) + % + + properties (SetObservable) + swobj + data + end + + methods + function obj = sw_fitpowder(swobj, data) + % constructor + obj.swobj = swobj; + obj.data = data; + end + + function result = fit(obj, fit_struct) + result = obj.swobj.fitpow(obj.data, fit_struct); + end + end +end \ No newline at end of file From bb876c1b04ec25ad1d1624aaff37f9645ff1b650 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 19 Mar 2024 13:42:03 +0000 Subject: [PATCH 06/43] Add dE option of sw_egrid as input argument to powspec --- swfiles/@spinw/powspec.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index 97b8659d7..fb00771c9 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -213,9 +213,9 @@ inpForm.defval = [inpForm.defval {false false 0 'ebin' 'Sperp' }]; inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 -4] [1 -5] }]; -inpForm.fname = [inpForm.fname {'fid'}]; -inpForm.defval = [inpForm.defval {-1 }]; -inpForm.size = [inpForm.size {[1 1]}]; +inpForm.fname = [inpForm.fname {'fid' , 'dE'}]; +inpForm.defval = [inpForm.defval {-1, []}]; +inpForm.size = [inpForm.size {[1 1], [-6, -7]}]; param = sw_readparam(inpForm, varargin{:}); @@ -322,7 +322,7 @@ specQ.obj = obj; % use edge grid by default specQ = sw_egrid(specQ,struct('Evect',param.Evect,'T',param.T,'binType',param.binType,... - 'imagChk',param.imagChk,'component',param.component),'noCheck'); + 'imagChk',param.imagChk,'component',param.component, 'dE', param.dE),'noCheck'); powSpec = squeeze(sum(reshape(specQ.swConv, max(1, nE), nk, nQ), 2) / param.nRand); fprintf0(fid,'Calculation finished.\n'); From 47848ac8dc81b7a204d6fc7124c0ed6c167bd9c8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 19 Mar 2024 13:42:50 +0000 Subject: [PATCH 07/43] Class interface to fit 2D and 1D powder spectra Currenlty uses lsqnonlin to do the minmisation - this is in the optimization toolbox so not availible to all. The plan is to replace this with a call to an optimiser of choice (including some packaged with spinw) --- swfiles/sw_fitpowder.m | 245 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 237 insertions(+), 8 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 6d32e1eb4..05545e0ce 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -1,24 +1,253 @@ classdef sw_fitpowder < handle & matlab.mixin.SetGet % class to fit powders - % % [spinw.copy], [spinw.struct], [Comparing handle and value classes](https://www.mathworks.com/help/matlab/matlab_oop/comparing-handle-and-value-classes.html) % properties (SetObservable) - swobj - data + swobj; + y + e + ebin_cens + modQ_cens + nQ = 10 % number |Q| bins to calc. in integration limits of 1D cut + ndim + ncuts = 0; + sw_instrument_args = struct() + powspec_args = struct('nRand', 1e3, 'T', 0, 'formfact', true, 'hermit', false, 'imagChk', false, 'fibo', true, 'binType', 'cbin', 'component', 'Sperp') + fit_func + cost_function = "Rsq"; % "Rsq" or "chisq" + background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) + fbg = @(en, modQ, slope_en, slope_modQ, intercept) slope_en*en(:) + slope_modQ*modQ(:)' + intercept + % parameters + nparams_model + nparams_bg + params + bounds end - methods - function obj = sw_fitpowder(swobj, data) + methods + function obj = sw_fitpowder(swobj, data, fit_func, model_params, varargin) % constructor obj.swobj = swobj; - obj.data = data; + obj.fit_func = fit_func; + bg_params = []; + scale = 1; + % set background startegy and func + for ikey = 1:2:numel(varargin) + switch varargin{ikey} + case 'backgroundStrategy' + obj.background_strategy = varargin{ikey+1}; + case 'initialBackgroundParameters' + bg_params = varargin{ikey+1}; + case 'scale' + scale = varargin{ikey+1}; + end + end + % initialise parameters + obj.add_data(data) + obj.initialise_parameters_and_bounds(model_params, bg_params, scale) end - function result = fit(obj, fit_struct) - result = obj.swobj.fitpow(obj.data, fit_struct); + function initialise_parameters_and_bounds(obj, model_params, bg_params, scale) + obj.nparams_model = numel(model_params); + obj.nparams_bg = obj.get_nparams_in_background_func(); + if isempty(bg_params) + % zero intialise + if obj.background_strategy == "independent" + nparams_bg_total = obj.ncuts * obj.nparams_bg; + else + nparams_bg_total = obj.nparams_bg; + end + bg_params = zeros(nparams_bg_total, 1); + end + obj.params = [model_params(:); bg_params(:); scale]; + obj.bounds = [-inf, inf].*ones(numel(obj.params),1); + obj.bounds(end,1) = 0; % set lower bound of scale to be 0 + end + + function fix_parameter(obj, iparams, varargin) + for ii = 1:numel(iparams) + iparam = obj.get_parameter_index(iparams(ii), varargin{:}); + obj.bounds(iparam, :) = obj.params(iparam); + end + end + + function set_parameter(obj, iparams, values, varargin) + for ii = 1:numel(iparams) + iparam = obj.get_parameter_index(iparams(ii), varargin{:}); + obj.params(iparam) = values(ii); + end + end + + function set_parameter_bounds(obj, iparams, varargin) + lb = []; + ub = []; + for ikey = 1:2:numel(varargin) + switch varargin{ikey} + case 'lb' + lb = varargin{ikey+1}; + case 'ub' + ub = varargin{ikey+1}; + end + end + for ii = 1:numel(iparams) + iparam = obj.get_parameter_index(iparams(ii), varargin{:}); + if ~isempty(lb) + obj.bounds(iparam, 1) = lb(ii); + end + if ~isempty(ub) + obj.bounds(iparam, 2) = ub(ii); + end + end + end + + function add_data(obj, data) + if numel(data) == 1 && numel(data.x) == 2 + % 2D + obj.ndim = 2; + obj.y = data.y'; % nE x n|Q| + obj.e = data.e'; + obj.ebin_cens = data.x{1}; + obj.modQ_cens = data.x{2}; + else + % 1D + obj.ndim = 1; + obj.ebin_cens = data(1).x; + obj.ncuts = numel(data); + for cut = data + obj.y = [obj.y cut.y(:)]; + obj.e = [obj.e cut.e(:)]; + obj.modQ_cens = [obj.modQ_cens, linspace(cut.qmin, cut.qmax, obj.nQ)]; % does this play nicely with sw_instrument Q convolution? + end + end + obj.powspec_args.Evect = obj.ebin_cens(:)'; % row vect + end + + function crop_energy_range(obj, emin, emax) + % crop data + ikeep = obj.ebin_cens >= emin & obj.ebin_cens <= emax; + obj.ebin_cens = obj.ebin_cens(ikeep); + obj.y = obj.y(ikeep, :); + obj.e = obj.e(ikeep, :); + end + + function bg = calc_background(obj, params) + % add background + istart = obj.nparams_model + 1; + if obj.background_strategy == "planar" + bg_params = num2cell(params(istart:istart+obj.nparams_bg-1)); + 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)); + for icut = 1:obj.ncuts + iend = istart + +obj.nparams_bg - 1; + bg_params = num2cell(params(istart:iend)); + bg(:,icut) = obj.fbg(obj.ebin_cens(:), bg_params{:}); + istart = iend + 1; + end + end + end + + function [ycalc, bg] = calc_spinwave_spec(obj, params) + % set spinw model interactions + obj.fit_func(obj.swobj, params(1:obj.nparams_model)'); + % simulate powder spectrum + spec = obj.swobj.powspec(obj.modQ_cens, obj.powspec_args); + spec = sw_instrument(spec, obj.sw_instrument_args); + % scale + ycalc = params(end)*spec.swConv; + % add background + bg = calc_background(obj, params); + ycalc = ycalc + bg; + if obj.ndim == 1 + % integrate nQ |Q| points for each cut + ycalc = obj.rebin_powspec_to_1D_cuts(ycalc); + end + end + + function resid_sq = calc_resid_sq(obj, params) + % evaluate fit function + [ycalc, ~] = calc_spinwave_spec(obj, params); + resid_sq = (obj.y - ycalc).^2; + if obj.cost_function == "chisq" + resid_sq = resid_sq./(obj.e.^2); + end + % exclude nans in both ycalc and input data + resid_sq = resid_sq(~isnan(resid_sq)); + end + + function result = fit(obj, options) + % call fit + [p, cost, resid_sq, exitflag, ~,~,jacobian] = lsqnonlin(@obj.calc_resid_sq,obj.params,obj.bounds(:,1),obj.bounds(:,2), options); + result.p = p; + result.exitflag = exitflag; + % calculate errors + ndof = numel(resid_sq) - numel(p); + result.reduced_cost = (cost/ndof); + cov = result.reduced_cost * inv(jacobian' * jacobian); + result.perr = diag(cov); + end + + function estimate_scale_factor(obj) + % set scale factor to 1 + params = obj.params; + params(end) = 1; + [ycalc, bg] = calc_spinwave_spec(obj, params); + if obj.ndim == 1 + % integrate nQ |Q| points for each cut + bg = obj.rebin_powspec_to_1D_cuts(bg); + end + scale = max(obj.y - bg, [], "all")/max(ycalc, [], "all"); + obj.params(end) = scale; + end + + end + + % private + methods (Hidden=true, Access = private) + 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(sum(ycalc, 2, 'omitnan')); + 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 + else + nbg_pars = nargin(obj.fbg) - 1; % dependent variable energy + end + end + + function iparam = get_parameter_index(obj, iparam, varargin) + icut = 0; + is_background_param = false; + for ikey = 1:2:numel(varargin) + switch varargin{ikey} + case 'background' + is_background_param = varargin{ikey+1}; + case 'cutIndex' + if obj.ndim == 1 && obj.background_strategy=="independent" + icut = varargin{ikey+1}; + if icut < 1 || icut > obj.ncuts + error("Invalid CutIndex supplied.") + end + else + error("Can only supply cutIndex when fitting 1D cuts with independent backgrounds") + end + end + end + % find index in vector of all parameters + if is_background_param + iparam = iparam + obj.nparams_model; + if obj.background_strategy == "independent" + nbg_pars = obj.get_nparams_in_background_func(); + iparam = iparam + (icut-1)*nbg_pars; + end + end end end end \ No newline at end of file From 1410321b456373228aad1aca539f2b0929c66bd0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Mar 2024 10:55:49 +0000 Subject: [PATCH 08/43] Change to support user supplied minimizer (assumed from ndbase) --- swfiles/sw_fitpowder.m | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 05545e0ce..db9138382 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -5,6 +5,7 @@ % properties (SetObservable) + % data swobj; y e @@ -13,13 +14,16 @@ nQ = 10 % number |Q| bins to calc. in integration limits of 1D cut ndim ncuts = 0; + % spinw funtion inputs sw_instrument_args = struct() powspec_args = struct('nRand', 1e3, 'T', 0, 'formfact', true, 'hermit', false, 'imagChk', false, 'fibo', true, 'binType', 'cbin', 'component', 'Sperp') + % functions fit_func cost_function = "Rsq"; % "Rsq" or "chisq" background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) fbg = @(en, modQ, slope_en, slope_modQ, intercept) slope_en*en(:) + slope_modQ*modQ(:)' + intercept - % parameters + % fit and parameters + optimizer = @ndbase.simplex nparams_model nparams_bg params @@ -167,27 +171,24 @@ function crop_energy_range(obj, emin, emax) end end - function resid_sq = calc_resid_sq(obj, params) + function resid_sq_sum = calc_cost_func(obj, params) % evaluate fit function [ycalc, ~] = calc_spinwave_spec(obj, params); - resid_sq = (obj.y - ycalc).^2; + resid = (obj.y - ycalc); if obj.cost_function == "chisq" - resid_sq = resid_sq./(obj.e.^2); + resid = resid./(obj.e); end % exclude nans in both ycalc and input data - resid_sq = resid_sq(~isnan(resid_sq)); + ikeep = ~isnan(resid); + resid_sq_sum = resid(ikeep)'*resid(ikeep); end - function result = fit(obj, options) - % call fit - [p, cost, resid_sq, exitflag, ~,~,jacobian] = lsqnonlin(@obj.calc_resid_sq,obj.params,obj.bounds(:,1),obj.bounds(:,2), options); - result.p = p; - result.exitflag = exitflag; - % calculate errors - ndof = numel(resid_sq) - numel(p); - result.reduced_cost = (cost/ndof); - cov = result.reduced_cost * inv(jacobian' * jacobian); - result.perr = diag(cov); + function result = fit(obj, varargin) + % setup cell for output of ndbase optimizer/minimizer + result = cell(1,nargout(obj.optimizer)); + [result{:}] = obj.optimizer([], @obj.calc_cost_func, obj.params, ... + 'lb', obj.bounds(:,1), ... + 'ub', obj.bounds(:,2)); end function estimate_scale_factor(obj) From b94e944b5f3f2c3cb4981580adec9fbb9caa11b0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Mar 2024 13:50:57 +0000 Subject: [PATCH 09/43] Add default fbg for 1D independent cuts and fix bug in bg subtraction --- swfiles/sw_fitpowder.m | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index db9138382..58e70fe93 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -21,7 +21,9 @@ fit_func cost_function = "Rsq"; % "Rsq" or "chisq" background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) - fbg = @(en, modQ, slope_en, slope_modQ, intercept) slope_en*en(:) + slope_modQ*modQ(:)' + intercept + fbg + 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; % fit and parameters optimizer = @ndbase.simplex nparams_model @@ -48,8 +50,12 @@ scale = varargin{ikey+1}; end end - % initialise parameters obj.add_data(data) + if obj.background_strategy == "planar" + obj.fbg = obj.fbg_planar; + else + obj.fbg = obj.fbg_indep; + end obj.initialise_parameters_and_bounds(model_params, bg_params, scale) end @@ -162,13 +168,21 @@ function crop_energy_range(obj, emin, emax) spec = sw_instrument(spec, obj.sw_instrument_args); % scale ycalc = params(end)*spec.swConv; - % add background + % calc background bg = calc_background(obj, params); - ycalc = ycalc + bg; - if obj.ndim == 1 + if obj.background_strategy == "planar" + ycalc = ycalc + bg; % add planar bg before any rebinning + if obj.ndim == 1 + % integrate nQ |Q| points for each cut + ycalc = obj.rebin_powspec_to_1D_cuts(ycalc); + end + else % integrate nQ |Q| points for each cut ycalc = obj.rebin_powspec_to_1D_cuts(ycalc); + % add background to individual cuts after rebin + ycalc = ycalc + bg; end + end function resid_sq_sum = calc_cost_func(obj, params) @@ -196,10 +210,6 @@ function estimate_scale_factor(obj) params = obj.params; params(end) = 1; [ycalc, bg] = calc_spinwave_spec(obj, params); - if obj.ndim == 1 - % integrate nQ |Q| points for each cut - bg = obj.rebin_powspec_to_1D_cuts(bg); - end scale = max(obj.y - bg, [], "all")/max(ycalc, [], "all"); obj.params(end) = scale; end From 6a4de70f1034f16a4e8fbce5be84f02c24988ef4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 20 Mar 2024 17:33:06 +0000 Subject: [PATCH 10/43] Fix bug in en bins passed to powspec when data cropped --- swfiles/sw_fitpowder.m | 1 + 1 file changed, 1 insertion(+) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 58e70fe93..a5381ebed 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -138,6 +138,7 @@ function crop_energy_range(obj, emin, emax) % crop data ikeep = obj.ebin_cens >= emin & obj.ebin_cens <= emax; obj.ebin_cens = obj.ebin_cens(ikeep); + obj.powspec_args.Evect = obj.ebin_cens(:)'; obj.y = obj.y(ikeep, :); obj.e = obj.e(ikeep, :); end From 90ee9cfdab579826246d004bd03062f2f18a56c0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 18 Apr 2024 09:50:22 +0100 Subject: [PATCH 11/43] Fix bug with estimating scale with planar bg and 1D cuts --- swfiles/sw_fitpowder.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index a5381ebed..05e3975af 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -211,7 +211,7 @@ function estimate_scale_factor(obj) params = obj.params; params(end) = 1; [ycalc, bg] = calc_spinwave_spec(obj, params); - scale = max(obj.y - bg, [], "all")/max(ycalc, [], "all"); + scale = max(obj.y - mean(bg(:)), [], "all")/max(ycalc, [], "all"); obj.params(end) = scale; end From 094eca9c81a2bcacc49a48a08892a2591ac9a27b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 18 Apr 2024 09:51:02 +0100 Subject: [PATCH 12/43] Pass varargin to optimizer --- swfiles/sw_fitpowder.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 05e3975af..bc76259a5 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -203,7 +203,8 @@ function crop_energy_range(obj, emin, emax) result = cell(1,nargout(obj.optimizer)); [result{:}] = obj.optimizer([], @obj.calc_cost_func, obj.params, ... 'lb', obj.bounds(:,1), ... - 'ub', obj.bounds(:,2)); + 'ub', obj.bounds(:,2), ... + varargin{:}); end function estimate_scale_factor(obj) From 738038e7e7119b60033ad214bfc71507f74ba4cc Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 18 Apr 2024 10:40:52 +0100 Subject: [PATCH 13/43] Reshape parameters and bounds to work with all ndbase optmisers And pass model parameters as a row vector as requried by matparser --- swfiles/sw_fitpowder.m | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index bc76259a5..290d7280b 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -163,7 +163,8 @@ function crop_energy_range(obj, emin, emax) function [ycalc, bg] = calc_spinwave_spec(obj, params) % set spinw model interactions - obj.fit_func(obj.swobj, params(1:obj.nparams_model)'); + % pass params as row-vector as expected by matparser + obj.fit_func(obj.swobj, reshape(params(1:obj.nparams_model), 1, [])); % simulate powder spectrum spec = obj.swobj.powspec(obj.modQ_cens, obj.powspec_args); spec = sw_instrument(spec, obj.sw_instrument_args); @@ -201,9 +202,10 @@ function crop_energy_range(obj, emin, emax) function result = fit(obj, varargin) % setup cell for output of ndbase optimizer/minimizer result = cell(1,nargout(obj.optimizer)); - [result{:}] = obj.optimizer([], @obj.calc_cost_func, obj.params, ... - 'lb', obj.bounds(:,1), ... - 'ub', obj.bounds(:,2), ... + % pass params and bounds as rows for ndbase optimisers + [result{:}] = obj.optimizer([], @obj.calc_cost_func, obj.params(:)', ... + 'lb', obj.bounds(:,1)', ... + 'ub', obj.bounds(:,2)', ... varargin{:}); end From a769f7c1ee1c28d12e0545857b0d28466cd05730 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 18 Apr 2024 16:46:43 +0100 Subject: [PATCH 14/43] Add methods to plot results --- swfiles/sw_fitpowder.m | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 290d7280b..166e3cac7 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -218,6 +218,46 @@ function estimate_scale_factor(obj) obj.params(end) = scale; end + function plot_result(obj, params, varargin) + [ycalc, ~] = obj.calc_spinwave_spec(params); + if obj.ndim == 1 + obj.plot_1d_cuts_on_data(ycalc, varargin{:}) + else + obj.plot_2d_contour_on_data(ycalc, varargin{:}) + end + end + + function plot_1d_cuts_on_data(obj, ycalc, varargin) + figure("color","white"); + modQs = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); + for icut = 1:obj.ncuts + ax = subplot(1, obj.ncuts, icut); + hold on; box on; + plot(ax, obj.ebin_cens, obj.y(:,icut), 'ok'); + plot(ax, obj.ebin_cens, ycalc(:,icut), '-r', varargin{:}); + xlim(ax, [0, 12]); + ylim(ax, [-2, 30]); + xlabel(ax, 'Energy (meV)') + ylabel(ax, 'Intensity'); + title(ax, num2str(modQs(icut), 2) + " $\AA^{-1}$", 'interpreter','latex') + end + end + + function plot_2d_contour_on_data(obj, ycalc, varargin) + % varargin passed to contour + figure("color","white"); + ax = subplot(1,1,1); + box on; hold on; + h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, obj.y); + h.AlphaData = obj.y > 0; % make empty bins transparent + contour(ax, obj.modQ_cens, obj.ebin_cens, ycalc, varargin{:}); + cbar = colorbar(ax); + cbar.Label.String = "Intensity"; + xlabel(ax, "$\left|Q\right| (\AA^{-1})$", 'interpreter','latex'); + ylabel(ax, "Energy (meV)"); + legend('Calculated'); + end + end % private From 13117e344a7ed737d635cd98a638574879439b63 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 26 Apr 2024 15:46:14 +0100 Subject: [PATCH 15/43] Implement caching of spinwave calc if params unchanged --- swfiles/sw_fitpowder.m | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 166e3cac7..672741932 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -30,6 +30,10 @@ nparams_bg params bounds + % cache calculation + do_cache = true + ycalc_cached = [] + model_params_cached = [] end methods @@ -161,15 +165,38 @@ function crop_energy_range(obj, emin, emax) end end + function set_caching(obj, do_cache) + obj.do_cache = do_cache; + if ~obj.do_cache + obj.clear_cache() + end + end + + function clear_cache(obj) + obj.ycalc_cached = []; + obj.model_params_cached = []; + end + function [ycalc, bg] = calc_spinwave_spec(obj, params) - % set spinw model interactions - % pass params as row-vector as expected by matparser - obj.fit_func(obj.swobj, reshape(params(1:obj.nparams_model), 1, [])); - % simulate powder spectrum - spec = obj.swobj.powspec(obj.modQ_cens, obj.powspec_args); - spec = sw_instrument(spec, obj.sw_instrument_args); + model_params = reshape(params(1:obj.nparams_model), 1, []); + if obj.do_cache && ~isempty(obj.model_params_cached) && all(abs(model_params - obj.model_params_cached) < 1e-10) + ycalc = obj.ycalc_cached; % do not repeat calc + else + % set spinw model interactions + % pass params as row-vector as expected by matparser + obj.fit_func(obj.swobj, model_params); + % simulate powder spectrum + spec = obj.swobj.powspec(obj.modQ_cens, obj.powspec_args); + spec = sw_instrument(spec, obj.sw_instrument_args); + ycalc = spec.swConv; + % cache if required + if obj.do_cache + obj.ycalc_cached = spec.swConv; + obj.model_params_cached = model_params; + end + end % scale - ycalc = params(end)*spec.swConv; + ycalc = params(end)*ycalc; % calc background bg = calc_background(obj, params); if obj.background_strategy == "planar" From 8e86112eb3118c49345cb6f47e2f3f77320f7524 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 26 Apr 2024 19:34:04 +0100 Subject: [PATCH 16/43] Add seperate setters for model, bg and scale parameters and bounds COuld possibly be made more concise, but at least is more user friendly and can adjust parameters accross multiple slices with one call now. --- swfiles/sw_fitpowder.m | 128 ++++++++++++++++++++++++++--------------- 1 file changed, 83 insertions(+), 45 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 672741932..d7c0c1b75 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -80,42 +80,98 @@ function initialise_parameters_and_bounds(obj, model_params, bg_params, scale) obj.bounds(end,1) = 0; % set lower bound of scale to be 0 end - function fix_parameter(obj, iparams, varargin) - for ii = 1:numel(iparams) - iparam = obj.get_parameter_index(iparams(ii), varargin{:}); - obj.bounds(iparam, :) = obj.params(iparam); + function fix_model_parameters(obj, iparams) + if any(iparams > obj.nparams_model) + error('sw_fitpowder', 'Parameter indices supplied must be within number of model parameters'); + end + for iparam = iparams + obj.bounds(iparam, :) = obj.params(iparam); end end - function set_parameter(obj, iparams, values, varargin) - for ii = 1:numel(iparams) - iparam = obj.get_parameter_index(iparams(ii), varargin{:}); - obj.params(iparam) = values(ii); + 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 + for iparam = obj.get_index_of_background_parameters(iparams_bg, icuts) + obj.bounds(iparam, :) = obj.params(iparam); end end - function set_parameter_bounds(obj, iparams, varargin) - lb = []; - ub = []; - for ikey = 1:2:numel(varargin) - switch varargin{ikey} - case 'lb' - lb = varargin{ikey+1}; - case 'ub' - ub = varargin{ikey+1}; + function fix_scale(obj) + obj.bounds(end, :) = obj.params(end); + end + + function set_model_parameters(obj, iparams, values) + if any(iparams > obj.nparams_model) + error('sw_fitpowder', 'Parameter indices supplied must be within number of model parameters'); + end + for ival = 1:numel(values) + obj.params(iparams(ival)) = values(ival); + end + 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 + for ival = 1:numel(values) + iparams = obj.get_index_of_background_parameters(iparams_bg(ival), icuts); + obj.params(iparams) = values(ival); + end + end + + function set_scale(obj, scale) + obj.params(end) = scale; + end + + function set_model_parameter_bounds(obj, iparams, lb, ub) + if any(iparams > obj.nparams_model) + error('sw_fitpowder', 'Parameter indices supplied must be within number of model parameters'); + end + for ibnd = 1:numel(iparams) + if ~isempty(lb) + obj.bounds(iparams(ibnd), 1) = lb(ibnd); + end + if ~isempty(ub) + obj.bounds(iparams(ibnd), 2) = ub(ibnd); end end - for ii = 1:numel(iparams) - iparam = obj.get_parameter_index(iparams(ii), varargin{:}); + 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 + for ibnd = 1:numel(iparams_bg) + iparams = obj.get_index_of_background_parameters(iparams_bg(ibnd), icuts); if ~isempty(lb) - obj.bounds(iparam, 1) = lb(ii); + obj.bounds(iparams, 1) = lb(ibnd); end if ~isempty(ub) - obj.bounds(iparam, 2) = ub(ii); + obj.bounds(iparams, 2) = ub(ibnd); end end end + function set_scale_bounds(obj, lb, ub) + if ~isempty(lb) + obj.bounds(end, 1) = lb; + end + if ~isempty(ub) + obj.bounds(end, 2) = ub; + end + end + function add_data(obj, data) if numel(data) == 1 && numel(data.x) == 2 % 2D @@ -284,7 +340,6 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) ylabel(ax, "Energy (meV)"); legend('Calculated'); end - end % private @@ -304,31 +359,14 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) end end - function iparam = get_parameter_index(obj, iparam, varargin) - icut = 0; - is_background_param = false; - for ikey = 1:2:numel(varargin) - switch varargin{ikey} - case 'background' - is_background_param = varargin{ikey+1}; - case 'cutIndex' - if obj.ndim == 1 && obj.background_strategy=="independent" - icut = varargin{ikey+1}; - if icut < 1 || icut > obj.ncuts - error("Invalid CutIndex supplied.") - end - else - error("Can only supply cutIndex when fitting 1D cuts with independent backgrounds") - end - end + function iparams = get_index_of_background_parameters(obj, iparams_bg, icuts) + if any(icuts == 0) && obj.background_strategy == "independent" + icuts = 1:obj.ncuts; % apply to all cuts end % find index in vector of all parameters - if is_background_param - iparam = iparam + obj.nparams_model; - if obj.background_strategy == "independent" - nbg_pars = obj.get_nparams_in_background_func(); - iparam = iparam + (icut-1)*nbg_pars; - end + iparams = []; + for icut = icuts + iparams = [iparams, iparams_bg + obj.nparams_model + (icut-1)*obj.nparams_bg]; end end end From 6622e09067d3385ab46c3f86e2f1107504a163f8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Sat, 27 Apr 2024 16:04:22 +0100 Subject: [PATCH 17/43] Add docs string to class --- swfiles/sw_fitpowder.m | 128 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 124 insertions(+), 4 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index d7c0c1b75..841b89c86 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -1,8 +1,128 @@ classdef sw_fitpowder < handle & matlab.mixin.SetGet - % class to fit powders - % - % [spinw.copy], [spinw.struct], [Comparing handle and value classes](https://www.mathworks.com/help/matlab/matlab_oop/comparing-handle-and-value-classes.html) - % +% Class to fit powder averaged spectra to inelastic neutron scattering data +% +% ### Syntax +% +% `fitpow = sw_fitpowder(swobj, data, fit_func, model_params, Name,Value)` +% +% ### Description +% +% Fits powder averaged spectra to constant-|Q| cuts or 2D |Q| vs en slices +% of inelastic neutron scattering data accounting for the instrument +% resolution +% +% ### Input Arguments +% +% `swobj` +% : spinwave object with magnetic structure defined +% +% `data` +% : Possible inputs depend on dimensionality of the data: +% * 2D data in either a IX_dataset_2d or struct containing fields `x`, `y` +% and `e`. The size of `y` and `e` should be (N|Q|, NEnergy) and `x` is +% a cell array of size (1,2) where x{1} contains a vector of energy +% transers and x{2} is vector of |Q| bin centers. +% * Vector of 1d datasets at constant |Q| - either IX_dataset_1d or +% structs with fields `x`, `y`, `e`, `qmin` and `qmax`, The field `x` +% is a vector of energy transfer bin centers, `qmin` and `qmax` are the +% limits of the integration in |Q|. +% +% `fit_func` +% : Function handle changing the interactions in the spinwave model. +% For example if the spinw model had matrices `J1`, `J2` and `D` then a +% `fit_func` could be +% ``` +% fit_func = @(obj, p) matparser(obj, 'param', p, 'mat', +% {'J1', 'J2', 'D(3,3)'}, 'init', true); +% ``` +% +% `model_params` +% : Vector of initial paramers to pass to fit_func +% +% ### Name-Value Pair Arguments +% +% `'backgroundStrategy'` +% : A string determining the type of background: +% * `planar` (2D planar background in |Q| and energy transfer) +% * `independent` (1D linear backgroudn as function of energy transfer) +% +% `'initialBackgroundParameters'` +% : vector containing parameters for the background (size depends on the +% dimensionality of the data and the background strategy) +% +% `'scale'` +% : Initial guess for multiplicative scale factor applied to spinwave calc +% +% ### Output Arguments +% +% `'result'` +% : cell array containgin output of sw_fitpowder.optimizer +% For ndbase optimizers in the spinw package then the reuslt can be +% unpacked as follows +% ``` +% [fitted_params, cost_val, stat] = result{:} +% ``` +% See docs for e.g. ndbase.simplex (default) for details +% +% ### Examples +% +% ``` +% >> % 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); +% >> +% >> % define resolution (from PyChop) +% >> Ei = 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 ); +% >> +% >> % fit powder +% >> fitpow = sw_fitpowder(mnf2, mnf2dat, fit_func, [J1 J2 D], 'backgroundStrategy', 'planar'); +% >> fitpow.crop_energy_range(2.0, 8.0); +% >> fitpow.powspec_args.dE = eres; +% >> fitpow.sw_instrument_args = struct('dQ', dQ, 'ThetaMin', 3.5, 'Ei', Ei); +% >> 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_parameter_bounds(3, 'lb', 0.0, 'background', true) % Set lb of constant bg = 0 +% >> fitpow.set_parameter(3, 0.02, 'background', true); % set constant background +% >> fitpow.fix_parameter(1:2, 'background', true); % fix slopes of background to 0 +% >> fitpow.cost_function = "Rsq"; % or "chisq" +% >> fitpow.optimizer = @ndbase.simplex; +% >> result = fitpow.fit('MaxIter', 1); % passes varargin to optimizer +% >> +% >> [pfit,cost_val,stat] = result{:}; +% >> fitpow.plot_result(pfit, 26, 'EdgeAlpha', 0.9, 'LineWidth', 2) +% ``` +% +% [spinw.spinwave] \| [sw_fitspec] +% properties (SetObservable) % data From c7e8704ac41256c91edd0504292cec2e062d2042 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Sat, 27 Apr 2024 16:10:14 +0100 Subject: [PATCH 18/43] Set some properties to be private --- swfiles/sw_fitpowder.m | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 841b89c86..9b4fae5a9 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -142,15 +142,17 @@ cost_function = "Rsq"; % "Rsq" or "chisq" background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) fbg - 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; % fit and parameters optimizer = @ndbase.simplex nparams_model nparams_bg params bounds - % cache calculation + end + + properties (SetAccess = private) + 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; do_cache = true ycalc_cached = [] model_params_cached = [] From b1fb8022c26fab152db92d41c4412e10ae5e366d Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 29 Apr 2024 09:28:08 +0100 Subject: [PATCH 19/43] Fix bug in getting bg param index for planar bg --- swfiles/sw_fitpowder.m | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 9b4fae5a9..9a5e63ea0 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -482,8 +482,12 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) end function iparams = get_index_of_background_parameters(obj, iparams_bg, icuts) - if any(icuts == 0) && obj.background_strategy == "independent" - icuts = 1:obj.ncuts; % apply to all cuts + if any(icuts == 0) + if obj.background_strategy == "independent" + icuts = 1:obj.ncuts; % apply to all cuts + else + icuts = 1; % will work for planar bg + end end % find index in vector of all parameters iparams = []; From 5596ea812c991456bdd4082f4d6001f5806487a3 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 29 Apr 2024 16:51:31 +0100 Subject: [PATCH 20/43] Remove Name-Val pairs from constructor and add bg estimation method --- swfiles/sw_fitpowder.m | 78 ++++++++++++++++++++---------------------- 1 file changed, 38 insertions(+), 40 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 9a5e63ea0..66917d7ff 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -39,20 +39,6 @@ % `model_params` % : Vector of initial paramers to pass to fit_func % -% ### Name-Value Pair Arguments -% -% `'backgroundStrategy'` -% : A string determining the type of background: -% * `planar` (2D planar background in |Q| and energy transfer) -% * `independent` (1D linear backgroudn as function of energy transfer) -% -% `'initialBackgroundParameters'` -% : vector containing parameters for the background (size depends on the -% dimensionality of the data and the background strategy) -% -% `'scale'` -% : Initial guess for multiplicative scale factor applied to spinwave calc -% % ### Output Arguments % % `'result'` @@ -159,46 +145,37 @@ end methods - function obj = sw_fitpowder(swobj, data, fit_func, model_params, varargin) + function obj = sw_fitpowder(swobj, data, fit_func, model_params, background_strategy) % constructor obj.swobj = swobj; obj.fit_func = fit_func; - bg_params = []; - scale = 1; - % set background startegy and func - for ikey = 1:2:numel(varargin) - switch varargin{ikey} - case 'backgroundStrategy' - obj.background_strategy = varargin{ikey+1}; - case 'initialBackgroundParameters' - bg_params = varargin{ikey+1}; - case 'scale' - scale = varargin{ikey+1}; - end - end obj.add_data(data) + if nargin < 5 + obj.background_strategy = "planar"; + else + obj.background_strategy = background_strategy; + end if obj.background_strategy == "planar" obj.fbg = obj.fbg_planar; else obj.fbg = obj.fbg_indep; end - obj.initialise_parameters_and_bounds(model_params, bg_params, scale) + obj.initialise_parameters_and_bounds(model_params) end - function initialise_parameters_and_bounds(obj, model_params, bg_params, scale) + function initialise_parameters_and_bounds(obj, model_params) obj.nparams_model = numel(model_params); obj.nparams_bg = obj.get_nparams_in_background_func(); - if isempty(bg_params) - % zero intialise - if obj.background_strategy == "independent" - nparams_bg_total = obj.ncuts * obj.nparams_bg; - else - nparams_bg_total = obj.nparams_bg; - end - bg_params = zeros(nparams_bg_total, 1); + % zero intialise background parameters + if obj.background_strategy == "independent" + nparams_bg_total = obj.ncuts * obj.nparams_bg; + else + nparams_bg_total = obj.nparams_bg; end - obj.params = [model_params(:); bg_params(:); scale]; - obj.bounds = [-inf, inf].*ones(numel(obj.params),1); + bg_params = zeros(nparams_bg_total, 1); + % last parameter is scale factor (default to 1) + obj.params = [model_params(:); bg_params(:); 1]; + obj.bounds = [-inf, inf].*ones(numel(obj.params), 1); obj.bounds(end,1) = 0; % set lower bound of scale to be 0 end @@ -423,6 +400,27 @@ function estimate_scale_factor(obj) obj.params(end) = scale; end + function estimate_constant_background(obj) + ysort = sort(obj.y(obj.y < mean(obj.y, 'omitnan')), 'descend'); + bg = ysort(int32(numel(ysort)/2)); % default to median + prev_skew = inf; + for ipt = 2:numel(ysort) + this_mean = mean(ysort(ipt:end)); + this_skew = mean((ysort(ipt:end) - this_mean).^3)/(std(ysort(ipt:end)).^3); + if this_skew < 0 + bg = this_mean; + break + else + prev_skew = this_skew; + end + end + if obj.background_strategy == "planar" && obj.ndim == 1 + bg = bg/obj.nQ; + end + % set constant background assuming last bg parameter + obj.set_bg_parameters(obj.nparams_bg, bg); % set constant background + end + function plot_result(obj, params, varargin) [ycalc, ~] = obj.calc_spinwave_spec(params); if obj.ndim == 1 From 4a3542f61a4d64e8ab36be204752338d53c0378e Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 29 Apr 2024 18:14:19 +0100 Subject: [PATCH 21/43] Support sqw, d1d and d2d objects and add validation --- swfiles/sw_fitpowder.m | 61 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 53 insertions(+), 8 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 66917d7ff..759445076 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -18,14 +18,23 @@ % % `data` % : Possible inputs depend on dimensionality of the data: -% * 2D data in either a IX_dataset_2d or struct containing fields `x`, `y` -% and `e`. The size of `y` and `e` should be (N|Q|, NEnergy) and `x` is -% a cell array of size (1,2) where x{1} contains a vector of energy -% transers and x{2} is vector of |Q| bin centers. -% * Vector of 1d datasets at constant |Q| - either IX_dataset_1d or -% structs with fields `x`, `y`, `e`, `qmin` and `qmax`, The field `x` -% is a vector of energy transfer bin centers, `qmin` and `qmax` are the -% limits of the integration in |Q|. +% * 2D data +% Either a struct containing fields `x`, `y` and `e` or a 2D HORACE +% object (sqw or d2d) +% where `y` is a matrix of intensities with shape (N|Q|, NEnergy) +% `e` is a matrix of errorsbars with shape (N|Q|, NEnergy) +% `x` is a cell array of size (1,2): +% x{1} contains a vector of energy bin centers +% x{2} contains a vector of |Q| bin centers. +% * Vector of 1d datasets at constant |Q| +% Either a struct containing fields `x`, `y` and `e` `qmin` and `qmax` +% or a 1D HORACE object (sqw or d1d) +% where `y` is a vector of intensities with shape (1, NEnergy) +% `e` is a vector of errorsbars with shape (1, NEnergy) +% `x` is a vector of energy bin centers +% x{2} contains a vector of |Q| bin centers. +% `qmin` is a scalar denoting the lower Q value of the cut +% `qmax` is a scalar denoting the upper Q value of the cut % % `fit_func` % : Function handle changing the interactions in the spinwave model. @@ -272,8 +281,14 @@ function set_scale_bounds(obj, lb, ub) end function add_data(obj, data) + if ~isa(data, "struct") + data = arrayfun(@convert_horace_to_struct, data); + end if numel(data) == 1 && numel(data.x) == 2 % 2D + assert(all(isfield(data, {'x', 'y', 'e'})), ... + 'spinw:fitpow:invalidinput', ... + 'Input cell does not have correct fields'); obj.ndim = 2; obj.y = data.y'; % nE x n|Q| obj.e = data.e'; @@ -285,6 +300,9 @@ function add_data(obj, data) obj.ebin_cens = data(1).x; obj.ncuts = numel(data); for cut = data + assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... + 'spinw:fitpow:invalidinput', ... + 'Input cell does not have correct fields'); obj.y = [obj.y cut.y(:)]; obj.e = [obj.e cut.e(:)]; obj.modQ_cens = [obj.modQ_cens, linspace(cut.qmin, cut.qmax, obj.nQ)]; % does this play nicely with sw_instrument Q convolution? @@ -494,4 +512,31 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) end end end + methods (Static=true, Hidden=true, Access = private) + function data_struct = convert_horace_to_struct(data) + if isa(data, "sqw") + ndim = data.dimensions; + assert(ndim < 3, 'spinw:fitpow:invalidinput', 'Input SQW object must be 1D or 2D'); + data_obj = data.data; + elseif isa(data, "d1d") + ndim = 1; + data_obj = data; + elseif is(data, "d2d") + ndim = 2; + data_obj = data; + else + error('spinw:fitpow:invalidinput', 'Input must be a Horace object (sqw, d1d or d2d)') + end + assert(strcmp(data_obj.ulabel{1}, '|Q|'), 'spinw:fitpow:invalidinput', 'Input Horace object is not a powder cut'); + % convert edges to bin centers + cens = cellfun(@(edges) (edges(1:end-1) + edges(2:end))/2, data_obj.p(1:ndim), 'UniformOutput', false); + data_struct = struct('x', cens, 'y', data_obj.s, 'e', data_obj.e); + if ndim == 1 + % add q integration range + data_struct.qmin = data_obj.iint(1,1); + data_struct.qmax = data_obj.iint(2,1); + end + end + + end end \ No newline at end of file From 068672cd642524f0491bd719a4cb53fe2957c1e5 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 30 Apr 2024 13:38:21 +0100 Subject: [PATCH 22/43] Add crop |Q| method for 2D data --- swfiles/sw_fitpowder.m | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 759445076..0aac56389 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -287,7 +287,7 @@ function add_data(obj, data) if numel(data) == 1 && numel(data.x) == 2 % 2D assert(all(isfield(data, {'x', 'y', 'e'})), ... - 'spinw:fitpow:invalidinput', ... + 'sw_fitpowder:invalidinput', ... 'Input cell does not have correct fields'); obj.ndim = 2; obj.y = data.y'; % nE x n|Q| @@ -301,7 +301,7 @@ function add_data(obj, data) obj.ncuts = numel(data); for cut = data assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... - 'spinw:fitpow:invalidinput', ... + 'sw_fitpowder:invalidinput', ... 'Input cell does not have correct fields'); obj.y = [obj.y cut.y(:)]; obj.e = [obj.e cut.e(:)]; @@ -320,6 +320,16 @@ function crop_energy_range(obj, emin, emax) obj.e = obj.e(ikeep, :); end + function crop_q_range(obj, qmin, qmax) + assert(obj.ndim==2, 'sw_fitpowder:invalidinput', ... + '|Q| range can only be cropped on 2D datasets'); + % crop data + ikeep = obj.modQ_cens >= qmin & obj.modQ_cens <= qmax; + obj.modQ_cens = obj.modQ_cens(ikeep); + obj.y = obj.y(:, ikeep); + obj.e = obj.e(:, ikeep); + end + function bg = calc_background(obj, params) % add background istart = obj.nparams_model + 1; From 908e7cc9c36ebffcf85e9c5f3fe4b4bb9e7dbee4 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 30 Apr 2024 15:15:26 +0100 Subject: [PATCH 23/43] Fix bug in modQ bin cens for adjacent cuts so no overlap --- swfiles/sw_fitpowder.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 0aac56389..606decd48 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -284,7 +284,7 @@ function add_data(obj, data) if ~isa(data, "struct") data = arrayfun(@convert_horace_to_struct, data); end - if numel(data) == 1 && numel(data.x) == 2 + if numel(data) == 1 && isa(data.x, "cell") % 2D assert(all(isfield(data, {'x', 'y', 'e'})), ... 'sw_fitpowder:invalidinput', ... @@ -305,7 +305,8 @@ function add_data(obj, data) 'Input cell does not have correct fields'); obj.y = [obj.y cut.y(:)]; obj.e = [obj.e cut.e(:)]; - obj.modQ_cens = [obj.modQ_cens, linspace(cut.qmin, cut.qmax, obj.nQ)]; % does this play nicely with sw_instrument Q convolution? + dQ = (cut.qmax - cut.qmin)/obj.nQ; + obj.modQ_cens = [obj.modQ_cens, (cut.qmin+dQ/2):dQ:cut.qmax]; end end obj.powspec_args.Evect = obj.ebin_cens(:)'; % row vect From 9e2b47b3119cf5c150c5896f15ed9e3124909c26 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 30 Apr 2024 15:16:11 +0100 Subject: [PATCH 24/43] Add option to specify nQ for 1D cuts in init --- swfiles/sw_fitpowder.m | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 606decd48..903b8dd21 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -154,10 +154,13 @@ end methods - function obj = sw_fitpowder(swobj, data, fit_func, model_params, background_strategy) + function obj = sw_fitpowder(swobj, data, fit_func, model_params, background_strategy, nQ) % constructor obj.swobj = swobj; obj.fit_func = fit_func; + if nargin == 6 + obj.nQ = nQ; % set this before add data + end obj.add_data(data) if nargin < 5 obj.background_strategy = "planar"; From 32c40472cc0b97307a6c7fefbffc33ee27019f25 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Tue, 30 Apr 2024 16:59:35 +0100 Subject: [PATCH 25/43] Add unit tests for sw_fitpowder --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 243 ++++++++++++++++++ 1 file changed, 243 insertions(+) create mode 100644 +sw_tests/+unit_tests/unittest_sw_fitpowder.m diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m new file mode 100644 index 000000000..d8b42f377 --- /dev/null +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -0,0 +1,243 @@ +classdef unittest_sw_fitpowder < sw_tests.unit_tests.unittest_super + % Runs through unit test for @spinw/spinwave.m + + properties + swobj = []; + data_1d_cuts = arrayfun(@(qmin) struct('x', 1:3, 'y', 1:3, ... + 'e', 1:3, 'qmin', qmin, ... + 'qmax',qmin+1), 3.5:4.5); + data_2d = struct('x', {{1:3, 4:5}}, 'y', [1:3; 1:3], ... + 'e', [1:3; 1:3]); + fit_func = @(obj, p) matparser(obj, 'param', p, 'mat', {'J_1'}, 'init', true); + j1 = 2; + default_fitpow = []; + default_fields = []; + end + + methods (TestClassSetup) + function setup_spinw_obj_and_expected_result(testCase) + % setup spinw object + testCase.swobj = sw_model('triAF', 1); + % subset of default fitpow object fields + testCase.default_fitpow = struct('y', testCase.data_2d.y', ... + 'e', testCase.data_2d.e', ... + 'ebin_cens', testCase.data_2d.x{1}, ... + 'modQ_cens', testCase.data_2d.x{2}, ... + 'params', [2;0;0;0;1], ... + 'bounds', [-Inf Inf; + -Inf Inf; + -Inf Inf; + -Inf Inf; + 0 Inf]); + testCase.default_fields = fieldnames(testCase.default_fitpow); + end + end + + methods + function verify_results(testCase, observed, expected, fieldnames) + if nargin < 4 + fieldnames = testCase.default_fields; + end + for ifld = 1:numel(fieldnames) + fld = fieldnames{ifld}; + testCase.verify_val(observed.(fld), expected.(fld)); + end + end + end + + methods (Test) + function test_init_data_2d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + testCase.verify_results(out, testCase.default_fitpow); + end + + function test_init_data_1d_planar_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.modQ_cens = 3.55:0.1:5.45; % integrtate over nQ pts + testCase.verify_results(out, expected_fitpow); + end + + function test_init_data_1d_indep_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent"); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.modQ_cens = 3.55:0.1:5.45; % integrtate over nQ pts + % add extra background param + 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_init_data_1d_specify_nQ(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "planar", 1); + testCase.verify_results(out, testCase.default_fitpow); + end + + function test_set_model_parameters(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, 5); + out.set_model_parameters(1, testCase.j1); + testCase.verify_results(out, testCase.default_fitpow); + end + + function test_set_model_parameter_bounds(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + lb = -5; + ub = 5; + out.set_model_parameter_bounds(1, -5, 5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.bounds(1, :) = [lb, ub]; + testCase.verify_results(out, expected_fitpow); + end + + function test_fix_model_parameters(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.fix_model_parameters(1); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.bounds(1, :) = [testCase.j1, testCase.j1]; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_parameters_planar_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + bg_pars = [-5, 5]; + out.set_bg_parameters(1:2, 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); + bg_pars = [-5, 5]; + out.set_bg_parameters(1:2, 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_set_bg_parameters_indep_bg_specify_icut(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(1:2, bg_pars, 2); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = [expected_fitpow.params(1:3); + bg_pars(:); 1]; + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + testCase.verify_results(out, expected_fitpow); + end + + function test_fix_bg_parameters_indep_bg_specify_icut(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + out.fix_bg_parameters(1:2, 2); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); + expected_fitpow.bounds = [expected_fitpow.bounds(1:3,:); + zeros(2,2); + expected_fitpow.bounds(end,:)]; + testCase.verify_results(out, expected_fitpow); + end + + function test_fix_bg_parameters_indep_bg_all_cuts(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + out.fix_bg_parameters(1:2); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = expected_fitpow.params([1:2,2:end],:); + expected_fitpow.bounds = [expected_fitpow.bounds(1,:); + zeros(4,2); + expected_fitpow.bounds(end,:)]; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_parameter_bounds_indep_bg_all_cuts(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + ub = 5; + out.set_bg_parameter_bounds(1:2, [], [ub, ub]); % lb unchanged + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = expected_fitpow.params([1:2,2:end]); + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + expected_fitpow.bounds(2:5, 2) = ub; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_bg_parameter_bounds_indep_bg_specify_icut(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + ub = 5; + out.set_bg_parameter_bounds(1:2, [], [ub, ub], 1); % lb unchanged + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params = expected_fitpow.params([1:2,2:end]); + expected_fitpow.bounds = expected_fitpow.bounds([1:2,2:end],:); + expected_fitpow.bounds(2:3, 2) = ub; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_scale(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + scale = 5; + out.set_scale(scale); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(end) = scale; + testCase.verify_results(out, expected_fitpow); + end + + function test_fix_scale(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.fix_scale(); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.bounds(end, :) = 1; + testCase.verify_results(out, expected_fitpow); + end + + function test_set_scale_bounds(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + lb = 0.5; + out.set_scale_bounds(0.5, []); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.bounds(end, 1) = lb; + testCase.verify_results(out, expected_fitpow); + end + + function test_crop_energy_range(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.crop_energy_range(2,5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.y = expected_fitpow.y(2:end,:); + expected_fitpow.e = expected_fitpow.e(2:end,:); + expected_fitpow.ebin_cens = expected_fitpow.ebin_cens(2:end); + testCase.verify_results(out, expected_fitpow); + testCase.verify_val(out.powspec_args.Evect, ... + expected_fitpow.ebin_cens) + end + + function test_crop_q_range(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.crop_q_range(4.5,5.5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.y = expected_fitpow.y(:, 2:end); + expected_fitpow.e = expected_fitpow.e(:, 2:end); + expected_fitpow.modQ_cens = expected_fitpow.modQ_cens(2:end); + testCase.verify_results(out, expected_fitpow); + end + end + +end From 3cf2890c93eb34ecda1e1713d4a84fffdda312f8 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 15:36:36 +0100 Subject: [PATCH 26/43] Add unit tests for estimating bg and scale --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index d8b42f377..1e9f4125d 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -34,13 +34,13 @@ function setup_spinw_obj_and_expected_result(testCase) end methods - function verify_results(testCase, observed, expected, fieldnames) + function verify_results(testCase, observed, expected, fieldnames, varargin) if nargin < 4 fieldnames = testCase.default_fields; end for ifld = 1:numel(fieldnames) fld = fieldnames{ifld}; - testCase.verify_val(observed.(fld), expected.(fld)); + testCase.verify_val(observed.(fld), expected.(fld), varargin{:}); end end end @@ -238,6 +238,28 @@ function test_crop_q_range(testCase) expected_fitpow.modQ_cens = expected_fitpow.modQ_cens(2:end); testCase.verify_results(out, expected_fitpow); end + + function test_estimate_constant_background(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.estimate_constant_background() + expected_fitpow = testCase.default_fitpow; + expected_fitpow.params(end-1) = 1; + testCase.verify_results(out, expected_fitpow); + end + + function test_estimate_scale_factor(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + 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; + testCase.verify_results(out, expected_fitpow, ... + testCase.default_fields, 'abs_tol', 1e-1); + end + end end From c3fae343b14fcd6fa39507af314caea2c91347cf Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 15:42:18 +0100 Subject: [PATCH 27/43] Add tests for calc_cost_func --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 1e9f4125d..ef19a17f4 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -260,6 +260,32 @@ function test_estimate_scale_factor(testCase) testCase.default_fields, 'abs_tol', 1e-1); end + function test_calc_cost_func_1d_planar_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "planar", 1); + out.powspec_args.dE = 0.1; % constant energy resolution + out.powspec_args.hermit = true; + cost = out.calc_cost_func(out.params); + testCase.verify_val(cost, 25.3, 'abs_tol', 0.1); + end + + function test_calc_cost_func_1d_indep_bg(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_1d_cuts, ... + testCase.fit_func, testCase.j1, "independent", 1); + out.powspec_args.dE = 0.1; % constant energy resolution + out.powspec_args.hermit = true; + cost = out.calc_cost_func(out.params); + testCase.verify_val(cost, 25.3, 'abs_tol', 0.1); + end + + function test_calc_cost_func_2d(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.powspec_args.dE = 0.1; % constant energy resolution + out.powspec_args.hermit = true; + cost = out.calc_cost_func(out.params); + testCase.verify_val(cost, 25.3, 'abs_tol', 0.1); + end end end From 83606d9d034890248e9ad8567a7258c7dff06d7a Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 16:24:22 +0100 Subject: [PATCH 28/43] Update doc strings/comments --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 1 - swfiles/sw_fitpowder.m | 22 ++++++++++++++----- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index ef19a17f4..f60c503e6 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -1,5 +1,4 @@ classdef unittest_sw_fitpowder < sw_tests.unit_tests.unittest_super - % Runs through unit test for @spinw/spinwave.m properties swobj = []; diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 903b8dd21..513ac9c33 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -48,10 +48,19 @@ % `model_params` % : Vector of initial paramers to pass to fit_func % +% `background_strategy` (optional) +% : A string determining the type of background: +% * `planar` (default) - 2D planar background in |Q| and energy transfer +% * `independent` - 1D linear background as function of energy transfer +% +% `nQ` (optional) +% : Scalar int correpsonding to number of Q-points to subdivide cuts into +% for averaging +% % ### Output Arguments % % `'result'` -% : cell array containgin output of sw_fitpowder.optimizer +% : cell array containing output of sw_fitpowder.optimizer % For ndbase optimizers in the spinw package then the reuslt can be % unpacked as follows % ``` @@ -99,15 +108,18 @@ % >> dQ = 2.35 * sqrt( (ki*a1)^2/12 + (ki*a2)^2/12 + (ki*a3)^2/12 + (ki*a4)^2/12 ); % >> % >> % fit powder -% >> fitpow = sw_fitpowder(mnf2, mnf2dat, fit_func, [J1 J2 D], 'backgroundStrategy', 'planar'); +% >> backgroundStrategy = 'planar'; % 'planar' or 'independent' +% >> nQ = 10; % Number of Q-points to subdivide cuts into for averaging +% >> +% >> fitpow = sw_fitpowder(mnf2, mnf2dat, fit_func, [J1 J2 D], backgroundStrategy, nQ); % >> fitpow.crop_energy_range(2.0, 8.0); % >> 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_parameter_bounds(3, 'lb', 0.0, 'background', true) % Set lb of constant bg = 0 -% >> fitpow.set_parameter(3, 0.02, 'background', true); % set constant background -% >> fitpow.fix_parameter(1:2, 'background', true); % fix slopes of background to 0 +% >> 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"; % or "chisq" % >> fitpow.optimizer = @ndbase.simplex; % >> result = fitpow.fit('MaxIter', 1); % passes varargin to optimizer From 2a00cc0feac5ec3574098914f09264456a98b834 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 16:24:52 +0100 Subject: [PATCH 29/43] Fix typo calling static method without obj. --- swfiles/sw_fitpowder.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 513ac9c33..116c6e27f 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -297,7 +297,7 @@ function set_scale_bounds(obj, lb, ub) function add_data(obj, data) if ~isa(data, "struct") - data = arrayfun(@convert_horace_to_struct, data); + data = arrayfun(@obj.convert_horace_to_struct, data); end if numel(data) == 1 && isa(data.x, "cell") % 2D From 0835cdbce2258a3b31aedc5b52b462608968a854 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 16:25:37 +0100 Subject: [PATCH 30/43] Set appropriate axes limits --- swfiles/sw_fitpowder.m | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 116c6e27f..21b41490e 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -482,8 +482,10 @@ function plot_1d_cuts_on_data(obj, ycalc, varargin) hold on; box on; plot(ax, obj.ebin_cens, obj.y(:,icut), 'ok'); plot(ax, obj.ebin_cens, ycalc(:,icut), '-r', varargin{:}); - xlim(ax, [0, 12]); - ylim(ax, [-2, 30]); + xlim(ax, [obj.ebin_cens(1), obj.ebin_cens(end)]); + ymin = min(min(obj.y(:,icut)), min(ycalc(:,icut))); + ymax = max(max(obj.y(:,icut)), max(ycalc(:,icut))); + ylim(ax, [ymin, ymax]); xlabel(ax, 'Energy (meV)') ylabel(ax, 'Intensity'); title(ax, num2str(modQs(icut), 2) + " $\AA^{-1}$", 'interpreter','latex') @@ -502,6 +504,8 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) cbar.Label.String = "Intensity"; xlabel(ax, "$\left|Q\right| (\AA^{-1})$", 'interpreter','latex'); ylabel(ax, "Energy (meV)"); + ylim(ax, [obj.ebin_cens(1), obj.ebin_cens(end)]); + xlim(ax, [obj.modQ_cens(1), obj.modQ_cens(end)]); legend('Calculated'); end end From 3a18555e7e674dac19f28d824713d599fe43ca01 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 16:40:36 +0100 Subject: [PATCH 31/43] Delete fitpow --- swfiles/@spinw/fitpow.m | 668 ---------------------------------------- 1 file changed, 668 deletions(-) delete mode 100644 swfiles/@spinw/fitpow.m diff --git a/swfiles/@spinw/fitpow.m b/swfiles/@spinw/fitpow.m deleted file mode 100644 index 09161dc81..000000000 --- a/swfiles/@spinw/fitpow.m +++ /dev/null @@ -1,668 +0,0 @@ -function fitsp = fitpow(obj, data, varargin) -% fits experimental powder spin wave data -% -% ### Syntax -% -% `fitsp = fitpow(obj, data, Name, Value)` -% -% ### Description -% `fitsp = fitpow(obj, data, Name, Value)` takes a set of energy cuts of powder -% spin wave spectrum and fits this. -% -% ### Input Arguments -% -% `obj` -% : [spinw] object. -% -% `data` -% : input data cuts which may be an array of any of: -% * a Horace 1D sqw or d1d object -% * a struct with fields 'x', 'y', 'e', 'qmin', 'qmax' -% * a struct with fields 'IX', 'qmin', 'qmax' -% * a cell array {x y e qmin qmax} -% * a cell array {IX_dataset_1d qmin qmax} -% where x is a vector of energy transfer, -% y is a vector of intensity, -% e is a vector of uncertainties (errorbar) -% qmin is a scalar denoting the lower Q value of the cut -% qmax is a scalar denoting the upper Q value of the cut -% IX / IX_dataset_1d is a Horace IX_dataset_1D object -% Instead of the pair of scalars qmin, qmax, you specify a -% 2-vector: [qmin qmax] under the key 'qrange', e.g. -% * a struct with fields 'x', 'y', 'e', 'qrange' -% * a cell array {x y e [qmin qmax]} -% -% ### Name-Value Pair Arguments -% -% `'func'` -% : Function to change the Hamiltonian in `obj`, it needs to have the -% following header: -% ``` -% obj = @func(obj, x) -% ``` -% -% `'xmin'` -% : Lower limit of the optimisation parameters, optional. -% -% `'xmax'` -% : Upper limit of the optimisation parameters, optional. -% -% `'x0'` -% : Starting value of the optimisation parameters. If empty -% or undefined, random values are used within the given limits. -% -% `'optimizer'` -% : String that determines the type of optimizer to use, possible values: -% * `'lm3'` (Default) Levenberg-Marquardt as implemented in Horace, see [ndbase.lm3], -% * `'pso'` Particle-swarm optimizer, see [ndbase.pso], -% * `'simplex'` Matlab built-in simplex optimizer, see [fminsearch](www.mathworks.ch/help/matlab/ref/fminsearch.html). -% -% `'nRun'` -% : Number of consecutive fitting runs, each result is saved in the -% output `fitsp.x` and `fitsp.R` arrays. If the Hamiltonian given by the -% random `x` parameters is incompatible with the ground state, -% those `x` values will be omitted and new random `x` values will be -% generated instead. Default value is 1. -% -% `'nMax'` -% : Maximum number of runs, including the ones that produce error -% (due to incompatible ground state). Default value is 1000. -% -% `'nQ'` -% : Number of |Q| points to use for each cut to simulate cut thickness -% Default is 10. -% -% `'hermit'` -% : Method for matrix diagonalization, for details see [spinw.spinwave]. -% -% `'formfact'` -% : Whether to include the form factor in the calculation (Default: true) -% -% `'formfactfun'` -% : The function used to calculate the form factor (Default: sw_mff) -% -% `'conv_func'` -% : The convolution function used to broaden the data in |Q| and E. -% Default: no convolution is applied. -% -% `'imagChk'` -% : Checks that the imaginary part of the spin wave dispersion is -% smaller than the energy bin size. -% If false, will not check -% If 'penalize' will apply a penalty to iterations that yield imaginary modes -% If true, will stop the fit if an iteration gives imaginary modes -% Default is `penalize`. -% -% `'intensity_scale'` -% : The scale factor for the intensity. Use a negative number for a fixed scale factor, -% positive for it to be fitted from the given initial value, -% and 0 for it to be fitted with an initial value guessed from the data. -% Default: -1 (fixed, no scale factor applied). -% -% `'background'` -% : The type of background to use. Either 'none', 'flat', 'linear', 'quadratic' -% -% `'background_par'` -% : The initial background parameters as a vector from highest order. -% E.g. if using ('background', 'quadratic'), then ('background_par', [a, b, c]) -% where the background used will be a*x.^2 + b*x + c. -% Default: [] (empty vector - initial parameters will be guessed from data). -% -% Parameters for monitoring the fitting: -% -% `'plot'` -% : If `true`, the data and fits are plotted as the fitting process progresses. -% Default is `true`. -% -% `'verbose'` -% : Whether to print information as the fit progresses -% Default is `true`. -% -% Optimizer options: -% -% `'TolX'` -% : Minimum change of` x` when convergence reached, default -% value is $10^{-4}$. -% -% `'TolFun'` -% : Minimum change of the R value when convergence reached, -% default value is $10^{-5}$. -% -% `'maxfunevals'` -% : Maximum number of function evaluations, default value is $10^3$. -% -% `'maxiter'` -% : Maximum number of iterations for the [ndbse.pso] optimizer. -% Default value is 20. -% -% ### Output Arguments -% -% Output `fitsp` is struct type with the following fields: -% * `obj` Copy of the input `obj`, with the best fitted -% Hamiltonian parameters. -% * `x` Final values of the fitted parameters, dimensions are -% $[n_{run}\times n_{par}]$. The rows of `x` are sorted according -% to increasing R values. -% * `redX2` Reduced $\chi^2_\eta$ value, goodness of the fit stored in a column -% vector with $n_{run}$ number of elements, sorted in increasing -% order. $\chi^2_\eta$ is defined as: -% -% $\begin{align} -% \chi^2_\eta &= \frac{\chi^2}{\eta},\\ -% \eta &= n-m+1, -% \end{align}$ -% where \\eta is the degree of freedom, $n$ number of -% observations and $m$ is the number of fitted parameters. -% -% * `exitflag` Exit flag of the `fminsearch` command. -% * `output` Output of the `fminsearch` command. -% -% {{note As a rule of thumb when the variance of the measurement error is -% known a priori, \\chi$^2_\eta$\\gg 1 indicates a poor model fit. A -% \\chi$^2_\eta$\\gg 1 indicates that the fit has not fully captured the -% data (or that the error variance has been underestimated). In principle, -% a value of \\chi$^2_\eta$= 1 indicates that the extent of the match -% between observations and estimates is in accord with the error variance. -% A \\chi$^2_\eta$ < 1 indicates that the model is 'over-fitting' the data: -% either the model is improperly fitting noise, or the error variance has -% been overestimated.}} -% -% Any other option used by [spinw.spinwave] function are also accepted. -% -% ### See Also -% -% [spinw.spinwave] \| [spinw.matparser] \| [sw_egrid] \| [sw_neutron] -% - -pref = swpref; -tid0 = pref.tid; - -if nargin < 2 - error('spinw:fitpow:MissingData','No input data cuts given') -end - -inpForm.fname = {'xmin' 'xmax' 'x0' 'func' 'fibo' 'nRand' 'verbose'}; -inpForm.defval = {[] [] [] [] true 100 true}; -inpForm.size = {[1 -3] [1 -4] [1 -5] [1 1] [1 1] [1 1] [1 1]}; -inpForm.soft = {true true true false false false false}; - -inpForm.fname = [inpForm.fname {'formfact' 'formfactfun' 'conv_func' 'nQ'}]; -inpForm.defval = [inpForm.defval {true @sw_mff @(s)s 10}]; -inpForm.size = [inpForm.size {[1 -2] [1 1] [1 1] [1 1]}]; -inpForm.soft = [inpForm.soft {false false false false}]; - -inpForm.fname = [inpForm.fname {'intensity_scale' 'background' 'background_par'}]; -inpForm.defval = [inpForm.defval {-1 'linear' []}]; -inpForm.size = [inpForm.size {[1 1] [1 -9] [1 -10]}]; -inpForm.soft = [inpForm.soft {false false true}]; - -inpForm.fname = [inpForm.fname {'tolx' 'tolfun' 'maxfunevals' 'nRun' 'plot'}]; -inpForm.defval = [inpForm.defval {1e-4 1e-5 1e3 1 true}]; -inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 1]}]; -inpForm.soft = [inpForm.soft {false false false false false}]; - -inpForm.fname = [inpForm.fname {'nMax' 'hermit' 'optimizer'}]; -inpForm.defval = [inpForm.defval {1e3 true 'lm3' }]; -inpForm.size = [inpForm.size {[1 1] [1 1] [1 -7] }]; -inpForm.soft = [inpForm.soft {false false false }]; - -inpForm.fname = [inpForm.fname {'maxiter' 'sw' 'optmem' 'tid' 'imagChk'}]; -inpForm.defval = [inpForm.defval {20 1 0 tid0 'penalize'}]; -inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1] [1 -8] }]; -inpForm.soft = [inpForm.soft {false false false false false }]; - -param = sw_readparam(inpForm, varargin{:}); - -% number of parameters (length of x) -nPar = max(max(length(param.xmin),length(param.xmax)),length(param.x0)); -nRun = param.nRun; - -% Parses the input data -[lindat, data] = parse_cuts(data); - -% setup parameters for spinw.spinwave() -param0 = param; -%param0.plot = false; -param0.tid = 0; - -% Parses the imagChk parameter - we need to set it to true/false for spinw.spinwave() -if strncmp(param.imagChk, 'penalize', 1) - param0.imagChk = true; - param0.penalize_imag = true; -else - param0.imagChk = logical(param.imagChk(1)); - param0.penalize_imag = false; -end - -optfunname = sprintf('ndbase.%s', param.optimizer); -assert(~isempty(which(optfunname)), 'spinw:fitpow:invalidoptimizer', 'Optimizer %s not recognised', param.optimizer); -optimizer = str2func(optfunname); - -%x = zeros(nRun,nPar); -redX2 = zeros(nRun,1); - -sw_timeit(0, 1, param.tid,'Fitting powder spin wave spectra'); - -idx = 1; - -while idx <= nRun - if ~isempty(param.x0) - x0 = param.x0; - elseif isempty(param.xmax) && isempty(param.xmin) - x0 = rand(1,nPar); - elseif isempty(param.xmax) - x0 = param.xmin + rand(1,nPar); - elseif isempty(param.xmin) - x0 = param.xmax - rand(1,nPar); - else - x0 = rand(1,nPar).*(param.xmax-param.xmin)+param.xmin; - end - xmin = param.xmin; - xmax = param.xmax; - - % Handle scale factor - if param.intensity_scale < 0 - param0.fit_scale = false; - else - if param.intensity_scale == 0 - sf = guess_intensity(obj, data, param0); - else - sf = param.intensity_scale; - end - param0.fit_scale = true; - x0 = [x0 sf]; - if ~isempty(xmin) - xmin = [xmin 0]; - end - if ~isempty(xmax) - xmax = [xmax 10*sf]; - end - end - - % Handle background - param0.bkgfun = []; - if ~strcmp(param.background, 'none') - xb = param.background_par; - if strcmp(param.background, 'flat') - param0.bkgfun = @(x, p) ones(numel(x),1)*p; - if isempty(xb) - xb = min(lindat.y); - end - xnx = 0; - xmx = max(lindat.y) / 2; - elseif strcmp(param.background, 'linear') - param0.bkgfun = @(x, p) x*p(1) + p(2); - if isempty(xb) - xb = guess_bkg(data, 'linear'); - end - xnx = [-abs(xb(1))*2, 0]; - xmx = [abs(xb(1))*2, max(lindat.y)/2]; - elseif strcmp(param.background, 'quadratic') - param0.bkgfun = @(x, p) x.^2*p(1) + x*p(2) + p(3); - if isempty(xb) - xb = guess_bkg(data, 'quadratic'); - end - xnx = [-abs(xb(1:2))*2, 0]; - xmx = [abs(xb(1:2))*2, max(lindat.y)/2]; - else - error('spinw:fitpow:unknownbackground', 'Background type %s not supported', param.background); - end - xb(isnan(xb)) = 0; - x0 = [x0 xb]; %#ok<*AGROW> - param0.nbkgpar = numel(xb); - if ~isempty(xmin) - xmin = [xmin xnx]; - end - if ~isempty(xmax) - xmax = [xmax xmx]; - end - end - - % Only cache powder calcs for gradient methods - if strncmp(param.optimizer, 'lm', 2) - pow_fitfun('clear_cache'); - else - pow_fitfun('no_cache'); - end - [x(idx,:),~, output(idx)] = optimizer(lindat, @(x,p)pow_fitfun(obj, lindat, data, param.func, p, param0), x0, ... - 'lb', xmin, 'ub', xmax, 'TolX', param.tolx, 'TolFun', param.tolfun, 'MaxIter', param.maxiter); - redX2(idx) = output.redX2; - - idx = idx + 1; - sw_timeit(idx/nRun*100,0,param.tid); -end - -sw_timeit(100,2,param.tid); - -% Sort results -[redX2, sortIdx] = sort(redX2); -x = x(sortIdx,:); - -% set the best fit to the spinw object -param.func(obj,x(1,1:nPar)); - -% Store all output in a struct variable. -fitsp.obj = copy(obj); -fitsp.x = x; -fitsp.redX2 = redX2; -%fitsp.exitflag = exitflag; -fitsp.output = output; - -end - -function out = horace1d_tostruct(dat) - assert(isa(dat, 'd1d') || isa(dat, 'sqw'), 'spinw:fitpow:invalidinput', 'Input must be a Horace 1D object, cell or struct'); - if isa(dat, 'sqw') - assert(dat.dimensions == 1, 'spinw:fitpow:invalidinput', 'Input SQW object must be 1D'); - dd = dat.data; - else - dd = dat; - end - assert(strcmp(dat.ulabel{1}, '|Q|'), 'spinw:fitpow:invalidinput', 'Input Horace object is not a powder cut'); - out = struct('x', (dd.p{1}(1:end-1) + dd.p{1}(2:end))/2, 'y', dd.s, 'e', dd.e, 'qmin', dd.iint(1,1), 'qmax', dd.iint(2,1)); -end - -function out = verify_cut_struct(dat) - if all(isfield(dat, {'x', 'y', 'e', 'qmin', 'qmax'})) - out = dat; - elseif all(isfield(dat, {'x', 'y', 'e', 'qrange'})) - out = struct('x', dat.x, 'y', dat.y, 'e', dat.e, 'qmin', dat.qrange(1), 'qmax', dat.qrange(2)); - elseif all(isfield(dat, {'IX', 'qmin', 'qmax'})) - out = struct('x', dat.IX.x, 'y', dat.IX.signal, 'e', dat.IX.error, 'qmin', dat.qmin, 'qmax', dat.qmax); - elseif all(isfield(dat, {'IX', 'qrange'})) - out = struct('x', dat.IX.x, 'y', dat.IX.signal, 'e', dat.IX.error, 'qmin', dat.qrange(1), 'qmax', dat.qrange(2)); - else - error('spinw:fitpow:invalidinput', 'Input structure must have fields (x,y,e) or (IX) and (qmin,qmax) or (qrange)'); - end -end - -function out = xyecell_tostruct(dat) - errmsg = 'Cell input should consist of IX_dataset objects or numeric arrays'; - if numel(dat) > 3 - assert(all(cellfun(@isnumeric, dat)), 'spinw:fitpow:invalidinput', errmsg); - xye = dat(1:3); - remd = dat(4:end); - elseif numel(dat) > 1 - assert(isa(dat{1}, 'IX_dataset_1d'), 'spinw:fitpow:invalidinput', errmsg); - xye = {dat{1}.x, dat{1}.signal, dat{1}.error}; - remd = dat(2:end); - end - if numel(remainder) == 2 - out = struct('x', xye{1}, 'y', xye{2}, 'e', xye{3}, 'qmin', remd{1}, 'qmax', remd{2}); - else - out = struct('x', xye{1}, 'y', xye{2}, 'e', xye{3}, 'qmin', remd{1}(1), 'qmax', remd{1}(2)); - end -end - -function [lindat, outstruct] = parse_cuts(input) - %outstruct = struct(); - if iscell(input) - for ii = 1:numel(input) - if isstruct(input{ii}) - outstruct(ii) = verify_cut_struct(input{ii}); - elseif iscell(input{ii}) - outstruct(ii) = xyecell_tostruct(input{ii}); - else - outstruct(ii) = horace1d_tostruct(input{ii}); - end - end - elseif isstruct(input) - for ii = 1:numel(input) - outstruct(ii) = verify_cut_struct(input(ii)); - end - else - for ii = 1:numel(input) - outstruct(ii) = horace1d_tostruct(input(ii)); - end - end - for ii = 1:numel(outstruct) - xx = outstruct(ii).x(:).'; - df = diff(xx) ./ 2; - outstruct(ii).evect = [xx(1)-df(1), xx(1:end-1) + df, xx(end) + df(end)]; - end - lindat = linearise_data(outstruct); -end - -function lindat = linearise_data(outstruct) - % Linearise the data - [x, y, e] = deal([], [], []); - is_same_e = true; - for ii = 1:numel(outstruct) - x = [x; ii+([1:numel(outstruct(ii).y)].'/1e9)]; %#ok - y = [y; outstruct(ii).y(:)]; - e = [e; outstruct(ii).e(:)]; - if ii == 1 - evec = outstruct(ii).x; - else - if ~is_same_e && (numel(evec) ~= numel(outstruct(ii).x) ... - || sum(abs(evec(:) - outstruc(ii).x(:))) > 1e-5) - is_same_e = false; - end - end - end - e(isnan(y)) = 1; - y(isnan(y)) = 0; - e(e==0) = 1; - lindat = struct('x', x, 'y', y, 'e', e, 'is_same_e', is_same_e); -end - -function out = is_same_data(d0, data) - out = numel(d0) == numel(data) && numel(d0(1).x) == numel(data(1).x) && ... - all(d0(1).x(:) == data(1).x(:)); -end - -function out = guess_bkg(data, bktype) - xx = []; yy = []; - for ii = 1:numel(data) - xx = [xx; data(ii).x(:)]; - yy = [yy; data(ii).y(:)]; - end - my = min(yy); - idx = find(yy < (max(yy)-my)/20+my); - [xx, isr] = sort(xx(idx)); - yy = yy(idx(isr)); - p2 = yy(end) / xx(end).^2; - p1 = (yy(end) - yy(1)) / (xx(end) - xx(1)); - p0 = min(yy); - if strcmp(bktype, 'linear') - out = fminsearch(@(p) sum(abs(yy - (p(1)*xx + p(2)))), [p1 p0]); - elseif strcmp(bktype, 'quadratic') - out = fminsearch(@(p) sum(abs(yy - (p(1)*xx.^2 + p(2)*xx + p(3)))), [p2 p1 p0]); - end -end - -function out = guess_intscale(swobj, data, param) - sfs = zeros(numel(data), 1); - for ii = 1:numel(data) - qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); - spec = swobj.powspec(qq, 'Evect', data(ii).evect, 'nRand', 100, 'fibo', param.fibo, ... - 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... - 'imagChk', param.imagChk); - spec = param.conv_func(spec); - yy = sum(spec.swConv, 2); - sfs(ii) = max(data(ii).y) / max(yy); - end - out = mean(sfs); -end - -function yCalc = pow_fitfun(obj, lindat, data, parfunc, x, param) -% calculates the powder spin wave spectrum to directly compare to data -% -% yCalc = SW_FITFUN(obj, data, swfunc, x, param) -% -% swobj spinw object -% data Structure contains the experimental data -% swfunc Function to change the Hamiltonian in obj, of the form: -% obj = swfunc(obj,x); -% x Actual parameter values. -% param Parameters for the spinwave function. - -persistent swp hf lines d0 useCache cache nCache iCache nEval; - -if ischar(obj) - if strcmp(obj, 'clear_cache') - cache = []; - nCache = []; - iCache = 1; - useCache = true; - elseif strcmp(obj, 'no_cache') - useCache = false; - end - nEval = 0; - return; -end - -if isempty(swp) - swp = swpref; -end -fid0 = swp.fid; -swp.fid = 0; -cleanupObj = onCleanup(@()swpref.setpref('fid', fid0)); - -% Strip out background and intensity parameters before feeding the SpinW model -x_rem = x; -if param.verbose - fprintf('Evaluation #%i with parameters:', nEval); fprintf('%f ', x_rem); - nEval = nEval + 1; -end -sf = 1.0; -if param.intensity_scale < 0 - sf = abs(param.intensity_scale); -else - if ~isempty(param.bkgfun) - bkgpar = x_rem((end-param.nbkgpar+1):end); - x_rem = x_rem(1:(end-param.nbkgpar)); - end - if param.fit_scale - sf = x_rem(end); - x_rem = x_rem(1:end-1); - end -end - -% Modify the SpinW model with new parameters -parfunc(obj, x_rem); - -% Other setup -fitstruc = data; -nD = numel(data); -nY = 0; - -if param.plot - if isempty(hf) || ~isvalid(hf) && (isempty(d0) || is_same_data(d0, data)) - hf = findobj('Tag', 'fitpow_plot'); - end - if isempty(hf) || ~isvalid(hf) - hf = figure('Tag', 'fitpow_plot'); - nR = ceil(nD / 4); - for ii = 1:nD - subplot(nR, ceil(nD / nR), ii); hold all; - errorbar(data(ii).x, data(ii).y, data(ii).e, '.k'); - title(sprintf('%4.2f < |Q| < %4.2f', data(ii).qmin, data(ii).qmax)); - box on; - end - lines = []; - d0 = data; - end - set(0, 'CurrentFigure', hf); -end - -not_cached = true; -if ~isempty(cache) && useCache - for ii = 1:numel(cache) - if x_rem == cache{ii}{1} - not_cached = false; - fitstruc = cache{ii}{2}; - break - end - end -end -if not_cached - if lindat.is_same_e - nY = nD * numel(data(1).x); - qq = zeros(nD * param.nQ); - for ii = 1:nD - i0 = (ii - 1) * param.nQ + 1; - qq(i0:(i0 + param.nQ - 1)) = linspace(data(ii).qmin, data(ii).qmax, param.nQ); - end - spec = obj.powspec(qq, 'Evect', data(1).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... - 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... - 'imagChk', param.imagChk); - spec = param.conv_func(spec); - for ii = 1:nD - i0 = (ii - 1) * param.nQ + 1; - fitstruc(ii).y = sum(spec.swConv(:, i0:(i0 + param.nQ - 1)), 2); - end - else - for ii = 1:nD - qq = linspace(data(ii).qmin, data(ii).qmax, param.nQ); - spec = obj.powspec(qq, 'Evect', data(ii).evect, 'nRand', param.nRand, 'fibo', param.fibo, ... - 'hermit', param.hermit, 'formfact', param.formfact, 'formfactfun', param.formfactfun, ... - 'imagChk', param.imagChk); - spec = param.conv_func(spec); - fitstruc(ii).y = sum(spec.swConv, 2); - nY = nY + numel(fitstruc(ii).y); - end - end - if useCache - if isempty(nCache) % Determines cache size from free memory up to max 1GB - nCache = min([numel(x_rem)*5, ceil((min([sw_freemem, 1e9]) / 4) / (numel(fitstruc) * numel(fitstruc(1).y) * 3))]); - end - cache{iCache} = {x_rem, fitstruc}; - iCache = iCache + 1; - if iCache > nCache - iCache = 1; - end - end -end - -yCalc = zeros(nY, 1); -i0 = 1; -for ii = 1:numel(data) - %if param.fit_separate - % sf = 1.0; - % ssf = fminsearch(@(x)sum(abs(data(ii).y(:) - x*fitstruc(ii).y(:)), 'omitnan'), max(data(ii).y) / max(fitstruc(ii).y)); - % if ssf < 1e5 - % fitstruc(ii).y = ssf * fitstruc(ii).y; - % end - %end - if sf ~= 1.0 - fitstruc(ii).y = fitstruc(ii).y * sf; - end - i1 = i0 + numel(fitstruc(ii).y) - 1; - if ~isempty(param.bkgfun) - yCalc(i0:i1) = fitstruc(ii).y + param.bkgfun(data(ii).x, bkgpar); - else - yCalc(i0:i1) = fitstruc(ii).y; - end - i0 = i1 + 1; -end - -%if param.fit_together -% sf = fminsearch(@(x)sum(abs(lindat.y(:) - x*yCalc(:)), 'omitnan'), max(l0.y) / max(yCalc)); -%end - -yCalc(isnan(yCalc)) = 0; -if ~param.imagChk - yCalc = abs(yCalc); -end - -if param.plot - nR = ceil(nD / 4); - i0 = 1; - for ii = 1:numel(data) - subplot(nR, ceil(nD / nR), ii); - try %#ok - delete(lines(ii)); - end - i1 = i0 + numel(fitstruc(ii).y) - 1; - lines(ii) = plot(fitstruc(ii).x, yCalc(i0:i1), '-r'); - i0 = i1 + 1; - end - drawnow; -end - -if param.verbose - fprintf('. chi2 = %d\n', sum((yCalc - lindat.y).^2 ./ lindat.e.^2) / numel(yCalc)); -end - -end From feb6cb3d37df7edbdfb4fbddafd1085cbda5c819 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 17:00:22 +0100 Subject: [PATCH 32/43] Add live plot - code form @mducle --- swfiles/sw_fitpowder.m | 33 ++++++++++++++++++++++++++------- 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 21b41490e..d836bd815 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -155,6 +155,7 @@ nparams_bg params bounds + liveplot = false end properties (SetAccess = private) @@ -416,6 +417,10 @@ function clear_cache(obj) function resid_sq_sum = calc_cost_func(obj, params) % evaluate fit function [ycalc, ~] = calc_spinwave_spec(obj, params); + if obj.liveplot + obj.plot_1d_or_2d(ycalc); + drawnow; + end resid = (obj.y - ycalc); if obj.cost_function == "chisq" resid = resid./(obj.e); @@ -426,6 +431,9 @@ function clear_cache(obj) end function result = fit(obj, varargin) + if obj.liveplot + figure("color","white"); + end % setup cell for output of ndbase optimizer/minimizer result = cell(1,nargout(obj.optimizer)); % pass params and bounds as rows for ndbase optimisers @@ -467,15 +475,15 @@ function estimate_constant_background(obj) function plot_result(obj, params, varargin) [ycalc, ~] = obj.calc_spinwave_spec(params); - if obj.ndim == 1 - obj.plot_1d_cuts_on_data(ycalc, varargin{:}) - else - obj.plot_2d_contour_on_data(ycalc, varargin{:}) - end + obj.plot_1d_or_2d(ycalc, varargin{:}); end function plot_1d_cuts_on_data(obj, ycalc, varargin) - figure("color","white"); + if ~obj.liveplot + figure("color","white"); + else + clf; + end modQs = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); for icut = 1:obj.ncuts ax = subplot(1, obj.ncuts, icut); @@ -494,7 +502,11 @@ function plot_1d_cuts_on_data(obj, ycalc, varargin) function plot_2d_contour_on_data(obj, ycalc, varargin) % varargin passed to contour - figure("color","white"); + if ~obj.liveplot + figure("color","white"); + else + clf; + end ax = subplot(1,1,1); box on; hold on; h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, obj.y); @@ -541,6 +553,13 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) iparams = [iparams, iparams_bg + obj.nparams_model + (icut-1)*obj.nparams_bg]; end end + function plot_1d_or_2d(obj, ycalc, varargin) + if obj.ndim == 1 + obj.plot_1d_cuts_on_data(ycalc, varargin{:}) + else + obj.plot_2d_contour_on_data(ycalc, varargin{:}) + end + end end methods (Static=true, Hidden=true, Access = private) function data_struct = convert_horace_to_struct(data) From ed634a9509d00a6b7a44c9c8b36cc3cc65a4f3ba Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 2 May 2024 09:33:45 +0100 Subject: [PATCH 33/43] Add liveplot_interval --- swfiles/sw_fitpowder.m | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index d836bd815..61a3f0d44 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -155,7 +155,7 @@ nparams_bg params bounds - liveplot = false + liveplot_interval = 0 end properties (SetAccess = private) @@ -164,6 +164,7 @@ do_cache = true ycalc_cached = [] model_params_cached = [] + liveplot_counter = 0 end methods @@ -417,9 +418,13 @@ function clear_cache(obj) function resid_sq_sum = calc_cost_func(obj, params) % evaluate fit function [ycalc, ~] = calc_spinwave_spec(obj, params); - if obj.liveplot - obj.plot_1d_or_2d(ycalc); - drawnow; + if obj.liveplot_interval > 0 + obj.liveplot_counter = obj.liveplot_counter + 1; + if obj.liveplot_counter == obj.liveplot_interval + obj.plot_1d_or_2d(ycalc); + drawnow; + obj.liveplot_counter = 0; + end end resid = (obj.y - ycalc); if obj.cost_function == "chisq" @@ -431,8 +436,11 @@ function clear_cache(obj) end function result = fit(obj, varargin) - if obj.liveplot + if obj.liveplot_interval > 0 figure("color","white"); + obj.liveplot_counter = 0; + else + clf; end % setup cell for output of ndbase optimizer/minimizer result = cell(1,nargout(obj.optimizer)); @@ -479,7 +487,7 @@ function plot_result(obj, params, varargin) end function plot_1d_cuts_on_data(obj, ycalc, varargin) - if ~obj.liveplot + if obj.liveplot_interval == 0 figure("color","white"); else clf; @@ -502,7 +510,7 @@ function plot_1d_cuts_on_data(obj, ycalc, varargin) function plot_2d_contour_on_data(obj, ycalc, varargin) % varargin passed to contour - if ~obj.liveplot + if obj.liveplot_interval == 0 figure("color","white"); else clf; From f5088ec997b1f4ce0552bfed8e25004118eb613f Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 2 May 2024 09:56:08 +0100 Subject: [PATCH 34/43] Make dE soft parameter in powspec --- swfiles/@spinw/powspec.m | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index fb00771c9..1521bb392 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -161,6 +161,10 @@ % * `1` Display the timing in the Command Window. % * `2` Show the timing in a separat pup-up window. % +% `dE` +% : Energy resolution (FWHM) can be function, or a numeric matrix that +% has length 1 or the number of energy bin centers. +% % The function accepts some parameters of [spinw.scga] with the most important % parameters are: % @@ -204,18 +208,22 @@ inpForm.fname = {'nRand' 'Evect' 'T' 'formfact' 'formfactfun' 'tid' 'nInt'}; inpForm.defval = {100 zeros(1,0) T0 false @sw_mff tid0 1e3 }; inpForm.size = {[1 1] [1 -1] [1 1] [1 -2] [1 1] [1 1] [1 1] }; +inpForm.soft = {false false false false false false false }; inpForm.fname = [inpForm.fname {'hermit' 'gtensor' 'title' 'specfun' 'imagChk'}]; inpForm.defval = [inpForm.defval {true false title0 @spinwave true }]; inpForm.size = [inpForm.size {[1 1] [1 1] [1 -3] [1 1] [1 1] }]; +inpForm.soft = [inpForm.soft {false false false false false}]; inpForm.fname = [inpForm.fname {'extrap' 'fibo' 'optmem' 'binType' 'component'}]; inpForm.defval = [inpForm.defval {false false 0 'ebin' 'Sperp' }]; inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 -4] [1 -5] }]; +inpForm.soft = [inpForm.soft {false false false false false}]; inpForm.fname = [inpForm.fname {'fid' , 'dE'}]; inpForm.defval = [inpForm.defval {-1, []}]; inpForm.size = [inpForm.size {[1 1], [-6, -7]}]; +inpForm.soft = [inpForm.soft {false true}]; param = sw_readparam(inpForm, varargin{:}); From 135ca599126ce6d138eb16e5be430bb12e6dbbc3 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 2 May 2024 10:11:16 +0100 Subject: [PATCH 35/43] Add fastmode and neutorn_output to powspec Set to default true in sw_fitpowder --- swfiles/@spinw/powspec.m | 25 ++++++++++++++++++++++--- swfiles/sw_fitpowder.m | 6 +++++- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index 1521bb392..93bead7d3 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -165,6 +165,17 @@ % : Energy resolution (FWHM) can be function, or a numeric matrix that % has length 1 or the number of energy bin centers. % +% `'neutron_output'` +% : If `true`, the spinwave will output only `Sperp`, the S(q,w) component +% perpendicular to Q that is measured by neutron scattering, and will +% *not* output the full Sab tensor. (Usually sw_neutron is used to +% calculate `Sperp`.) Default value is `false`. +% +% `'fastmode'` +% : If `true`, will set `'neutron_output', true`, `'fitmode', true`, +% `'sortMode', false`, and will only output intensity for positive energy +% (neutron energy loss) modes. Default value is `false`. +% % The function accepts some parameters of [spinw.scga] with the most important % parameters are: % @@ -225,6 +236,12 @@ inpForm.size = [inpForm.size {[1 1], [-6, -7]}]; inpForm.soft = [inpForm.soft {false true}]; +inpForm.fname = [inpForm.fname {'neutron_output' 'fastmode'}]; +inpForm.defval = [inpForm.defval {false false }]; +inpForm.size = [inpForm.size {[1 1] [1 1] }]; +inpForm.soft = [inpForm.soft {false false}]; + + param = sw_readparam(inpForm, varargin{:}); if param.fid == -1 @@ -312,21 +329,23 @@ warnState = warning('off','sw_readparam:UnreadInput'); specQ = param.specfun(obj,hkl,varargin{:}); warning(warnState); + specQ = sw_neutron(specQ,'pol',false); case 1 % @spinwave specQ = spinwave(obj,hkl,struct('fitmode',true,'notwin',true,... 'Hermit',param.hermit,'formfact',param.formfact,... 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... - 'optmem',param.optmem,'tid',param.tid,'fid',0),'noCheck'); + 'optmem',param.optmem,'tid',param.tid,'fid',0, ... + 'neutron_output', param.neutron_output, 'fastmode', ... + param.fastmode),'noCheck'); case 2 % @scga specQ = scga(obj,hkl,struct('fitmode',true,'formfact',param.formfact,... 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... 'fid',0,'lambda',specQ.lambda,'nInt',param.nInt,'T',param.T,... 'plot',false),'noCheck'); + specQ = sw_neutron(specQ,'pol',false); end - -specQ = sw_neutron(specQ,'pol',false); specQ.obj = obj; % use edge grid by default specQ = sw_egrid(specQ,struct('Evect',param.Evect,'T',param.T,'binType',param.binType,... diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 61a3f0d44..4148f3126 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -143,7 +143,11 @@ ncuts = 0; % spinw funtion inputs sw_instrument_args = struct() - powspec_args = struct('nRand', 1e3, 'T', 0, 'formfact', true, 'hermit', false, 'imagChk', false, 'fibo', true, 'binType', 'cbin', 'component', 'Sperp') + powspec_args = struct('nRand', 1e3, 'T', 0, 'formfact', true, ... + 'hermit', false, 'imagChk', false, ... + 'fibo', true, 'binType', 'cbin', ... + 'component', 'Sperp', 'neutron_output', true, ... + 'fastmode', true); % functions fit_func cost_function = "Rsq"; % "Rsq" or "chisq" From 589d70d69f058db5fd6fd9f6d0d5ab07e69b6e92 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 2 May 2024 13:58:32 +0100 Subject: [PATCH 36/43] Fix bugs loading horace objects --- swfiles/sw_fitpowder.m | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 4148f3126..cfa2f0bd1 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -582,15 +582,18 @@ function plot_1d_or_2d(obj, ycalc, varargin) elseif isa(data, "d1d") ndim = 1; data_obj = data; - elseif is(data, "d2d") + elseif isa(data, "d2d") ndim = 2; data_obj = data; else error('spinw:fitpow:invalidinput', 'Input must be a Horace object (sqw, d1d or d2d)') end assert(strcmp(data_obj.ulabel{1}, '|Q|'), 'spinw:fitpow:invalidinput', 'Input Horace object is not a powder cut'); - % convert edges to bin centers - cens = cellfun(@(edges) (edges(1:end-1) + edges(2:end))/2, data_obj.p(1:ndim), 'UniformOutput', false); + % convert edges to bin centers (need energy first if 2D) + cens = cellfun(@(edges) (edges(1:end-1) + edges(2:end))/2, data_obj.p(ndim:-1:1), 'UniformOutput', false); + if ndim == 2 + cens = {cens}; % otherwise MATLAB makes multiple structs + end data_struct = struct('x', cens, 'y', data_obj.s, 'e', data_obj.e); if ndim == 1 % add q integration range From bd31eb1f0f1371cc2a6a76067551743a572b6f5a Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 2 May 2024 14:49:13 +0100 Subject: [PATCH 37/43] Fix bug where not running sw_neutron if neutron_output=false --- swfiles/@spinw/powspec.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index 93bead7d3..79f644dbd 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -329,7 +329,6 @@ warnState = warning('off','sw_readparam:UnreadInput'); specQ = param.specfun(obj,hkl,varargin{:}); warning(warnState); - specQ = sw_neutron(specQ,'pol',false); case 1 % @spinwave specQ = spinwave(obj,hkl,struct('fitmode',true,'notwin',true,... @@ -344,7 +343,9 @@ 'formfactfun',param.formfactfun,'gtensor',param.gtensor,... 'fid',0,'lambda',specQ.lambda,'nInt',param.nInt,'T',param.T,... 'plot',false),'noCheck'); - specQ = sw_neutron(specQ,'pol',false); +end +if funIdx ~= 1 || (~param.neutron_output && ~param.fastmode) + specQ = sw_neutron(specQ,'pol',false); end specQ.obj = obj; % use edge grid by default From eefb777883216b7f40d84f27465f30dab8660ed9 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 11:30:45 +0100 Subject: [PATCH 38/43] Add function to exclude energy range ref #185 --- +sw_tests/+unit_tests/unittest_sw_fitpowder.m | 13 +++++++++++++ swfiles/sw_fitpowder.m | 16 ++++++++++++---- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index f60c503e6..ef752aaf0 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -227,6 +227,19 @@ function test_crop_energy_range(testCase) expected_fitpow.ebin_cens) end + function test_exclude_energy_range(testCase) + out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... + testCase.fit_func, testCase.j1); + out.exclude_energy_range(1.5,2.5); + expected_fitpow = testCase.default_fitpow; + expected_fitpow.y = expected_fitpow.y([1,end],:); + expected_fitpow.e = expected_fitpow.e([1,end],:); + expected_fitpow.ebin_cens = expected_fitpow.ebin_cens([1,end]); + testCase.verify_results(out, expected_fitpow); + testCase.verify_val(out.powspec_args.Evect, ... + expected_fitpow.ebin_cens) + end + function test_crop_q_range(testCase) out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... testCase.fit_func, testCase.j1); diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index cfa2f0bd1..f158f31e7 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -336,10 +336,12 @@ function add_data(obj, data) function crop_energy_range(obj, emin, emax) % crop data ikeep = obj.ebin_cens >= emin & obj.ebin_cens <= emax; - obj.ebin_cens = obj.ebin_cens(ikeep); - obj.powspec_args.Evect = obj.ebin_cens(:)'; - obj.y = obj.y(ikeep, :); - obj.e = obj.e(ikeep, :); + obj.apply_energy_mask(ikeep); + end + + function exclude_energy_range(obj, elo, ehi) + ikeep = obj.ebin_cens < elo | obj.ebin_cens > ehi; + obj.apply_energy_mask(ikeep); end function crop_q_range(obj, qmin, qmax) @@ -572,6 +574,12 @@ function plot_1d_or_2d(obj, ycalc, varargin) obj.plot_2d_contour_on_data(ycalc, varargin{:}) end end + function apply_energy_mask(obj, ikeep) + obj.ebin_cens = obj.ebin_cens(ikeep); + obj.powspec_args.Evect = obj.ebin_cens(:)'; + obj.y = obj.y(ikeep, :); + obj.e = obj.e(ikeep, :); + end end methods (Static=true, Hidden=true, Access = private) function data_struct = convert_horace_to_struct(data) From 1fadf9c3f79d24e68e5d3c342eb1dd4063caff3d Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 12:00:14 +0100 Subject: [PATCH 39/43] Fix bug with empty plot created when liveplot_interval=0 --- swfiles/sw_fitpowder.m | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index f158f31e7..401dc1082 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -445,8 +445,6 @@ function clear_cache(obj) if obj.liveplot_interval > 0 figure("color","white"); obj.liveplot_counter = 0; - else - clf; end % setup cell for output of ndbase optimizer/minimizer result = cell(1,nargout(obj.optimizer)); @@ -493,11 +491,6 @@ function plot_result(obj, params, varargin) end function plot_1d_cuts_on_data(obj, ycalc, varargin) - if obj.liveplot_interval == 0 - figure("color","white"); - else - clf; - end modQs = mean(reshape(obj.modQ_cens, [], obj.ncuts), 1); for icut = 1:obj.ncuts ax = subplot(1, obj.ncuts, icut); @@ -515,12 +508,6 @@ function plot_1d_cuts_on_data(obj, ycalc, varargin) end function plot_2d_contour_on_data(obj, ycalc, varargin) - % varargin passed to contour - if obj.liveplot_interval == 0 - figure("color","white"); - else - clf; - end ax = subplot(1,1,1); box on; hold on; h = imagesc(ax, obj.modQ_cens, obj.ebin_cens, obj.y); @@ -568,6 +555,11 @@ function plot_2d_contour_on_data(obj, ycalc, varargin) end end function plot_1d_or_2d(obj, ycalc, varargin) + if obj.liveplot_interval == 0 + figure("color","white"); + else + clf; + end if obj.ndim == 1 obj.plot_1d_cuts_on_data(ycalc, varargin{:}) else From 4c9a49e79a98d3ccee799e74c787f1c4eb53aa0b Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 14:06:45 +0100 Subject: [PATCH 40/43] Remove duplicate code setting parameter bounds --- swfiles/sw_fitpowder.m | 50 +++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 401dc1082..0ae8e50a5 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -213,9 +213,7 @@ function fix_model_parameters(obj, iparams) if any(iparams > obj.nparams_model) error('sw_fitpowder', 'Parameter indices supplied must be within number of model parameters'); end - for iparam = iparams - obj.bounds(iparam, :) = obj.params(iparam); - end + obj.fix_parameters( iparams) end function fix_bg_parameters(obj, iparams_bg, icuts) @@ -225,13 +223,11 @@ function fix_bg_parameters(obj, iparams_bg, icuts) if nargin < 3 icuts = 0; end - for iparam = obj.get_index_of_background_parameters(iparams_bg, icuts) - obj.bounds(iparam, :) = obj.params(iparam); - end + obj.fix_parameters(obj.get_index_of_background_parameters(iparams_bg, icuts)); end function fix_scale(obj) - obj.bounds(end, :) = obj.params(end); + obj.fix_parameters(numel(obj.params)) end function set_model_parameters(obj, iparams, values) @@ -265,12 +261,7 @@ function set_model_parameter_bounds(obj, iparams, lb, ub) error('sw_fitpowder', 'Parameter indices supplied must be within number of model parameters'); end for ibnd = 1:numel(iparams) - if ~isempty(lb) - obj.bounds(iparams(ibnd), 1) = lb(ibnd); - end - if ~isempty(ub) - obj.bounds(iparams(ibnd), 2) = ub(ibnd); - end + obj.set_bounds(iparams(ibnd), lb, ub, ibnd); end end @@ -283,22 +274,12 @@ function set_bg_parameter_bounds(obj, iparams_bg, lb, ub, icuts) end for ibnd = 1:numel(iparams_bg) iparams = obj.get_index_of_background_parameters(iparams_bg(ibnd), icuts); - if ~isempty(lb) - obj.bounds(iparams, 1) = lb(ibnd); - end - if ~isempty(ub) - obj.bounds(iparams, 2) = ub(ibnd); - end + obj.set_bounds(iparams, lb, ub, ibnd); end end function set_scale_bounds(obj, lb, ub) - if ~isempty(lb) - obj.bounds(end, 1) = lb; - end - if ~isempty(ub) - obj.bounds(end, 2) = ub; - end + obj.set_bounds(size(obj.bounds,1), lb, ub); end function add_data(obj, data) @@ -572,6 +553,25 @@ function apply_energy_mask(obj, ikeep) obj.y = obj.y(ikeep, :); obj.e = obj.e(ikeep, :); end + function set_bounds(obj, iparams, lb, ub, ibnd) + if ~isempty(lb) + if nargin == 5 + lb = lb(ibnd); + end + obj.bounds(iparams, 1) = lb; + end + if ~isempty(ub) + if nargin == 5 + ub = ub(ibnd); + end + obj.bounds(iparams, 2) = ub; + end + end + function fix_parameters(obj, iparams) + for iparam = iparams + obj.bounds(iparam, :) = obj.params(iparam); + end + end end methods (Static=true, Hidden=true, Access = private) function data_struct = convert_horace_to_struct(data) From 599d53416e0c7ed19a9818b7c98ffca70103d958 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 14:07:45 +0100 Subject: [PATCH 41/43] Fix typo - stop printing variables in lm3 if not converged --- swfiles/+ndbase/lm3.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/swfiles/+ndbase/lm3.m b/swfiles/+ndbase/lm3.m index 43f682acb..4c8ab8895 100644 --- a/swfiles/+ndbase/lm3.m +++ b/swfiles/+ndbase/lm3.m @@ -357,8 +357,8 @@ tmp = repmat(1./sqrt(diag(cov)),[1,NpFree]); cor = tmp.*cov.*tmp'; else - sigP = [] - iter=0 + sigP = []; + iter=0; chisqr = chi2Best/nnorm; ok = true; warning('WARNING: Convergence not achieved') From 6136a142c6c99a3db7cec88519102f3e097467af Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 17:23:05 +0100 Subject: [PATCH 42/43] Clear cache on changing data --- swfiles/sw_fitpowder.m | 3 +++ 1 file changed, 3 insertions(+) diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 0ae8e50a5..6a71810f5 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -286,6 +286,7 @@ function add_data(obj, data) if ~isa(data, "struct") data = arrayfun(@obj.convert_horace_to_struct, data); end + obj.clear_cache(); if numel(data) == 1 && isa(data.x, "cell") % 2D assert(all(isfield(data, {'x', 'y', 'e'})), ... @@ -333,6 +334,7 @@ function crop_q_range(obj, qmin, qmax) obj.modQ_cens = obj.modQ_cens(ikeep); obj.y = obj.y(:, ikeep); obj.e = obj.e(:, ikeep); + obj.clear_cache(); end function bg = calc_background(obj, params) @@ -552,6 +554,7 @@ function apply_energy_mask(obj, ikeep) obj.powspec_args.Evect = obj.ebin_cens(:)'; obj.y = obj.y(ikeep, :); obj.e = obj.e(ikeep, :); + obj.clear_cache(); end function set_bounds(obj, iparams, lb, ub, ibnd) if ~isempty(lb) From d72280cd6e9fdbf9d02d75d764a7ab00217744de Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 20 May 2024 17:17:42 +0100 Subject: [PATCH 43/43] Fix bug with lm3 not fixing parameters - thanks @mducle for code --- swfiles/+ndbase/lm3.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/swfiles/+ndbase/lm3.m b/swfiles/+ndbase/lm3.m index 4c8ab8895..80be94852 100644 --- a/swfiles/+ndbase/lm3.m +++ b/swfiles/+ndbase/lm3.m @@ -426,11 +426,11 @@ del=-min_abs_del; end end - if dp(j)>=0 + if dp(j) > 0 ppos=p; ppos(j)=p(j)+del; jac(:,j)=(func(dat.x,ppos)-f)/del; - else + elseif dp(j) < 0 ppos=p; ppos(j)=p(j)+del; pneg=p; pneg(j)=p(j)-del; jac(:,j)=(func(dat.x,ppos)-func(dat.x,pneg))/(2*del);