1% STK_PARAM_ESTIM_OPTIMIZE [STK internal]
2%
3% INTERNAL FUNCTION WARNING:
4%
5%    This function is currently considered as internal.  API-breaking
6%    changes are very likely to happen in future releases.
7
8% Copyright Notice
9%
10%    Copyright (C) 2015-2019 CentraleSupelec
11%    Copyright (C) 2014 Ashwin Ravisankar
12%    Copyright (C) 2011-2014 SUPELEC
13%
14%    Authors:  Julien Bect        <julien.bect@centralesupelec.fr>
15%              Emmanuel Vazquez   <emmanuel.vazquez@centralesupelec.fr>
16%              Ashwin Ravisankar  <ashwinr1993@gmail.com>
17
18% Copying Permission Statement
19%
20%    This file is part of
21%
22%            STK: a Small (Matlab/Octave) Toolbox for Kriging
23%               (http://sourceforge.net/projects/kriging)
24%
25%    STK is free software: you can redistribute it and/or modify it under
26%    the terms of the GNU General Public License as published by the Free
27%    Software Foundation,  either version 3  of the License, or  (at your
28%    option) any later version.
29%
30%    STK is distributed  in the hope that it will  be useful, but WITHOUT
31%    ANY WARRANTY;  without even the implied  warranty of MERCHANTABILITY
32%    or FITNESS  FOR A  PARTICULAR PURPOSE.  See  the GNU  General Public
33%    License for more details.
34%
35%    You should  have received a copy  of the GNU  General Public License
36%    along with STK.  If not, see <http://www.gnu.org/licenses/>.
37
38function [model_opt, info] = stk_param_estim_optim ...
39    (model0, xi, zi, criterion, covparam_select, noiseparam_select)
40
41select = [covparam_select; noiseparam_select];
42
43% Starting point
44v0 = stk_get_optimizable_parameters (model0);
45w0 = v0(select);
46
47% Bounds
48% FIXME: this could (should) be implemented directly for models
49[covparam_lb, covparam_ub] = stk_param_getdefaultbounds (model0.covariance_type, model0.param, xi, zi);
50[covparam_lb, covparam_ub] = select_bounds (covparam_lb, covparam_ub, covparam_select);
51[noiseparam_lb, noiseparam_ub] = stk_param_getdefaultbounds_lnv (model0, model0.lognoisevariance, xi, zi);
52[noiseparam_lb, noiseparam_ub] = select_bounds (noiseparam_lb, noiseparam_ub, noiseparam_select);
53lb = [covparam_lb; noiseparam_lb];
54ub = [covparam_ub; noiseparam_ub];
55
56% Define objective function
57f = @(v)(crit_wrapper (model0, v, xi, zi, criterion, covparam_select, noiseparam_select));
58
59bounds_available = (~ isempty (lb)) && (~ isempty (ub));
60
61% Sanity check (part 1)
62crit0 = f (w0);
63if ~ (isscalar (crit0) && isfinite (crit0))
64    errmsg = '*** PANIC: crit0 is not a finite scalar value. ***';
65    stk_error (errmsg, 'OptimizationFailure');
66end
67
68if bounds_available
69    A = stk_options_get ('stk_param_estim', 'minimize_box');
70    [w_opt, crit_opt] = stk_minimize_boxconstrained (A, f, w0, lb, ub);
71else
72    A = stk_options_get ('stk_param_estim', 'minimize_unc');
73    [w_opt, crit_opt] = stk_minimize_unconstrained (A, f, w0);
74end
75
76% Sanity check (part 2)
77if crit0 < crit_opt
78    s1 = '*** PANIC: Something went SERIOUSLY WRONG during the optimization ***';
79    s2 = sprintf ('crit0 = %f,  crit_opt = %f:  crit0 < crit_opt', crit0, crit_opt);
80    stk_error (sprintf ('%s\n%s\n', s1, s2), 'OptimizationFailure');
81end
82
83% Create outputs
84v_opt = v0;
85v_opt(select) = w_opt;
86model_opt = stk_set_optimizable_parameters (model0, v_opt);
87
88% Create 'info' structure, if requested
89if nargout > 1
90    info.criterion = criterion;
91    info.crit_opt = crit_opt;
92    info.starting_point = w0;
93    info.final_point = w_opt;
94    info.lower_bounds = lb;
95    info.upper_bounds = ub;
96    info.param_select = covparam_select;
97    info.noiseparam_select = noiseparam_select;
98end
99
100end % function
101
102%#ok<*CTCH,*LERR,*SPWRN,*WNTAG>
103
104
105function [C, dC] = crit_wrapper ...
106    (model, w, xi, zi, criterion, covparam_select, noise_select)
107
108v = stk_get_optimizable_parameters (model);
109v([covparam_select; noise_select]) = w;
110model = stk_set_optimizable_parameters (model, v);
111
112if nargout == 1
113
114    % Compute only the value of the criterion
115    C = criterion (model, xi, zi);
116
117elseif any (noise_select)
118
119    % Compute the value of the criterion and the gradients
120    % FIXME: We might be computing a lot of derivatives that we don't really need...
121    [C, dC_param, dC_lnv] = criterion (model, xi, zi);
122
123    dC = [dC_param(covparam_select); dC_lnv(noise_select)];
124
125else
126
127    % Compute the value of the criterion and the gradients
128    % FIXME: We might be computing a lot of derivatives that we don't really need...
129    [C, dC_param] = criterion (model, xi, zi);
130
131    dC = dC_param(covparam_select);
132
133end
134
135end % function
136
137
138function [lb, ub] = select_bounds (lb, ub, select)
139
140if ~ isempty (lb)
141    lb = lb(select);
142end
143
144if ~ isempty (ub)
145    ub = ub(select);
146end
147
148end % function
149