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..ef752aaf0 --- /dev/null +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -0,0 +1,303 @@ +classdef unittest_sw_fitpowder < sw_tests.unit_tests.unittest_super + + 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, varargin) + if nargin < 4 + fieldnames = testCase.default_fields; + end + for ifld = 1:numel(fieldnames) + fld = fieldnames{ifld}; + testCase.verify_val(observed.(fld), expected.(fld), varargin{:}); + 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_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); + 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 + + 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 + + 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 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/+ndbase/lm3.m b/swfiles/+ndbase/lm3.m index 9cd3a8abc..80be94852 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') @@ -424,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); diff --git a/swfiles/@spinw/powspec.m b/swfiles/@spinw/powspec.m index 6f6e9fdb4..79f644dbd 100644 --- a/swfiles/@spinw/powspec.m +++ b/swfiles/@spinw/powspec.m @@ -161,6 +161,21 @@ % * `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. +% +% `'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: % @@ -204,18 +219,28 @@ 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}]; + +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}]; -inpForm.fname = [inpForm.fname {'fid'}]; -inpForm.defval = [inpForm.defval {-1 }]; -inpForm.size = [inpForm.size {[1 1]}]; param = sw_readparam(inpForm, varargin{:}); @@ -246,7 +271,6 @@ end nQ = numel(hklA); -powSpec = zeros(max(1,nE),nQ); fprintf0(fid,'Calculating powder spectra...\n'); @@ -260,6 +284,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 +304,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 +318,40 @@ 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 -sw_timeit(100,2,param.tid); +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, ... + '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'); +end +if funIdx ~= 1 || (~param.neutron_output && ~param.fastmode) + specQ = sw_neutron(specQ,'pol',false); +end +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, 'dE', param.dE),'noCheck'); +powSpec = squeeze(sum(reshape(specQ.swConv, max(1, nE), nk, nQ), 2) / param.nRand); fprintf0(fid,'Calculation finished.\n'); @@ -377,4 +404,4 @@ F = num(end-1); F1 = num(end-2); end -end \ No newline at end of file +end diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m new file mode 100644 index 000000000..6a71810f5 --- /dev/null +++ b/swfiles/sw_fitpowder.m @@ -0,0 +1,609 @@ +classdef sw_fitpowder < handle & matlab.mixin.SetGet +% 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 +% 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. +% 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 +% +% `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 containing 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 +% >> 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_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 +% >> +% >> [pfit,cost_val,stat] = result{:}; +% >> fitpow.plot_result(pfit, 26, 'EdgeAlpha', 0.9, 'LineWidth', 2) +% ``` +% +% [spinw.spinwave] \| [sw_fitspec] +% + + properties (SetObservable) + % data + swobj; + y + e + ebin_cens + modQ_cens + 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', 'neutron_output', true, ... + 'fastmode', true); + % functions + fit_func + cost_function = "Rsq"; % "Rsq" or "chisq" + background_strategy = "planar" % "planar" or "independent" (1D only - fbg = @(en, p1, p2, ..., pN) + fbg + % fit and parameters + optimizer = @ndbase.simplex + nparams_model + nparams_bg + params + bounds + liveplot_interval = 0 + 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 = [] + liveplot_counter = 0 + end + + methods + 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"; + 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) + end + + function initialise_parameters_and_bounds(obj, model_params) + obj.nparams_model = numel(model_params); + obj.nparams_bg = obj.get_nparams_in_background_func(); + % 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 + 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 + + 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 + obj.fix_parameters( iparams) + end + + 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 + obj.fix_parameters(obj.get_index_of_background_parameters(iparams_bg, icuts)); + end + + function fix_scale(obj) + obj.fix_parameters(numel(obj.params)) + 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) + obj.set_bounds(iparams(ibnd), lb, ub, ibnd); + end + 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); + obj.set_bounds(iparams, lb, ub, ibnd); + end + end + + function set_scale_bounds(obj, lb, ub) + obj.set_bounds(size(obj.bounds,1), lb, ub); + end + + 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'})), ... + 'sw_fitpowder:invalidinput', ... + 'Input cell does not have correct fields'); + 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 + assert(all(isfield(cut, {'x', 'y', 'e', 'qmin','qmax'})), ... + 'sw_fitpowder:invalidinput', ... + 'Input cell does not have correct fields'); + obj.y = [obj.y cut.y(:)]; + obj.e = [obj.e cut.e(:)]; + 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 + end + + function crop_energy_range(obj, emin, emax) + % crop data + ikeep = obj.ebin_cens >= emin & obj.ebin_cens <= emax; + 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) + 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); + obj.clear_cache(); + 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 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) + 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)*ycalc; + % calc background + bg = calc_background(obj, params); + 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) + % evaluate fit function + [ycalc, ~] = calc_spinwave_spec(obj, params); + 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" + resid = resid./(obj.e); + end + % exclude nans in both ycalc and input data + ikeep = ~isnan(resid); + resid_sq_sum = resid(ikeep)'*resid(ikeep); + end + + function result = fit(obj, varargin) + if obj.liveplot_interval > 0 + figure("color","white"); + obj.liveplot_counter = 0; + end + % setup cell for output of ndbase optimizer/minimizer + result = cell(1,nargout(obj.optimizer)); + % 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 + + function estimate_scale_factor(obj) + % set scale factor to 1 + params = obj.params; + params(end) = 1; + [ycalc, bg] = calc_spinwave_spec(obj, params); + scale = max(obj.y - mean(bg(:)), [], "all")/max(ycalc, [], "all"); + 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); + obj.plot_1d_or_2d(ycalc, varargin{:}); + end + + function plot_1d_cuts_on_data(obj, ycalc, varargin) + 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, [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') + end + end + + function plot_2d_contour_on_data(obj, ycalc, varargin) + 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)"); + ylim(ax, [obj.ebin_cens(1), obj.ebin_cens(end)]); + xlim(ax, [obj.modQ_cens(1), obj.modQ_cens(end)]); + legend('Calculated'); + 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 iparams = get_index_of_background_parameters(obj, iparams_bg, icuts) + 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 = []; + for icut = icuts + iparams = [iparams, iparams_bg + obj.nparams_model + (icut-1)*obj.nparams_bg]; + 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 + 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, :); + obj.clear_cache(); + 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) + 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 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 (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 + 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 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 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 fc34ceed7..46daab4a7 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); @@ -109,6 +110,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