From 5b7258535514557fd0132aec574dff3b1d6cacc5 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Fri, 15 Nov 2024 19:27:35 +0000 Subject: [PATCH] Add vary arg to simplex and lm4 (and tests) --- .../+unit_tests/unittest_ndbase_optimisers.m | 14 +++++++++++++ swfiles/+ndbase/cost_function_wrapper.m | 11 ++++++---- swfiles/+ndbase/lm4.m | 16 +++++++++++++-- swfiles/+ndbase/simplex.m | 20 +++++++++++++------ 4 files changed, 49 insertions(+), 12 deletions(-) diff --git a/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m b/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m index 62bcaabf..89974dfd 100644 --- a/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m +++ b/+sw_tests/+unit_tests/unittest_ndbase_optimisers.m @@ -82,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]); @@ -89,5 +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 \ No newline at end of file diff --git a/swfiles/+ndbase/cost_function_wrapper.m b/swfiles/+ndbase/cost_function_wrapper.m index b0ef6218..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 diff --git a/swfiles/+ndbase/lm4.m b/swfiles/+ndbase/lm4.m index 9100e248..924f2cd0 100644 --- a/swfiles/+ndbase/lm4.m +++ b/swfiles/+ndbase/lm4.m @@ -25,7 +25,7 @@ % 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 minimsie e.g. chi-squared). If +% 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). % @@ -89,6 +89,11 @@ % 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` @@ -138,9 +143,16 @@ 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); +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/+ndbase/simplex.m b/swfiles/+ndbase/simplex.m index b99e4e91..d18f557f 100644 --- a/swfiles/+ndbase/simplex.m +++ b/swfiles/+ndbase/simplex.m @@ -95,7 +95,7 @@ % 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 minimsie e.g. chi-squared). If +% 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). % @@ -138,10 +138,15 @@ % 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 @@ -161,13 +166,16 @@ 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'}]; - inpForm.defval = [inpForm.defval {false}]; - inpForm.size = [inpForm.size {[1 1]}]; + 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, 'resid_handle', param.resid_handle); + 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);