Skip to content

Commit

Permalink
Add support for function handles returning residal arrays
Browse files Browse the repository at this point in the history 
In ndbase.lm4 and simplex optimisers. Also added tests
This allows support for ND fitting problems.
  • Loading branch information
@RichardWaiteSTFC
RichardWaiteSTFC committed Nov 4, 2024
1 parent dd978e2 commit 44a6290
Show file tree
Hide file tree
Showing 5 changed files with 71 additions and 15 deletions.
22 changes: 22 additions & 0 deletions +sw_tests/+unit_tests/unittest_ndbase_cost_function_wrapper.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,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
9 changes: 9 additions & 0 deletions +sw_tests/+unit_tests/unittest_ndbase_optimisers.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,15 @@ function test_optimise_data_struct(testCase, optimiser, poly_func)
testCase.verify_val(cost_val, 0, 'abs_tol', 1e-6);
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', 1e-3);
testCase.verify_val(cost_val, 0, 'abs_tol', 1e-6);
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);
Expand Down
17 changes: 14 additions & 3 deletions swfiles/+ndbase/cost_function_wrapper.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,11 @@
% `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
Expand All @@ -76,14 +81,20 @@
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;
obj.calc_resid = [];
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
if ~isfield(options.data,'e') || isempty(options.data.e) || ~any(options.data.e(:))
Expand Down
28 changes: 17 additions & 11 deletions swfiles/+ndbase/lm4.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,9 @@
% Options can be given using the modified output of optimset() or as option
% name string option value pairs.
%
% resid_handle Boolean scalar - if true and `dat` is empty then function
% function handle returns array of residuals, if false (default)
% then function handle returns scalar cost function.
% dp Vector with N or 1 element, defines the fractional increment of
% 'p' when calculating the Jacobian matrix dFunc(x,p)/dp:
% dp(j)>0 central differences calculated,
Expand Down Expand Up @@ -94,46 +97,49 @@
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'}];
inpForm.defval = [inpForm.defval {10 0.3}];
inpForm.size = [inpForm.size {[1 1] [1 1]}];
inpForm.fname = [inpForm.fname {'nu_up', 'nu_dn', 'resid_handle'}];
inpForm.defval = [inpForm.defval {10 0.3, false}];
inpForm.size = [inpForm.size {[1 1] [1 1], [1 1]}];

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);
% transform starting values into their unconstrained surrogates.
p0_free = cost_func_wrap.get_free_parameters(p0);

if isempty(dat)
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;
ndof = 1;
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;
ndof = numel(dat.x) - cost_func_wrap.get_num_free_parameters() + 1;
% 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)
% All parameters fixed, evaluate cost at initial guess
% don't use p0 as could contain fixed params outside bounds
% Re-calc bound params as p0 as could have fixed params outside bounds
pOpt = cost_func_wrap.get_bound_parameters(p0_free);
fVal = cost_func_wrap.eval_cost_function(p0_free)/ndof;
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;
[cost_val, resids] = eval_cost_func(cost_func_wrap, p);
[hess, jac] = calc_hessian_and_jacobian_scalar(cost_func_wrap, p, diff_step, cost_val, resids);
exit_flag = 0;
for niter = 1:param.MaxIter
Expand Down
10 changes: 9 additions & 1 deletion swfiles/+ndbase/simplex.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,10 @@
% the weighted least square deviation from data). Default value is
% $10^{-3}$.
%
% `'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.
%
% ### See Also
%
Expand All @@ -149,10 +153,14 @@
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'}];
inpForm.defval = [inpForm.defval {false}];
inpForm.size = [inpForm.size {[1 1]}];

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);

% transform starting values into their unconstrained surrogates.
p0_free = cost_func_wrap.get_free_parameters(p0);
Expand Down

0 comments on commit 44a6290

Please sign in to comment.