diff --git a/+sw_tests/+unit_tests/unittest_ndbase_cost_function_wrapper.m b/+sw_tests/+unit_tests/unittest_ndbase_cost_function_wrapper.m index bca8e859..56322cdb 100644 --- a/+sw_tests/+unit_tests/unittest_ndbase_cost_function_wrapper.m +++ b/+sw_tests/+unit_tests/unittest_ndbase_cost_function_wrapper.m @@ -44,8 +44,8 @@ function test_init_with_fcost_lower_bound_only(testCase) % note first param outside bounds cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'lb', [3, 1]); [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); - testCase.verify_val(pfree, [0, 3.8730], 'abs_tol', 1e-4); - testCase.verify_val(pbound, [3, 4], 'abs_tol', 1e-4); + testCase.verify_val(pfree, [1.1180, 3.8730], 'abs_tol', 1e-4); + testCase.verify_val(pbound, [3.5, 4], 'abs_tol', 1e-4); testCase.verify_val(cost_val, testCase.fcost(pbound), 'abs_tol', 1e-4); end @@ -53,8 +53,8 @@ function test_init_with_fcost_upper_bound_only(testCase) % note second param outside bounds cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'ub', [3, 1]); [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); - testCase.verify_val(pfree, [1.7320, 0], 'abs_tol', 1e-4); - testCase.verify_val(pbound, [2, 1], 'abs_tol', 1e-4); + testCase.verify_val(pfree, [1.7320, 1.1180], 'abs_tol', 1e-4); + testCase.verify_val(pbound, [2, 0.5], 'abs_tol', 1e-4); testCase.verify_val(cost_val, testCase.fcost(pbound), 'abs_tol', 1e-4); end @@ -62,8 +62,8 @@ function test_init_with_fcost_both_bounds(testCase) % note second param outside bounds cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'lb', [1, 2], 'ub', [3, 2.5]); [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); - testCase.verify_val(pfree, [0, 1.5708], 'abs_tol', 1e-4); - testCase.verify_val(pbound, [2, 2.5], 'abs_tol', 1e-4); + testCase.verify_val(pfree, [0, 0], 'abs_tol', 1e-4); + testCase.verify_val(pbound, [2, 2.25], 'abs_tol', 1e-4); testCase.verify_val(cost_val, testCase.fcost(pbound), 'abs_tol', 1e-4); end @@ -80,22 +80,31 @@ function test_init_with_fcost_both_bounds_with_fixed_param(testCase) end - function test_init_with_fcost_both_bounds_with_fixed_param_using_ifix(testCase) + function test_init_with_fcost_both_bounds_fixed_invalid_param_using_ifix(testCase) % note second param outside bounds cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'lb', [1, 2], 'ub', [3, 2.5], 'ifix', [2]); [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); testCase.verify_val(pfree, 0, 'abs_tol', 1e-4); % only first param free - testCase.verify_val(pbound, [2, 2.5], 'abs_tol', 1e-4); + testCase.verify_val(pbound, [2, 2.25], 'abs_tol', 1e-4); testCase.verify_val(cost_val, testCase.fcost(pbound), 'abs_tol', 1e-4); testCase.verify_val(cost_func_wrap.ifixed, 2); testCase.verify_val(cost_func_wrap.ifree, 1); - testCase.verify_val(cost_func_wrap.pars_fixed, 2.5); + testCase.verify_val(cost_func_wrap.pars_fixed, 2.25); + end + + function test_init_with_fcost_both_bounds_fixed_param_using_ifix(testCase) + % note second param outside bounds + cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'lb', [1, 2], 'ub', [3, 6], 'ifix', [2]); + [pfree, pbound, ~] = testCase.get_pars_and_cost_val(cost_func_wrap); + testCase.verify_val(pfree, 0, 'abs_tol', 1e-4); % only first param free + testCase.verify_val(pbound, testCase.params, 'abs_tol', 1e-4); + testCase.verify_val(cost_func_wrap.pars_fixed, testCase.params(2)); end function test_init_with_fcost_no_bounds_with_fixed_param_using_ifix(testCase) % note second param outside bounds cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'ifix', [2]); - [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); + [pfree, pbound, ~] = testCase.get_pars_and_cost_val(cost_func_wrap); testCase.verify_val(pfree, testCase.params(1), 'abs_tol', 1e-4); % only first param free testCase.verify_val(pbound, testCase.params, 'abs_tol', 1e-4); testCase.verify_val(cost_func_wrap.ifixed, 2); @@ -130,6 +139,28 @@ function test_incompatible_bounds(testCase) @() ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'lb', [1,1,], 'ub', [0,0]), ... 'ndbase:cost_function_wrapper:WrongInput'); end + + function test_init_with_resid_handle(testCase) + x = 1:3; + y = polyval(testCase.params, x); + cost_func_wrap = ndbase.cost_function_wrapper(@(p) y - polyval(p, x), testCase.params, 'resid_handle', true); + [~, ~, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); + testCase.verify_val(cost_val, 0, 'abs_tol', 1e-4); + end + + function test_init_with_fcost_all_params_fixed(testCase) + % note second param outside bounds + ifixed = 1:2; + cost_func_wrap = ndbase.cost_function_wrapper(testCase.fcost, testCase.params, 'ifix', ifixed); + [pfree, pbound, cost_val] = testCase.get_pars_and_cost_val(cost_func_wrap); + testCase.verifyEmpty(pfree); + testCase.verify_val(pbound, testCase.params); + testCase.verify_val(cost_val, testCase.fcost(pbound), 'abs_tol', 1e-4); + testCase.verify_val(cost_func_wrap.ifixed, ifixed); + testCase.verifyEmpty(cost_func_wrap.ifree); + testCase.verify_val(cost_func_wrap.pars_fixed, testCase.params); + end + end end \ No newline at end of file diff --git a/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m b/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m index f50f0a05..89974dfd 100644 --- a/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m +++ b/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m @@ -8,7 +8,7 @@ end properties (TestParameter) - optimiser = {@ndbase.simplex}; + optimiser = {@ndbase.simplex, @ndbase.lm4}; poly_func = {@(x, p) polyval(p, x), '@(x, p) polyval(p, x)'} end @@ -18,14 +18,23 @@ function test_optimise_data_struct(testCase, optimiser, poly_func) dat = struct('x', 1:3, 'e', ones(1,3)); dat.y = polyval(linear_pars, dat.x); [pars_fit, cost_val, ~] = optimiser(dat, poly_func, [-1,-1]); - testCase.verify_val(pars_fit, linear_pars, 'abs_tol', 1e-3); - testCase.verify_val(cost_val, 2.5e-7, 'abs_tol', 1e-8); + testCase.verify_val(pars_fit, linear_pars, 'abs_tol', 2e-4); + testCase.verify_val(cost_val, 0, 'abs_tol', 5e-7); + end + + function test_optimise_residual_array_lm(testCase, optimiser) + linear_pars = [2, 1]; + x = 1:3; + y = polyval(linear_pars, x); + [pars_fit, cost_val, ~] = optimiser([], @(p) y - polyval(p, x), [-1,-1], 'resid_handle', true); + testCase.verify_val(pars_fit, linear_pars, 'abs_tol', 2e-4); + testCase.verify_val(cost_val, 0, 'abs_tol', 5e-7); end function test_optimise_rosen_free(testCase, optimiser) [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [-1,-1]); testCase.verify_val(pars_fit, testCase.rosenbrock_minimum, 'abs_tol', 1e-3); - testCase.verify_val(cost_val, 0, 'abs_tol', 1e-6); + testCase.verify_val(cost_val, 0, 'abs_tol', 2e-7); end function test_optimise_rosen_lower_bound_minimum_accessible(testCase, optimiser) @@ -48,10 +57,9 @@ function test_optimise_rosen_upper_bound_minimum_accessible(testCase, optimiser) end function test_optimise_rosen_upper_bound_minimum_not_accessible(testCase, optimiser) - % note intital guess is outside bounds [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [-1,-1], 'ub', [0, inf]); testCase.verify_val(pars_fit, [0, 0], 'abs_tol', 1e-3); - testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); + testCase.verify_val(cost_val, 1, 'abs_tol', 1e-4); end function test_optimise_rosen_both_bounds_minimum_accessible(testCase, optimiser) @@ -62,7 +70,7 @@ function test_optimise_rosen_both_bounds_minimum_accessible(testCase, optimiser) function test_optimise_rosen_both_bounds_minimum_not_accessible(testCase, optimiser) % note intital guess is outside bounds - [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [-1,-1], 'lb', [-0.5, -0.5], 'ub', [0, 0]); + [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [-1,-1], 'lb', [-2, -2], 'ub', [0, 0]); testCase.verify_val(pars_fit, [0, 0], 'abs_tol', 1e-3); testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); end @@ -74,6 +82,13 @@ function test_optimise_rosen_parameter_fixed_minimum_not_accessible(testCase, op testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); end + function test_optimise_rosen_parameter_fixed_minimum_not_accessible_with_vary_arg(testCase, optimiser) + % note intital guess is outside bounds + [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [0,-1], 'lb', [nan, -0.5], 'ub', [nan, 0], 'vary', [false, true]); + testCase.verify_val(pars_fit, [0, 0], 'abs_tol', 1e-3); + testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); + end + function test_optimise_rosen_parameter_all_fixed(testCase, optimiser) % note intital guess is outside bounds [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [-1,-1], 'lb', [0, 0], 'ub', [0, 0]); @@ -81,6 +96,12 @@ function test_optimise_rosen_parameter_all_fixed(testCase, optimiser) testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); end + function test_optimise_rosen_parameter_all_fixed_with_vary_arg(testCase, optimiser) + % note intital guess is outside bounds + [pars_fit, cost_val, ~] = optimiser([], testCase.rosenbrock, [0, 0], 'vary', [false, false]); + testCase.verify_val(pars_fit, [0, 0], 'abs_tol', 1e-3); + testCase.verify_val(cost_val, 1, 'abs_tol', 1e-6); + end + end -end -% do parameterised test with all minimisers \ No newline at end of file +end \ No newline at end of file diff --git a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m index 420d9bde..925ce96d 100644 --- a/+sw_tests/+unit_tests/unittest_sw_fitpowder.m +++ b/+sw_tests/+unit_tests/unittest_sw_fitpowder.m @@ -14,6 +14,10 @@ default_modQ_cens_1d = 3.55:0.1:5.45; % integrate over nQ pts end + properties (TestParameter) + fit_params = {{}, {'resid_handle', true}}; + end + methods (TestClassSetup) function setup_spinw_obj_and_expected_result(testCase) % setup spinw object @@ -321,13 +325,13 @@ function test_estimate_constant_background(testCase) testCase.verify_results(out, expected_fitpow); end - function test_fit_background(testCase) + function test_fit_background(testCase, fit_params) out = sw_fitpowder(testCase.swobj, testCase.data_2d, ... testCase.fit_func, testCase.j1); out.y(1) = 10; % higher so other bins are background out.fix_bg_parameters(1:2); % fix slopes of background to 0 out.set_bg_parameters(3, 1.5); % initial guess - out.fit_background() + out.fit_background(fit_params{:}) expected_fitpow = testCase.default_fitpow; expected_fitpow.y(1) = 10; expected_fitpow.ibg = [3;6;2;5;4]; diff --git a/swfiles/+ndbase/cost_function_wrapper.m b/swfiles/+ndbase/cost_function_wrapper.m index 326eb8c5..cc0bb3ec 100644 --- a/swfiles/+ndbase/cost_function_wrapper.m +++ b/swfiles/+ndbase/cost_function_wrapper.m @@ -17,12 +17,15 @@ % ### Input Arguments % % `func` -% : Function handle or char with one of the following definition: -% * `R2 = func(p)` if `dat` is empty, +% : Function handle with one of the following definition: +% * `R = func(p)` if `dat` is empty, % * `y = func(x,p)` if `dat` is a struct. % Here `x` is a vector of $N$ independent variables, `p` are the -% $M$ parameters to be optimized and `y` is the simulated model, `R2` -% is the value to minimize. +% $M$ parameters to be optimized and `y` is the simulated model. +% If `resid_handle` argument is false (default) then the function returns +% a scalar (the cost function to minimise e.g. chi-squared). If +% `resid_handle` is true then the function returns a vector of residuals +% (not the residuals squared). % % `parameters` % : Vector of doubles @@ -51,10 +54,16 @@ % `ifix` % : Optional vector of ints corresponding of indices of parameters to fix % (overides bounds if provided) +% +% `resid_handle` +% : Boolean scalar - if true and `dat` is empty then fucntion handle +% returns array of residuals, if false (default) then function handle +% returns a scalar cost function. properties (SetObservable) % data cost_func + calc_resid free_to_bound_funcs bound_to_free_funcs ifixed @@ -75,23 +84,30 @@ options.ub double = [] options.data struct = struct() options.ifix = [] + options.resid_handle = false end if ischar(fhandle) fhandle = str2func(fhandle); % convert to fuction handle end if isempty(fieldnames(options.data)) - % fhandle calculates cost_val - obj.cost_func = fhandle; + if options.resid_handle + obj.calc_resid = @(p) reshape(fhandle(p), [], 1); + obj.cost_func = @(p) sum(obj.calc_resid(p).^2); + else + % fhandle calculates cost_val + obj.cost_func = fhandle; + obj.calc_resid = []; + end else % fhandle calculates fit/curve function - calc_resdid = @(p) fhandle(options.data.x(:), p) - options.data.y(:); if ~isfield(options.data,'e') || isempty(options.data.e) || ~any(options.data.e(:)) warning("ndbase:cost_function_wrapper:InvalidWeights",... "Invalid weights provided - unweighted residuals will be used.") - obj.cost_func = @(p) sum(calc_resdid(p).^2); + obj.calc_resid = @(p) fhandle(options.data.x(:), p) - options.data.y(:); else - obj.cost_func = @(p) sum((calc_resdid(p)./options.data.e(:)).^2); + obj.calc_resid = @(p) (fhandle(options.data.x(:), p) - options.data.y(:))./options.data.e(:); end + obj.cost_func = @(p) sum(obj.calc_resid(p).^2); end % validate size of bounds @@ -128,27 +144,30 @@ function init_bound_parameter_transforms(obj, pars, lb, ub, ifix) has_ub = ~isempty(ub) && isfinite(ub(ipar)); is_fixed = any(uint8(ifix) == ipar); if has_lb && has_ub - % both bounds specified and parameter not fixed obj.free_to_bound_funcs{ipar} = @(p) obj.free_to_bound_has_lb_and_ub(p, lb(ipar), ub(ipar)); obj.bound_to_free_funcs{ipar} = @(p) obj.bound_to_free_has_lb_and_ub(p, lb(ipar), ub(ipar)); % check if fixed bounds_equal = abs(ub(ipar) - lb(ipar)) < max(abs(ub(ipar)), 1)*obj.fix_tol; is_fixed = is_fixed || bounds_equal; + if is_fixed + pars(ipar) = obj.reset_bound_par_if_invalid_has_lb_and_ub(pars(ipar), lb(ipar), ub(ipar)); + end elseif has_lb obj.free_to_bound_funcs{ipar} = @(p) obj.free_to_bound_has_lb(p, lb(ipar)); obj.bound_to_free_funcs{ipar} = @(p) obj.bound_to_free_has_lb(p, lb(ipar)); + if is_fixed + pars(ipar) = obj.reset_bound_par_if_invalid_has_lb(pars(ipar), lb(ipar)); + end elseif has_ub obj.free_to_bound_funcs{ipar} = @(p) obj.free_to_bound_has_ub(p, ub(ipar)); obj.bound_to_free_funcs{ipar} = @(p) obj.bound_to_free_has_ub(p, ub(ipar)); + if is_fixed + pars(ipar) = obj.reset_bound_par_if_invalid_has_ub(pars(ipar), ub(ipar)); + end end % check fixed parameters if is_fixed obj.ifixed = [obj.ifixed, ipar]; - if has_lb && pars(ipar) < lb(ipar) - pars(ipar) = lb(ipar); - elseif has_ub && pars(ipar) > ub(ipar) - pars(ipar) = ub(ipar); - end obj.pars_fixed = [obj.pars_fixed, pars(ipar)]; end @@ -189,40 +208,58 @@ function init_bound_parameter_transforms(obj, pars, lb, ub, ifix) cost_val = obj.cost_func(pars); end + function resid = eval_resid(obj, pars_excl_fixed) + pars = obj.get_bound_parameters(pars_excl_fixed); + resid = obj.calc_resid(pars); + end + function nfree = get_num_free_parameters(obj) nfree = numel(obj.ifree); end end % private methods (Static=true, Hidden=true, Access = private) - function par_bound = free_to_bound_has_lb(par, lb) - par_bound = lb - 1 + sqrt(par^2 + 1); - end function par = bound_to_free_has_lb(par_bound, lb) - if par_bound < lb - par_bound = lb; - end + par_bound = ndbase.cost_function_wrapper.reset_bound_par_if_invalid_has_lb(par_bound, lb); par = sqrt((par_bound - lb + 1)^2 - 1); end - function par_bound = free_to_bound_has_ub(par, ub) - par_bound = ub + 1 - sqrt(par^2 + 1); - end function par = bound_to_free_has_ub(par_bound, ub) - if par_bound > ub - par_bound = ub; - end + par_bound = ndbase.cost_function_wrapper.reset_bound_par_if_invalid_has_ub(par_bound, ub); par = sqrt((ub - par_bound + 1)^2 - 1); end + function par = bound_to_free_has_lb_and_ub(par_bound, lb, ub) + par_bound = ndbase.cost_function_wrapper.reset_bound_par_if_invalid_has_lb_and_ub(par_bound, lb, ub); + par = asin((2*(par_bound - lb)/(ub-lb)) - 1); + end function par_bound = free_to_bound_has_lb_and_ub(par, lb, ub) par_bound = lb + (sin(par) + 1)*(ub-lb)/2; end - function par = bound_to_free_has_lb_and_ub(par_bound, lb, ub) + function par_bound = free_to_bound_has_ub(par, ub) + par_bound = ub + 1 - sqrt(par^2 + 1); + end + function par_bound = free_to_bound_has_lb(par, lb) + par_bound = lb - 1 + sqrt(par^2 + 1); + end + function par_bound = reset_bound_par_if_invalid_has_lb_and_ub(par_bound, lb, ub) + if par_bound < lb || par_bound > ub + warning("ndbase:cost_function_wrapper:InvalidParameter",... + "A parameter is outside bounds set - parameter will be reset.") + par_bound = (lb + ub) / 2; + end + end + function par_bound = reset_bound_par_if_invalid_has_lb(par_bound, lb) if par_bound < lb - par_bound = lb; - elseif par_bound > ub - par_bound = ub; + warning("ndbase:cost_function_wrapper:InvalidParameter",... + "A parameter is outside bounds set - parameter will be reset.") + par_bound = par_bound + max(abs(lb), 1) / 2; + end + end + function par_bound = reset_bound_par_if_invalid_has_ub(par_bound, ub) + if par_bound > ub + warning("ndbase:cost_function_wrapper:InvalidParameter",... + "A parameter is outside bounds set - parameter will be reset.") + par_bound = par_bound - max(abs(ub), 1) / 2; end - par = asin((2*(par_bound - lb)/(ub-lb)) - 1); end end end @@ -230,4 +267,5 @@ function init_bound_parameter_transforms(obj, pars, lb, ub, ifix) function isFunctionHandleOrChar(func) assert(isa(func, 'function_handle') || ischar(func)); -end \ No newline at end of file +end + diff --git a/swfiles/+ndbase/estimate_hessian.m b/swfiles/+ndbase/estimate_hessian.m index bc35034f..dd0e4548 100644 --- a/swfiles/+ndbase/estimate_hessian.m +++ b/swfiles/+ndbase/estimate_hessian.m @@ -68,6 +68,10 @@ % considered appropriate. Step size optimisation will be terminated once % this condition is met. Default is sqrt(epse) ~ 1e-8. % +% `cost_val` (optional) +% : Initial cost function evaluated at input parameters to avoid +% additional function call if already available. +% % ### Output Arguments % % `hessian` @@ -79,7 +83,8 @@ % * `cost_val`: Cost function evaluted at input parameters % * `step_size`: Vector of steps in finite-difference method for each % parameter. -% +% * `jacobian`: Vector of jacobian at input parameters (dcost/dp) + % ### Examples %``` % >> % make data @@ -108,6 +113,7 @@ options.ivary double = 0 options.niter (1,1) double = 16 options.cost_tol (1,1) double = sqrt(eps) + options.cost_val (1,1) double = nan end optimise_step = true; min_step = sqrt(eps); @@ -119,7 +125,7 @@ step_size = options.step; if numel(step_size)==1 % interpret as fractional step size - step_size = params .* step_size; + step_size = abs(params) .* step_size; end end if options.ivary @@ -140,7 +146,11 @@ % calculate jacobian at param hessian = zeros(npar); - initial_cost = fcost(params); + if isfinite(options.cost_val) + initial_cost = options.cost_val; + else + initial_cost = fcost(params); + end jac_one_step = zeros(1, npar); cost_one_step = zeros(1, npar); for irow = 1:npar @@ -183,5 +193,6 @@ end out_struct.cost_val = initial_cost; out_struct.step_size = step_size; + out_struct.jacobian = jac_one_step; varargout = {out_struct}; end diff --git a/swfiles/+ndbase/lm4.m b/swfiles/+ndbase/lm4.m new file mode 100644 index 00000000..ace81059 --- /dev/null +++ b/swfiles/+ndbase/lm4.m @@ -0,0 +1,297 @@ +function [pOpt, fVal, stat] = lm4(dat, func, p0, varargin) +% optimization of parameters using the Levenberg Marquardt method +% +% ### Syntax +% +% `[pOpt,fVal,stat] = ndbase.simplex([],func,p0,Name,Value)` +% +% `[pOpt,fVal,stat] = ndbase.simplex(dat,func,p0,Name,Value)` +% +% ### Input Arguments +% +% `dat` +% : Either empty or contains data to be fitted stored in a structure with +% fields: +% * `dat.x` vector of $N$ independent variables, +% * `dat.y` vector of $N$ data values to be fitted, +% * `dat.e` vector of $N$ standard deviation (positive numbers) +% used to weight the fit. If zero or missing +% `1/dat.y^2` will be assigned to each point. +% +% `func` +% : Function handle with one of the following definition: +% * `R = func(p)` if `dat` is empty, +% * `y = func(x,p)` if `dat` is a struct. +% Here `x` is a vector of $N$ independent variables, `p` are the +% $M$ parameters to be optimized and `y` is the simulated model. +% If `resid_handle` argument is false (default) then the function returns +% a scalar (the cost funciton to minimise e.g. chi-squared). If +% `resid_handle` is true then the function returns a vector of residuals +% (not the residuals squared). +% +% `p0` +% : vector of initial parameter guesses - starting point for the +% optimisation. +% +% ### Name-Value Pair Arguments +% +% `'lb'` +% : Vector with $M$ elements, lower boundary of the parameters. Default +% value is -inf (i.e. unbounded). +% +% `'ub'` +% : Vector with $M$ elements, upper boundary of the parameters. Default +% value is inf (i.e. unbounded). +% +% `'resid_handle'` +% : Boolean scalar - if true and `dat` is empty then 'func' fucntion handle +% returns array of residuals, if false (default) then function handle +% returns a scalar cost function. +% +% `'diff_step'` +% : Vector with $M$ or 1 element, defines the fractional increment of +% a parameter when calculating the Jacobians using the forward finite +% difference (default is 1e-7). +% +% `'MaxIter'` +% : Maximum number of iterations, default value is $100M$. +% +% `'gtol'` +% : Convergence tolerance on gradient vector of cost-function wrt change +% in parameters (default 1e-8). +% +% `'ftol'` +% : Convergence tolerance on change in cost function (default 1e-8). +% +% `'ptol'` +% : Convergence tolerance on relative length of parameter step vector wrt +% parameter vector (default 1e-8). +% +% `'lambda0'` +% : Initial Levenberg–Marquardt damping parameter which determines step +% size and angle of trajectory between gradient-descent (initially) and +% Gauss-Newton (where cost-function surface becomes increasingly +% parabolic). Initially `large` lambda values correponds to smaller steps +% along gradient descent trajectory. Decreasing `'lambda0'` may lead to +% faster convergence if the solution is near the minimum but will be less +% reliable further from the minimum where the second-order expansion is +% less valid. Default value is 1e-2. +% +% `'nu_up'` +% : Factor by which to scale the damping parameter (lambda) if parameter +% step increases the cost value. Typically > 1 (default is 5) - ie. +% increasing lambda to produce smaller steps such that the +% first-order approx. is more valid. +% +% `'nu_dn'` +% : Factor by which to scale the damping parameter (lambda) if parameter +% step decreases the cost value. Typically < 1 (default is 0.3) - i.e. +% decreasing lambda to speed-up convergence by taking larger steps as +% the cost-function surface becomes increasingly parabolic. +% +% `'vary'` +% : Boolean vector with $M$ elements (one per parameter), if an element is +% false, the corresponding parameter will be fixed (not optimised). +% Default is true for all parameters. +% +% ### Output +% +% `pOpt` +% : Value of the $M$ optimal parameters. +% +% `fVal` +% : Value of the cost function calculated at the optimal parameters +% +% `stat` +% : Structure storing the detailed output of the calculation with +% the following fields: +% p Least-squares optimal estimate of the parameter values. +% redX2 Reduced Chi squared statistic, its value +% should be close to 1. If the value is larger, the +% model is not a good description of the data. If the +% value is smaller, the model is overparameterized +% and fitting the statistical error of the data. +% sigP Asymptotic standard error of the parameters. +% sigY Asymptotic standard error of the curve-fit (to implement). +% corrP Correlation matrix of the parameters (to implement). +% Rsq R-squared cofficient of multiple determination (to implement). +% cvgHst Convergence history (to implement).. +% exitFlag The reason, why the code stopped: +% 0 maximum number of iterations reached +% 1 convergence step-size (see `ptol`) +% 2 convergence in cost (see `ftol`) +% 3 convergence in gradient (see `gtol`) +% msg String, one of the above messages. +% nIter The number of iterations executed during the fit. +% param Input parameters passed to ndbase.lm4 +% algorithm name of algorithm (Levenberg Marquardt); +% func Function handle corresponding to cost-function +% +% See also ndbase.simplex. + +nparams = numel(p0); + +inpForm.fname = {'diff_step' 'lb' 'ub' 'MaxIter' }; +inpForm.defval = {1e-7 -inf(1,nparams) inf(1,nparams) 100*nparams}; +inpForm.size = {[1 -1] [1 nparams] [1 nparams] [1 1]}; + +inpForm.fname = [inpForm.fname {'gtol' 'ftol' 'ptol' 'lambda0'}]; +inpForm.defval = [inpForm.defval {1e-8 1e-8 1e-8 1e-2}]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [1 1]}]; + +inpForm.fname = [inpForm.fname {'nu_up', 'nu_dn', 'resid_handle'}]; +inpForm.defval = [inpForm.defval {5 0.3, false}]; +inpForm.size = [inpForm.size {[1 1] [1 1], [1 1]}]; + +inpForm.fname = [inpForm.fname {'vary'}]; +inpForm.defval = [inpForm.defval {true(1, nparams)}]; +inpForm.size = [inpForm.size {[1 nparams]}]; + +param = sw_readparam(inpForm, varargin{:}); + +cost_func_wrap = ndbase.cost_function_wrapper(func, p0, "data", dat, ... + 'lb', param.lb, 'ub', param.ub, ... + 'resid_handle', param.resid_handle, ... + 'ifix', find(~param.vary)); +% transform starting values into their unconstrained surrogates. +p0_free = cost_func_wrap.get_free_parameters(p0); + +if isempty(dat) && ~param.resid_handle + % minimising scalar - empty resid + eval_cost_func = @eval_cost_scalar; + calc_hessian_and_jacobian = @calc_hessian_and_jacobian_scalar; + diff_step = param.diff_step; % always interpreted as fractional step even if 1 parameter +else + % minimising sum square residuals + eval_cost_func = @eval_cost_resids; + calc_hessian_and_jacobian = @calc_hessian_and_jacobian_resid; + % get absolute diff_step for each parameter + diff_step = abs(p0_free).*param.diff_step; + min_step = sqrt(eps); + diff_step(abs(diff_step) < min_step) = min_step; +end + +% eval at starting guess +[cost_val, resids] = eval_cost_func(cost_func_wrap, p0_free(:)); +if isempty(resids) + ndof = 1; % minimising scalar function +else + ndof = numel(resids) - numel(p0_free) + 1; +end +if isempty(p0_free) + % Re-calc bound params as p0 as could have fixed params outside bounds + pOpt = cost_func_wrap.get_bound_parameters(p0_free); + fVal = cost_val/ndof; + message = 'Parameters are fixed, no optimisation'; + perr = zeros(size(pOpt)); + niter = 0; + exit_flag = 0; +else + p = p0_free(:); + lambda = param.lambda0; + [hess, jac] = calc_hessian_and_jacobian_scalar(cost_func_wrap, p, diff_step, cost_val, resids); + exit_flag = 0; + for niter = 1:param.MaxIter + dp = calc_parameter_step(hess, jac, lambda); + if norm(dp) < param.ptol*(param.ptol + norm(p)) + message = "step size below tolerance ptol"; + exit_flag = 1; + break + end + new_p = p - dp; + new_cost_val = cost_func_wrap.eval_cost_function(new_p); + dcost = new_cost_val - cost_val; + if dcost < 0 + lambda = param.nu_dn*lambda; % decrease lambda + p = new_p; + [cost_val, resids] = eval_cost_func(cost_func_wrap, p); + [hess, jac] = calc_hessian_and_jacobian(cost_func_wrap, p, diff_step, cost_val, resids); + if abs(dcost) < param.ftol*ndof + message = "change in reduced chi-sq below tolerance ftol"; + exit_flag = 2; + break + elseif norm(jac) < param.gtol + message = "change in gradient vector below tolerance gtol"; + exit_flag = 3; + break + end + else + % increase lambda so take path closer to steepest descent + % don't recalc hess or jac + lambda = param.nu_up*lambda; + end + end + % collect output + pOpt = cost_func_wrap.get_bound_parameters(p); + fVal = cost_val / ndof; + if exit_flag > 0 + % converged on solution - calculate errors + cov = pinv(hess) * 2.0 * fVal; + perr = sqrt(diag(cov)); + else + message = "Failed to converge in MaxIter"; + perr = zeros(size(pOpt)); + end +end + +stat.Rsq = []; +stat.sigY = []; +stat.corrP = []; +stat.cvgHst = []; +stat.algorithm = 'Levenberg Marquardt'; +stat.func = cost_func_wrap.cost_func; +stat.param = param; +stat.param.Np = cost_func_wrap.get_num_free_parameters(); +stat.msg = message; +stat.iterations = niter; +stat.exitFlag = exit_flag; +stat.p = pOpt; +stat.redX2 = fVal; +stat.sigP = perr; + +end + +function dp = calc_parameter_step(hess, jac, lambda) + damped_hess = hess + lambda*diag(diag(hess)); % LM + if rcond(damped_hess) < eps + dp = pinv(damped_hess)*jac; % nearly singular + else + try + dp = damped_hess\jac; + catch + dp = pinv(damped_hess)*jac; + end + end +end + +function [hess, jac] = calc_hessian_and_jacobian_scalar(cost_func_wrap, p, diff_step, cost_val, ~) + % last input ignored (empty resids vector for consistent API) + [hess, extra_info] = ndbase.estimate_hessian(@cost_func_wrap.eval_cost_function, p, 'cost_val', cost_val, 'step', diff_step); + jac = extra_info.jacobian(:); +end + +function [cost_val, resids] = eval_cost_scalar(cost_func_wrap, p) + cost_val = cost_func_wrap.eval_cost_function(p); + resids = []; +end + +function [hess, jac] = calc_hessian_and_jacobian_resid(cost_func_wrap, p, diff_step, ~, resids) + % evaluate jacobian of residuals using finite difference + jac_resids = ones([numel(resids), numel(p)]); + for ipar = 1:numel(p) + p(ipar) = p(ipar) + diff_step(ipar); + resids_one_step = cost_func_wrap.eval_resid(p); + jac_resids(:,ipar) = (resids_one_step - resids)/diff_step(ipar); + p(ipar) = p(ipar) - diff_step(ipar); + end + % eval jacobian of cost function + jac = 2*(jac_resids')*resids; + % approx. hessian of cost fucntion + hess = 2*(jac_resids')*jac_resids; +end + +function [cost_val, resids] = eval_cost_resids(cost_func_wrap, p) + resids = cost_func_wrap.eval_resid(p); + cost_val = sum(resids(:).^2); +end + diff --git a/swfiles/+ndbase/simplex.m b/swfiles/+ndbase/simplex.m index 4eed1054..c0f8059a 100644 --- a/swfiles/+ndbase/simplex.m +++ b/swfiles/+ndbase/simplex.m @@ -90,11 +90,18 @@ % % `func` % : Function handle with one of the following definition: -% * `R2 = func(p)` if `dat` is empty, +% * `R = func(p)` if `dat` is empty, % * `y = func(x,p)` if `dat` is a struct. % Here `x` is a vector of $N$ independent variables, `p` are the -% $M$ parameters to be optimized and `y` is the simulated model, `R2` -% is the value to minimize. +% $M$ parameters to be optimized and `y` is the simulated model. +% If `resid_handle` argument is false (default) then the function returns +% a scalar (the cost funciton to minimise e.g. chi-squared). If +% `resid_handle` is true then the function returns a vector of residuals +% (not the residuals squared). +% +% `p0` +% : vector of initial parameter guesses - starting point for the +% optimisation. % % ### Name-Value Pair Arguments % @@ -126,11 +133,20 @@ % the weighted least square deviation from data). Default value is % $10^{-3}$. % +% `'resid_handle'` +% : Boolean scalar - if true and `dat` is empty then 'func' fucntion handle +% returns array of residuals, if false (default) then function handle +% returns a scalar cost function. +% +% `'vary'` +% : Boolean vector with $M$ elements (one per parameter), if an element is +% false, the corresponding parameter will be fixed (not optimised). +% Default is true for all parameters. % % ### See Also % % [ndbase.lm] \| [ndbase.pso] - +% % Original author: John D'Errico % E-mail: woodchips@rochester.rr.com % Release: 4 @@ -149,10 +165,17 @@ inpForm.defval = {'off' 1e-3 1e-3 100*Np [] [] }; inpForm.size = {[1 -1] [1 1] [1 1] [1 1] [-5 -2] [-3 -4] }; inpForm.soft = {false false false false true true }; + + inpForm.fname = [inpForm.fname {'resid_handle', 'vary'}]; + inpForm.defval = [inpForm.defval {false, true(1, Np)}]; + inpForm.size = [inpForm.size {[1 1], [1, Np]}]; param = sw_readparam(inpForm, varargin{:}); - cost_func_wrap = ndbase.cost_function_wrapper(func, p0, "data", dat, 'lb', param.lb, 'ub', param.ub); + cost_func_wrap = ndbase.cost_function_wrapper(func, p0, "data", dat, ... + 'lb', param.lb, 'ub', param.ub, ... + 'resid_handle', param.resid_handle, ... + 'ifix', find(~param.vary)); % transform starting values into their unconstrained surrogates. p0_free = cost_func_wrap.get_free_parameters(p0); diff --git a/swfiles/sw_fitpowder.m b/swfiles/sw_fitpowder.m index 8a3fda59..e98efd89 100644 --- a/swfiles/sw_fitpowder.m +++ b/swfiles/sw_fitpowder.m @@ -452,7 +452,7 @@ function clear_cache(obj) end - function resid_sq_sum = calc_cost_func(obj, params) + function resid = calc_residuals(obj, params) % evaluate fit function [ycalc, ~] = obj.calc_spinwave_spec(params); if obj.liveplot_interval > 0 @@ -463,16 +463,26 @@ function clear_cache(obj) obj.liveplot_counter = 0; end end - resid_sq_sum = obj.eval_cost(obj.y, obj.e, ycalc); + resid = obj.eval_residuals(obj.y, obj.e, ycalc); end - function resid_sq_sum = calc_cost_func_of_background(obj, bg_params) + function resid_sq_sum = calc_cost_func(obj, params) + resid = obj.calc_residuals(params); + resid_sq_sum = resid'*resid; + end + + function resid = calc_residuals_of_background(obj, bg_params) bg = obj.calc_background(bg_params); if obj.ndim == 1 && obj.background_strategy=="planar" % integrate nQ |Q| points for each cut bg = obj.rebin_powspec_to_1D_cuts(bg); end - resid_sq_sum = obj.eval_cost(obj.y(obj.ibg), obj.e(obj.ibg), bg(obj.ibg)); + resid = obj.eval_residuals(obj.y(obj.ibg), obj.e(obj.ibg), bg(obj.ibg)); + end + + function resid_sq_sum = calc_cost_func_of_background(obj, bg_params) + resid = obj.calc_residuals_of_background(bg_params); + resid_sq_sum = resid'*resid; end function varargout = fit(obj, varargin) @@ -483,7 +493,12 @@ function clear_cache(obj) % setup cell for output of ndbase optimizer/minimizer varargout = cell(1,nargout(obj.optimizer)); % pass params and bounds as rows for ndbase optimisers - [varargout{:}] = obj.optimizer([], @obj.calc_cost_func, obj.params(:)', ... + if any(cellfun(@(elem) elem =="resid_handle", varargin(1:2:end))) + fobj = @obj.calc_residuals; % minimise using residual array + else + fobj = @obj.calc_cost_func; % minimise scalar + end + [varargout{:}] = obj.optimizer([], fobj, obj.params(:)', ... 'lb', obj.bounds(:,1)', ... 'ub', obj.bounds(:,2)', ... varargin{:}); @@ -534,7 +549,13 @@ function fit_background(obj, varargin) error('spinw:sw_fitpowder:fit_background', ... 'Not enough points to fit the function.'); end - [bg_params, ~, stat] = ndbase.simplex([], @obj.calc_cost_func_of_background, ... + if any(cellfun(@(elem) elem =="resid_handle", varargin(1:2:end))) + fobj = @obj.calc_residuals_of_background; % minimise using residual array + else + fobj = @obj.calc_cost_func_of_background; % minimise scalar + end + + [bg_params, ~, stat] = ndbase.simplex([], fobj, ... obj.params(ibg_par)', ... 'lb', obj.bounds(ibg_par,1)', ... 'ub', obj.bounds(ibg_par,2)', ... @@ -662,14 +683,13 @@ function plot_1d_cuts_of_2d_data(obj, qmins, qmaxs, params) % private methods (Hidden=true, Access = private) - function resid_sq_sum = eval_cost(obj, y, e, ycalc) + function resid = eval_residuals(obj, y, e, ycalc) resid = (y - ycalc); if obj.cost_function == "chisq" resid = resid./e; end % exclude nans in both ycalc and input data - ikeep = isfinite(resid) & e > obj.zero_abs_tol; - resid_sq_sum = resid(ikeep)'*resid(ikeep); + resid = resid(isfinite(resid) & e > obj.zero_abs_tol); end function cut = cut_2d_data(obj, qmin, qmax)