octave_packages/optim-1.2.0/optim_problems.m

   1 %% Copyright (C) 2007 Paul Kienzle (sort-based lookup in ODE solver)
   2 %% Copyright (C) 2009 Thomas Treichl <thomas.treichl@gmx.net> (ode23 code)
   3 %% Copyright (C) 2010 Olaf Till <i7tiol@t-online.de>
   4 %%
   5 %% This program is free software; you can redistribute it and/or modify it under
   6 %% the terms of the GNU General Public License as published by the Free Software
   7 %% Foundation; either version 3 of the License, or (at your option) any later
   8 %% version.
   9 %%
  10 %% This program is distributed in the hope that it will be useful, but WITHOUT
  11 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 %% details.
  14 %%
  15 %% You should have received a copy of the GNU General Public License along with
  16 %% this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 %% Problems for testing optimizers. Documentation is in the code.
  19
  20 function ret = optim_problems ()
  21
  22   %% As little external code as possible is called. This leads to some
  23   %% duplication of external code. The advantages are that thus these
  24   %% problems do not change with evolving external code, and that
  25   %% optimization results in Octave can be compared with those in Matlab
  26   %% without influence of differences in external code (e.g. ODE
  27   %% solvers). Even calling 'interp1 (..., ..., ..., 'linear')' is
  28   %% avoided by using an internal subfunction, although this is possibly
  29   %% too cautious.
  30   %%
  31   %% For cross-program comparison of optimizers, the code of these
  32   %% problems is intended to be Matlab compatible.
  33   %%
  34   %% External data may be loaded, which should be supplied in the
  35   %% 'private/' subdirectory. Use the variable 'ddir', which contains
  36   %% the path to this directory.
  37
  38   %% Note: The difficulty of problems with dynamic models often
  39   %% decisively depends on the the accuracy of the used ODE(DAE)-solver.
  40
  41   %% Description of the returned structure
  42   %%
  43   %% According to 3 classes of problems, there are (or should be) three
  44   %% fields: 'curve' (curve fitting), 'general' (general optimization),
  45   %% and 'zero' (zero finding). The subfields are labels for the
  46   %% particular problems.
  47   %%
  48   %% Under the label fields, there are subfields mostly identical
  49   %% between the 3 classes of problems (some may contain empty values):
  50   %%
  51   %% .f: handle of an internally defined objective function (argument:
  52   %% column vector of parameters), meant for minimization, or to a
  53   %% 'model function' (arguments: independents, column vector of
  54   %% parameters) in the case of curve fitting, where .f should return a
  55   %% matrix of equal dimensions as .data.y below.
  56   %%
  57   %% .dfdp: handle of internally defined function for jacobian of
  58   %% objective function or 'model function', respectively.
  59   %%
  60   %% .init_p: initial parameters, column vector
  61   %%
  62   %% possibly .init_p_b: two column matrix of ranges to choose initial
  63   %% parameters from
  64   %%
  65   %% possibly .init_p_f: handle of internally defined function which
  66   %% returns a column vector of initial parameters unique to the index
  67   %% given as function argument; given '0' as function argument,
  68   %% .init_p_f returns the maximum index
  69   %%
  70   %% .result.p: parameters of best known result
  71   %%
  72   %% possibly .result.obj: value of objective function for .result.p (or
  73   %% sum of squared residuals in curve fitting).
  74   %%
  75   %% .data.x: matrix of independents (curve fitting)
  76   %%
  77   %% .data.y: matrix of observations, dimensions may be independent of
  78   %% .data.x (curve fitting)
  79   %%
  80   %% .data.wt: matrix of weights, same dimensions as .data.y (curve
  81   %% fitting)
  82   %%
  83   %% .data.cov: covariance matrix of .data.y(:) (not necessarily a
  84   %% diagonal matrix, which could be expressed in .data.wt)
  85   %%
  86   %% .strict_inequc, .non_strict_inequc, .equc: 'strict' inequality
  87   %% constraints (<, >), 'non-strict' inequality constraints (<=, >=),
  88   %% and equality constraints, respectively. Subfields are: .bounds
  89   %% (except in equality constraints): two-column matrix of ranges;
  90   %% .linear: cell-array {m, v}, meaning linear constraints m.' *
  91   %% parameters + v >|>=|== 0; .general: handle of internally defined
  92   %% function h with h (p) >|>=|== 0; possibly .general_dcdp: handle of
  93   %% internally defined function (argument: parameters) returning the
  94   %% jacobian of the constraints given in .general. For the sake of
  95   %% optimizers which can exploit this, the function in subfield
  96   %% .general may accept a logical index vector as an optional second
  97   %% argument, returning only the indexed constraint values.
  98
  99
 100   %% Please keep the following list of problems current.
 101   %%
 102   %% .curve.p_1, .curve.p_2, .curve.p_3_d2: from 'Comparison of gradient
 103   %% methods for the solution of nonlinear parameter estimation
 104   %% problems' (1970), Yonathan Bard, Siam Journal on Numerical Analysis
 105   %% 7(1), 157--186. The numbering of problems is the same as in the
 106   %% article. Since Bard used strict bounds, testing optimizers which
 107   %% used penalization for bounds, the bounds are changed here to allow
 108   %% testing with non-strict bounds (<= or >=). .curve.p_3_d2 involves
 109   %% dynamic modeling. These are not necessarily difficult problems.
 110   %%
 111   %% .curve.p_3_d2_noweights: problem .curve.p_3_d2 equivalently
 112   %% re-formulated without weights.
 113   %%
 114   %% .curve.p_r: A seemingly more difficult 'real life' problem with
 115   %% dynamic modeling. To assess optimizers, .init_p_f should be used
 116   %% with 1:64. There should be two groups of results, indicating the
 117   %% presence of two local minima. Olaf Till <olaf.till@uni-jena.de>
 118   %%
 119   %% .....schittkowski_...: Klaus Schittkowski: 'More test examples for
 120   %% nonlinear programming codes.' Lecture Notes in Economics and
 121   %% Mathematical Systems 282, Berlin 1987. The published problems are
 122   %% numbered from 201 to 395 and may appear here under the fields
 123   %% .curve, .general, or .zero.
 124   %%
 125   %% .general.schittkowski_281: 10 parameters, unconstrained.
 126   %%
 127   %% .general.schittkowski_289: 30 parameters, unconstrained.
 128   %%
 129   %% .general.schittkowski_327 and
 130   %%
 131   %% .curve.schittkowski_327: Two parameters, one general inequality
 132   %% constraint, two bounds. The best solution given in the publication
 133   %% seems not very good (it probably has been achieved with general
 134   %% minimization, not curve fitting) and has been replaced here by a
 135   %% better (leasqr).
 136   %%
 137   %% .curve.schittkowski_372 and
 138   %%
 139   %% .general.schittkowski_372: 9 parameters, 12 general inequality
 140   %% constraints, 6 bounds. Infeasible initial parameters
 141   %% (.curve.schittkowski_372.init_p_f(1) provides a set of more or less
 142   %% feasible parameters). leasqr sticks at the (feasible) initial
 143   %% values. sqp has no problems.
 144   %%
 145   %% .curve.schittkowski_373: 9 parameters, 6 equality constraints.
 146   %% Infeasible initial parameters (.curve.schittkowski_373.init_p_f(1)
 147   %% provides a set of more or less feasible parameters). leasqr sticks
 148   %% at the (feasible) initial values. sqp has no problems.
 149   %%
 150   %% .general.schittkowski_391: 30 parameters, unconstrained. The best
 151   %% solution given in the publication seems not very good, obviously
 152   %% the used routine had not managed to get very far from the starting
 153   %% parameters; it has been replaced here by a better (Octaves
 154   %% fminunc). The result still varies widely (without much changes in
 155   %% objective function) with changes of starting values. Maybe not a
 156   %% very good test problem, no well defined minimum ...
 157
 158   %% needed for some anonymous functions
 159   if (exist ('ifelse') ~= 5)
 160     ifelse = @ scalar_ifelse;
 161   end
 162
 163   if (~exist ('OCTAVE_VERSION'))
 164     NA = NaN;
 165   end
 166
 167   %% determine the directory of this functions file
 168   fdir = fileparts (mfilename ('fullpath'));
 169   %% data directory
 170   ddir = sprintf ('%s%sprivate%s', fdir, filesep, filesep);
 171
 172   ret.curve.p_1.dfdp = [];
 173   ret.curve.p_1.init_p = [1; 1; 1];
 174   ret.curve.p_1.data.x = cat (2, ...
 175                               (1:15).', ...
 176                               (15:-1:1).', ...
 177                               [(1:8).'; (7:-1:1).']);
 178   ret.curve.p_1.data.y = [.14; .18; .22; .25; .29; .32; .35; .39; ...
 179                           .37; .58; .73; .96; 1.34; 2.10; 4.39];
 180   ret.curve.p_1.data.wt = [];
 181   ret.curve.p_1.data.cov = [];
 182   ret.curve.p_1.result.p = [.08241040; 1.133033; 2.343697];
 183   ret.curve.p_1.strict_inequc.bounds = [0, 100; 0, 100; 0, 100];
 184   ret.curve.p_1.strict_inequc.linear = [];
 185   ret.curve.p_1.strict_inequc.general = [];
 186   ret.curve.p_1.non_strict_inequc.bounds = ...
 187       [eps, 100; eps, 100; eps, 100];
 188   ret.curve.p_1.non_strict_inequc.linear = [];
 189   ret.curve.p_1.non_strict_inequc.general = [];
 190   ret.curve.p_1.equc.linear = [];
 191   ret.curve.p_1.equc.general = [];
 192   ret.curve.p_1.f = @ f_1;
 193
 194   ret.curve.p_2.dfdp = [];
 195   ret.curve.p_2.init_p = [0; 0; 0; 0; 0];
 196   ret.curve.p_2.data.x = [.871, .643, .550; ...
 197                           .228, .669, .854; ...
 198                           .528, .229, .170; ...
 199                           .110, .354, .337; ...
 200                           .911, .056, .493; ...
 201                           .476, .154, .918; ...
 202                           .655, .421, .077; ...
 203                           .649, .140, .199; ...
 204                           .995, .045,   NA; ...
 205                           .130, .016, .195; ...
 206                           .823, .690, .690; ...
 207                           .768, .992, .389; ...
 208                           .203, .740, .120; ...
 209                           .302, .519, .221; ...
 210                           .991, .450, .249; ...
 211                           .224, .030, .502; ...
 212                           .428, .127, .772; ...
 213                           .552, .494, .110; ...
 214                           .461, .824, .714; ...
 215                           .799, .494, .295];
 216   ret.curve.p_2.data.y = zeros (20, 3);
 217   ret.curve.p_2.data.wt = [];
 218   ret.curve.p_2.data.cov = [];
 219   ret.curve.p_2.data.misc = [4.36, 5.21, 5.35; ...
 220                              4.99, 3.30, 3.10; ...
 221                              1.67,   NA, 2.75; ...
 222                              2.17, 1.48, 1.49; ...
 223                              2.98, 4.69, 4.23; ...
 224                              4.46, 3.87, 3.15; ...
 225                              1.79, 3.18, 3.57; ...
 226                              1.71, 3.13, 3.07; ...
 227                              3.07, 5.01, 4.58; ...
 228                              0.94, 0.93, 0.74; ...
 229                              4.97, 5.37, 5.35; ...
 230                              4.32, 4.85, 5.46; ...
 231                              2.17, 1.78, 2.43; ...
 232                              2.22, 2.18, 2.44; ...
 233                              2.88, 4.90, 5.11; ...
 234                              2.29, 1.94, 1.46; ...
 235                              3.76, 3.39, 2.71; ...
 236                              1.99, 2.93, 3.31; ...
 237                              4.95, 4.08, 4.19; ...
 238                              2.96, 4.26, 4.48];
 239   ret.curve.p_2.result.p = [.9925145; 2.005293; 3.999732; ...
 240                             2.680371; .4977683]; % from maximum
 241                                 % likelyhood optimization
 242   ret.curve.p_2.strict_inequc.bounds = [];
 243   ret.curve.p_2.strict_inequc.linear = [];
 244   ret.curve.p_2.strict_inequc.general = [];
 245   ret.curve.p_2.non_strict_inequc.bounds = [];
 246   ret.curve.p_2.non_strict_inequc.linear = [];
 247   ret.curve.p_2.non_strict_inequc.general = [];
 248   ret.curve.p_2.equc.linear = [];
 249   ret.curve.p_2.equc.general = [];
 250   ret.curve.p_2.f = @ (x, p) f_2 (x, p, ret.curve.p_2.data.misc);
 251
 252
 253
 254   ret.curve.p_3_d2.dfdp = [];
 255   ret.curve.p_3_d2.init_p = [.01; .01; .001; .001; .02; .001];
 256   ret.curve.p_3_d2.data.x = [0; 12.5; 25; 37.5; 50; ...
 257                              62.5; 75; 87.5; 100];
 258   ret.curve.p_3_d2.data.y=[1       1       0       0       0      ; ...
 259                            .945757 .961201 .494861 .154976 .111485; ...
 260                            .926486 .928762 .690492 .314501 .236263; ...
 261                            .917668 .915966 .751806 .709300 .311747; ...
 262                            .928987 .917542 .771559 1.19224 .333096; ...
 263                            .927782 .920075 .780903 1.68815 .340324; ...
 264                            .925304 .912330 .790539 2.19539 .356787; ...
 265                            .925083 .917684 .783933 2.74211 .358283; ...
 266                            .917277 .907529 .779259 3.20025 .361969];
 267   ret.curve.p_3_d2.data.y(:, 3) = ...
 268       ret.curve.p_3_d2.data.y(:, 3) / 10;
 269   ret.curve.p_3_d2.data.y(:, 4:5) = ...
 270       ret.curve.p_3_d2.data.y(:, 4:5) / 1000;
 271   ret.curve.p_3_d2.data.wt = repmat ([.1, .1, 1, 10, 100], 9, 1);
 272   ret.curve.p_3_d2.data.cov = [];
 273   ret.curve.p_3_d2.result.p = [.6358247e-2; ...
 274                                .6774551e-1; ...
 275                                .5914274e-4; ...
 276                                .4944010e-3; ...
 277                                .1018828; ...
 278                                .4210526e-3];
 279   ret.curve.p_3_d2.strict_inequc.bounds = [0, 1; ...
 280                                            0, 1; ...
 281                                            0, .1; ...
 282                                            0, .1; ...
 283                                            0, 2; ...
 284                                            0, .1];
 285   ret.curve.p_3_d2.strict_inequc.linear = [];
 286   ret.curve.p_3_d2.strict_inequc.general = [];
 287   ret.curve.p_3_d2.non_strict_inequc.bounds = [eps, 1; ...
 288                                                eps, 1; ...
 289                                                eps, .1; ...
 290                                                eps, .1; ...
 291                                                eps, 2; ...
 292                                                eps, .1];
 293   ret.curve.p_3_d2.non_strict_inequc.linear = [];
 294   ret.curve.p_3_d2.non_strict_inequc.general = [];
 295   ret.curve.p_3_d2.equc.linear = [];
 296   ret.curve.p_3_d2.equc.general = [];
 297   ret.curve.p_3_d2.f = @ f_3;
 298
 299   ret.curve.p_3_d2_noweights = ret.curve.p_3_d2;
 300   ret.curve.p_3_d2_noweights.data.wt = [];
 301   ret.curve.p_3_d2_noweights.data.y(:, 1:2) = ...
 302       ret.curve.p_3_d2_noweights.data.y(:, 1:2) * .1;
 303   ret.curve.p_3_d2_noweights.data.y(:, 4) = ...
 304       ret.curve.p_3_d2_noweights.data.y(:, 4) * 10;
 305   ret.curve.p_3_d2_noweights.data.y(:, 5) = ...
 306       ret.curve.p_3_d2_noweights.data.y(:, 5) * 100;
 307   ret.curve.p_3_d2_noweights.f = @ f_3_noweights;
 308
 309   ret.curve.p_r.dfdp = [];
 310   ret.curve.p_r.init_p = [.3; .03; .003; .7; 1000; .0205];
 311   ret.curve.p_r.init_p_b = [.3, .5; ...
 312                             .03, .05; ...
 313                             .003, .005; ...
 314                             .7, .9; ...
 315                             1000, 1300; ...
 316                             .0205, .023];
 317   ret.curve.p_r.init_p_f = @ (id) pc2 (ret.curve.p_r.init_p_b, id);
 318   hook.ns = [84; 84; 85; 86; 84; 84; 84; 84];
 319   xb = [0.2, 0.8640; ...
 320         0.2, 0.5320; ...
 321         0.2, 0.4856; ...
 322         0.2, 0.4210; ...
 323         0.2, 0.3328; ...
 324         0.2, 0.2996; ...
 325         0.2, 0.2664; ...
 326         0.2, 0.2498];
 327   ns = cat (1, 0, cumsum (hook.ns));
 328   x = zeros (ns(end), 1);
 329   for id = 1:8
 330     x(ns(id) + 1 : ns(id + 1)) = ...
 331         linspace (xb(id, 1), xb(id, 2), hook.ns(id)).';
 332   end
 333   hook.t = x;
 334   ret.curve.p_r.data.x = x;
 335   ret.curve.p_r.data.y = ...
 336       load (sprintf ('%soptim_problems_p_r_y.data', ddir));
 337   ret.curve.p_r.data.wt = [];
 338   ret.curve.p_r.data.cov = [];
 339   ret.curve.p_r.result.p = [4.742909e-01; ...
 340                             3.837951e-02; ...
 341                             3.652570e-03; ...
 342                             7.725986e-01; ...
 343                             1.180967e+03; ...
 344                             2.107000e-02];
 345   ret.curve.p_r.result.obj = 0.2043396;
 346   ret.curve.p_r.strict_inequc.bounds = [];
 347   ret.curve.p_r.strict_inequc.linear = [];
 348   ret.curve.p_r.strict_inequc.general = [];
 349   ret.curve.p_r.non_strict_inequc.bounds = [];
 350   ret.curve.p_r.non_strict_inequc.linear = [];
 351   ret.curve.p_r.non_strict_inequc.general = [];
 352   ret.curve.p_r.equc.linear = [];
 353   ret.curve.p_r.equc.general = [];
 354   hook.mc = [2.0019999999999999e-01, 1.9939999999999999e-01, ...
 355              1.9939999999999999e-01, 1.9780000000000000e-01, ...
 356              2.0080000000000001e-01, 1.9960000000000000e-01, ...
 357              1.9960000000000000e-01, 1.9980000000000001e-01; ...
 358              ...
 359              2.0060000000000000e-01, 2.0160000000000000e-01, ...
 360              2.0200000000000001e-01, 2.0200000000000001e-01, ...
 361              2.0180000000000001e-01, 2.0899999999999999e-01, ...
 362              2.0860000000000001e-01, 2.0820000000000000e-01; ...
 363              ...
 364              2.1999144799999999e-02, 2.1998803099999999e-02, ...
 365              2.2000449599999999e-02, 2.2000024399999998e-02, ...
 366              2.1998160999999999e-02, 2.1999289000000002e-02, ...
 367              2.1998038800000001e-02, 2.2000270999999998e-02; ...
 368              ...
 369              -6.8806551999999986e-03, -1.3768898999999999e-02, ...
 370              -1.6065479000000001e-02, -2.0657919600000001e-02, ...
 371              -3.4479971099999999e-02, -4.5934394099999998e-02, ...
 372              -6.9011619100000005e-02, -9.1971348400000000e-02; ...
 373              ...
 374              2.3383865100000002e-02, 2.4768462500000001e-02, ...
 375              2.5231915899999999e-02, 2.6155515300000001e-02, ...
 376              2.8933514200000000e-02, 3.1235568599999999e-02, ...
 377              3.5874086299999997e-02, 4.0490560699999997e-02; ...
 378              ...
 379              -1.8240616806039459e+05, -1.6895474269973661e+03, ...
 380              -8.1072652464694931e+02, -7.0113302985566395e+02, ...
 381              1.0929964862867249e+04, 3.5665776039585688e+02, ...
 382              5.7400262910547769e+02, 9.1737316974342252e+02; ...
 383              ...
 384              1.0965398741890911e+05, 1.0131334821116490e+03, ...
 385              4.8504892529762208e+02, 4.1801020186158411e+02, ...
 386              -6.6178457662355086e+03, -2.2103886018172699e+02, ...
 387              -3.5529578864017282e+02, -5.6690686490678263e+02; ...
 388              ...
 389              -2.1972917026209168e+04, -2.0250659086265861e+02, ...
 390              -9.6733175964156985e+01, -8.3069683020988421e+01, ...
 391              1.3356173243752210e+03, 4.5610806266307627e+01, ...
 392              7.3229009073208331e+01, 1.1667126232349770e+02; ...
 393              ...
 394              1.4676952576063929e+03, 1.3514357622838521e+01, ...
 395              6.4524906786197480e+00, 5.5245948033669476e+00, ...
 396              -8.9827382090060922e+01, -3.1118708128841241e+00, ...
 397              -5.0039950796246986e+00, -7.9749636293721071e+00];
 398   ret.curve.p_r.f = @ (x, p) f_r (x, p, hook);
 399
 400   ret.general.schittkowski_281.dfdp = ...
 401       @ (p) schittkowski_281_dfdp (p);
 402   ret.general.schittkowski_281.init_p = zeros (10, 1);
 403   ret.general.schittkowski_281.result.p = ones (10, 1); % 'theoretically'
 404   ret.general.schittkowski_281.result.obj = 0; % 'theoretically'
 405   ret.general.schittkowski_281.strict_inequc.bounds = [];
 406   ret.general.schittkowski_281.strict_inequc.linear = [];
 407   ret.general.schittkowski_281.strict_inequc.general = [];
 408   ret.general.schittkowski_281.non_strict_inequc.bounds = [];
 409   ret.general.schittkowski_281.non_strict_inequc.linear = [];
 410   ret.general.schittkowski_281.non_strict_inequc.general = [];
 411   ret.general.schittkowski_281.equc.linear = [];
 412   ret.general.schittkowski_281.equc.general = [];
 413   ret.general.schittkowski_281.f = ...
 414       @ (p) (sum (((1:10).') .^ 3 .* (p - 1) .^ 2)) ^ (1 / 3);
 415
 416   ret.general.schittkowski_289.dfdp = ...
 417       @ (p) exp (- sum (p .^ 2) / 60) / 30 * p;
 418   ret.general.schittkowski_289.init_p = [-1.03; 1.07; -1.10; 1.13; ...
 419                                          -1.17; 1.20; -1.23; 1.27; ...
 420                                          -1.30; 1.33; -1.37; 1.40; ...
 421                                          -1.43; 1.47; -1.50; 1.53; ...
 422                                          -1.57; 1.60; -1.63; 1.67; ...
 423                                          -1.70; 1.73; -1.77; 1.80; ...
 424                                          -1.83; 1.87; -1.90; 1.93; ...
 425                                          -1.97; 2.00];
 426   ret.general.schittkowski_289.result.p = zeros (30, 1); % 'theoretically'
 427   ret.general.schittkowski_289.result.obj = 0; % 'theoretically'
 428   ret.general.schittkowski_289.strict_inequc.bounds = [];
 429   ret.general.schittkowski_289.strict_inequc.linear = [];
 430   ret.general.schittkowski_289.strict_inequc.general = [];
 431   ret.general.schittkowski_289.non_strict_inequc.bounds = [];
 432   ret.general.schittkowski_289.non_strict_inequc.linear = [];
 433   ret.general.schittkowski_289.non_strict_inequc.general = [];
 434   ret.general.schittkowski_289.equc.linear = [];
 435   ret.general.schittkowski_289.equc.general = [];
 436   ret.general.schittkowski_289.f = @ (p) 1 - exp (- sum (p .^ 2) / 60);
 437
 438   ret.curve.schittkowski_327.dfdp = ...
 439       @ (x, p) [1 + exp(-p(2) * (x - 8)), ...
 440                 (p(1) + .49) * (8 - x) .* exp (-p(2) * (x - 8))];
 441   ret.curve.schittkowski_327.init_p = [.42; 5];
 442   ret.curve.schittkowski_327.data.x = ...
 443       [8; 8; 10; 10; 10; 10; 12; 12; 12; 12; 14; 14; 14; 16; 16; 16; ...
 444        18; 18; 20; 20; 20; 22; 22; 22; 24; 24; 24; 26; 26; 26; 28; ...
 445        28; 30; 30; 30; 32; 32; 34; 36; 36; 38; 38; 40; 42];
 446   ret.curve.schittkowski_327.data.y= ...
 447       [.49; .49; .48; .47; .48; .47; .46; .46; .45; .43; .45; .43; ...
 448        .43; .44; .43; .43; .46; .45; .42; .42; .43; .41; .41; .40; ...
 449        .42; .40; .40; .41; .40; .41; .41; .40; .40; .40; .38; .41; ...
 450        .40; .40; .41; .38; .40; .40; .39; .39];
 451   ret.curve.schittkowski_327.data.wt = [];
 452   ret.curve.schittkowski_327.data.cov = [];
 453   %% This result was given by Schittkowski. No constraint is active
 454   %% here. The second parameter is unchanged from initial value.
 455   %%
 456   %% ret.curve.schittkowski_327.result.p = [.4219; 5];
 457   %% ret.curve.schittkowski_327.result.obj = .0307986;
 458   %%
 459   %% This is the result of leasqr of Octave Forge. The general
 460   %% constraint is active here. Both parameters are different from
 461   %% initial value. The value of the objective function is better.
 462   %%
 463   ret.curve.schittkowski_327.result.p = [.4199227; 1.2842958];
 464   ret.curve.schittkowski_327.result.obj = .0284597;
 465   ret.curve.schittkowski_327.strict_inequc.bounds = [];
 466   ret.curve.schittkowski_327.strict_inequc.linear = [];
 467   ret.curve.schittkowski_327.strict_inequc.general = [];
 468   ret.curve.schittkowski_327.non_strict_inequc.bounds = [.4, Inf; ...
 469                                                          .4, Inf];
 470   ret.curve.schittkowski_327.non_strict_inequc.linear = [];
 471   ret.curve.schittkowski_327.non_strict_inequc.general = ...
 472       @ (p, varargin) apply_idx_if_given ...
 473       (-.09 - p(1) * p(2) + .49 * p(2), varargin{:});
 474   ret.curve.schittkowski_327.equc.linear = [];
 475   ret.curve.schittkowski_327.equc.general = [];
 476   ret.curve.schittkowski_327.f = ...
 477       @ (x, p) p(1) + (.49 - p(1)) * exp (-p(2) * (x - 8));
 478
 479   ret.general.schittkowski_327.init_p = [.42; 5];
 480   ret.general.schittkowski_327.data.x = ...
 481       [8; 8; 10; 10; 10; 10; 12; 12; 12; 12; 14; 14; 14; 16; 16; 16; ...
 482        18; 18; 20; 20; 20; 22; 22; 22; 24; 24; 24; 26; 26; 26; 28; ...
 483        28; 30; 30; 30; 32; 32; 34; 36; 36; 38; 38; 40; 42];
 484   ret.general.schittkowski_327.data.y= ...
 485       [.49; .49; .48; .47; .48; .47; .46; .46; .45; .43; .45; .43; ...
 486        .43; .44; .43; .43; .46; .45; .42; .42; .43; .41; .41; .40; ...
 487        .42; .40; .40; .41; .40; .41; .41; .40; .40; .40; .38; .41; ...
 488        .40; .40; .41; .38; .40; .40; .39; .39];
 489   x = ret.general.schittkowski_327.data.x;
 490   y = ret.general.schittkowski_327.data.y;
 491   ret.general.schittkowski_327.dfdp = ...
 492       @ (p) cat (2, ...
 493                  2 * sum ((exp (-p(2 * x - 8)) - 1) * ...
 494                           (y + (p(1) - .49) * ...
 495                            exp (-p(2) * (x - 8)) - p1)), ...
 496                  2 * (p(1) - .49) * ...
 497                  sum ((8 - x) * exp (-p(2 * x - 8)) * ...
 498                       (y + (p(1) - .49) * ...
 499                        exp (-p(2) * (x - 8)) - p1)));
 500   %% This result was given by Schittkowski. No constraint is active
 501   %% here. The second parameter is unchanged from initial value.
 502   %%
 503   %% ret.general.schittkowski_327.result.p = [.4219; 5];
 504   %% ret.general.schittkowski_327.result.obj = .0307986;
 505   %%
 506   %% This is the result of leasqr of Octave Forge. The general
 507   %% constraint is active here. Both parameters are different from
 508   %% initial value. The value of the objective function is better. sqp
 509   %% gives a similar result.
 510   ret.general.schittkowski_327.result.p = [.4199227; 1.2842958];
 511   ret.general.schittkowski_327.result.obj = .0284597;
 512   ret.general.schittkowski_327.strict_inequc.bounds = [];
 513   ret.general.schittkowski_327.strict_inequc.linear = [];
 514   ret.general.schittkowski_327.strict_inequc.general = [];
 515   ret.general.schittkowski_327.non_strict_inequc.bounds = [.4, Inf; ...
 516                                                          .4, Inf];
 517   ret.general.schittkowski_327.non_strict_inequc.linear = [];
 518   ret.general.schittkowski_327.non_strict_inequc.general = ...
 519       @ (p, varargin) apply_idx_if_given ...
 520       (-.09 - p(1) * p(2) + .49 * p(2), varargin{:});
 521   ret.general.schittkowski_327.equc.linear = [];
 522   ret.general.schittkowski_327.equc.general = [];
 523   ret.general.schittkowski_327.f = ...
 524       @ (p) sumsq (y - p(1) - (.49 - p(1)) * exp (-p(2) * (x - 8)));
 525
 526   ret.curve.schittkowski_372.dfdp = ...
 527       @ (x, p) cat (2, zeros (6, 3), eye (6));
 528   %% given by Schittkowski, not feasible
 529   ret.curve.schittkowski_372.init_p = [300; -100; -.1997; -127; ...
 530                                        -151; 379; 421; 460; 426];
 531   %% computed with sqp and a constant objective function, (almost)
 532   %% feasible
 533   ret.curve.schittkowski_372.init_p_f = @ (id) ...
 534       ifelse (id == 0, 1, [2.951277e+02; ...
 535                            -1.058720e+02; ...
 536                            -9.535824e-02; ...
 537                            2.421108e+00; ...
 538                            3.191822e+00; ...
 539                            3.790000e+02; ...
 540                            4.210000e+02; ...
 541                            4.600000e+02; ...
 542                            4.260000e+02]);
 543   ret.curve.schittkowski_372.data.x = (1:6).'; % any different numbers
 544   ret.curve.schittkowski_372.data.y= zeros (6, 1);
 545   ret.curve.schittkowski_372.data.wt = [];
 546   ret.curve.schittkowski_372.data.cov = [];
 547   %% recomputed with sqp (i.e. not with curve fitting)
 548   ret.curve.schittkowski_372.result.p = [5.2330557804078126e+02; ...
 549                                          -1.5694790476454301e+02; ...
 550                                          -1.9966450018535931e-01; ...
 551                                          2.9607990282984435e+01; ...
 552                                          8.6615541706550545e+01; ...
 553                                          4.7326722338555498e+01; ...
 554                                          2.6235616534580515e+01; ...
 555                                          2.2915996663200740e+01; ...
 556                                          3.9470733973874445e+01];
 557   ret.curve.schittkowski_372.result.obj = 13390.1;
 558   ret.curve.schittkowski_372.strict_inequc.bounds = [];
 559   ret.curve.schittkowski_372.strict_inequc.linear = [];
 560   ret.curve.schittkowski_372.strict_inequc.general = [];
 561   ret.curve.schittkowski_372.non_strict_inequc.bounds = [-Inf, Inf; ...
 562                                                          -Inf, Inf; ...
 563                                                          -Inf, Inf; ...
 564                                                          0, Inf; ...
 565                                                          0, Inf; ...
 566                                                          0, Inf; ...
 567                                                          0, Inf; ...
 568                                                          0, Inf; ...
 569                                                          0, Inf];
 570   ret.curve.schittkowski_372.non_strict_inequc.linear = [];
 571   ret.curve.schittkowski_372.non_strict_inequc.general = ...
 572       @ (p, varargin) apply_idx_if_given ...
 573       (cat (1, p(1) + p(2) * exp (-5 * p(3)) + p(4) - 127, ...
 574             p(1) + p(2) * exp (-3 * p(3)) + p(5) - 151, ...
 575             p(1) + p(2) * exp (-p(3)) + p(6) - 379, ...
 576             p(1) + p(2) * exp (p(3)) + p(7) - 421, ...
 577             p(1) + p(2) * exp (3 * p(3)) + p(8) - 460, ...
 578             p(1) + p(2) * exp (5 * p(3)) + p(9) - 426, ...
 579             -p(1) - p(2) * exp (-5 * p(3)) + p(4) + 127, ...
 580             -p(1) - p(2) * exp (-3 * p(3)) + p(5) + 151, ...
 581             -p(1) - p(2) * exp (-p(3)) + p(6) + 379, ...
 582             -p(1) - p(2) * exp (p(3)) + p(7) + 421, ...
 583             -p(1) - p(2) * exp (3 * p(3)) + p(8) + 460, ...
 584             -p(1) - p(2) * exp (5 * p(3)) + p(9) + 426), ...
 585        varargin{:});
 586   ret.curve.schittkowski_372.equc.linear = [];
 587   ret.curve.schittkowski_372.equc.general = [];
 588   ret.curve.schittkowski_372.f = @ (x, p) p(4:9);
 589
 590   ret.curve.schittkowski_373.dfdp = ...
 591       @ (x, p) cat (2, zeros (6, 3), eye (6));
 592   %% not feasible
 593   ret.curve.schittkowski_373.init_p = [300; -100; -.1997; -127; ...
 594                                        -151; 379; 421; 460; 426];
 595   %% feasible
 596   ret.curve.schittkowski_373.init_p_f = @ (id) ...
 597       ifelse (id == 0, 1, [2.5722721227695763e+02; ...
 598                            -1.5126681606092043e+02; ...
 599                            8.3101871447778766e-02; ...
 600                            -3.0390506000425454e+01; ...
 601                            1.1661334225083069e+01; ...
 602                            2.6097719374430665e+02; ...
 603                            3.2814725183082305e+02; ...
 604                            3.9686840023267564e+02; ...
 605                            3.9796353824451995e+02]);
 606   ret.curve.schittkowski_373.data.x = (1:6).'; % any different numbers
 607   ret.curve.schittkowski_373.data.y= zeros (6, 1);
 608   ret.curve.schittkowski_373.data.wt = [];
 609   ret.curve.schittkowski_373.data.cov = [];
 610   ret.curve.schittkowski_373.result.p = [523.31; ...
 611                                          -156.95; ...
 612                                          -.2; ...
 613                                          29.61; ...
 614                                          -86.62; ...
 615                                          47.33; ...
 616                                          26.24; ...
 617                                          22.92; ...
 618                                          -39.47];
 619   ret.curve.schittkowski_373.result.obj = 13390.1;
 620   ret.curve.schittkowski_373.strict_inequc.bounds = [];
 621   ret.curve.schittkowski_373.strict_inequc.linear = [];
 622   ret.curve.schittkowski_373.strict_inequc.general = [];
 623   ret.curve.schittkowski_373.non_strict_inequc.bounds = [];
 624   ret.curve.schittkowski_373.non_strict_inequc.linear = [];
 625   ret.curve.schittkowski_373.non_strict_inequc.general = [];
 626   ret.curve.schittkowski_373.equc.linear = [];
 627   ret.curve.schittkowski_373.equc.general =  ...
 628       @ (p, varargin) apply_idx_if_given ...
 629       (cat (1, p(1) + p(2) * exp (-5 * p(3)) + p(4) - 127, ...
 630             p(1) + p(2) * exp (-3 * p(3)) + p(5) - 151, ...
 631             p(1) + p(2) * exp (-p(3)) + p(6) - 379, ...
 632             p(1) + p(2) * exp (p(3)) + p(7) - 421, ...
 633             p(1) + p(2) * exp (3 * p(3)) + p(8) - 460, ...
 634             p(1) + p(2) * exp (5 * p(3)) + p(9) - 426), ...
 635        varargin{:});
 636   ret.curve.schittkowski_373.f = @ (x, p) p(4:9);
 637
 638   ret.general.schittkowski_372.dfdp = ...
 639       @ (p) cat (2, zeros (1, 3), 2 * p(4:9));
 640   %% not feasible
 641   ret.general.schittkowski_372.init_p = [300; -100; -.1997; -127; ...
 642                                        -151; 379; 421; 460; 426];
 643   %% recomputed with sqp
 644   ret.general.schittkowski_372.result.p = [5.2330557804078126e+02; ...
 645                                          -1.5694790476454301e+02; ...
 646                                          -1.9966450018535931e-01; ...
 647                                          2.9607990282984435e+01; ...
 648                                          8.6615541706550545e+01; ...
 649                                          4.7326722338555498e+01; ...
 650                                          2.6235616534580515e+01; ...
 651                                          2.2915996663200740e+01; ...
 652                                          3.9470733973874445e+01];
 653   ret.general.schittkowski_372.result.obj = 13390.1;
 654   ret.general.schittkowski_372.strict_inequc.bounds = [];
 655   ret.general.schittkowski_372.strict_inequc.linear = [];
 656   ret.general.schittkowski_372.strict_inequc.general = [];
 657   ret.general.schittkowski_372.non_strict_inequc.bounds = [-Inf, Inf; ...
 658                                                          -Inf, Inf; ...
 659                                                          -Inf, Inf; ...
 660                                                          0, Inf; ...
 661                                                          0, Inf; ...
 662                                                          0, Inf; ...
 663                                                          0, Inf; ...
 664                                                          0, Inf; ...
 665                                                          0, Inf];
 666   ret.general.schittkowski_372.non_strict_inequc.linear = [];
 667   ret.general.schittkowski_372.non_strict_inequc.general = ...
 668       @ (p, varargin) apply_idx_if_given ...
 669       (cat (1, p(1) + p(2) * exp (-5 * p(3)) + p(4) - 127, ...
 670             p(1) + p(2) * exp (-3 * p(3)) + p(5) - 151, ...
 671             p(1) + p(2) * exp (-p(3)) + p(6) - 379, ...
 672             p(1) + p(2) * exp (p(3)) + p(7) - 421, ...
 673             p(1) + p(2) * exp (3 * p(3)) + p(8) - 460, ...
 674             p(1) + p(2) * exp (5 * p(3)) + p(9) - 426, ...
 675             -p(1) - p(2) * exp (-5 * p(3)) + p(4) + 127, ...
 676             -p(1) - p(2) * exp (-3 * p(3)) + p(5) + 151, ...
 677             -p(1) - p(2) * exp (-p(3)) + p(6) + 379, ...
 678             -p(1) - p(2) * exp (p(3)) + p(7) + 421, ...
 679             -p(1) - p(2) * exp (3 * p(3)) + p(8) + 460, ...
 680             -p(1) - p(2) * exp (5 * p(3)) + p(9) + 426), ...
 681        varargin{:});
 682   ret.general.schittkowski_372.equc.linear = [];
 683   ret.general.schittkowski_372.equc.general = [];
 684   ret.general.schittkowski_372.f = @ (p) sumsq (p(4:9));
 685
 686   ret.general.schittkowski_391.dfdp = [];
 687   ret.general.schittkowski_391.init_p = ...
 688       -2.8742711 * alpha_391 (zeros (30, 1), 1:30);
 689   %% computed with fminunc (Octave)
 690   ret.general.schittkowski_391.result.p = [-1.1986682e+18; ...
 691                                            -1.1474574e+07; ...
 692                                            -1.3715802e+07; ...
 693                                            -1.0772255e+07; ...
 694                                            -1.0634232e+07; ...
 695                                            -1.0622915e+07; ...
 696                                            -8.8775399e+06; ...
 697                                            -8.8201496e+06; ...
 698                                            -9.7729975e+06; ...
 699                                            -1.0431808e+07; ...
 700                                            -1.0415089e+07; ...
 701                                            -1.0350400e+07; ...
 702                                            -1.0325094e+07; ...
 703                                            -1.0278561e+07; ...
 704                                            -1.0275751e+07; ...
 705                                            -1.0276546e+07; ...
 706                                            -1.0292584e+07; ...
 707                                            -1.0289350e+07; ...
 708                                            -1.0192566e+07; ...
 709                                            -1.0058577e+07; ...
 710                                            -1.0096341e+07; ...
 711                                            -1.0242386e+07; ...
 712                                            -1.0615831e+07; ...
 713                                            -1.1142096e+07; ...
 714                                            -1.1617283e+07; ...
 715                                            -1.2005738e+07; ...
 716                                            -1.2282117e+07; ...
 717                                            -1.2301260e+07; ...
 718                                            -1.2051365e+07; ...
 719                                            -1.1704693e+07];
 720   ret.general.schittkowski_391.result.obj = -5.1615468e+20;
 721   ret.general.schittkowski_391.strict_inequc.bounds = [];
 722   ret.general.schittkowski_391.strict_inequc.linear = [];
 723   ret.general.schittkowski_391.strict_inequc.general = [];
 724   ret.general.schittkowski_391.non_strict_inequc.bounds = [];
 725   ret.general.schittkowski_391.non_strict_inequc.linear = [];
 726   ret.general.schittkowski_391.non_strict_inequc.general = [];
 727   ret.general.schittkowski_391.equc.linear = [];
 728   ret.general.schittkowski_391.equc.general = [];
 729   ret.general.schittkowski_391.f = @ (p) sum (alpha_391 (p, 1:30));
 730
 731   function ret = f_1 (x, p)
 732
 733     ret = p(1) + x(:, 1) ./ (p(2) * x(:, 2) + p(3) * x(:, 3));
 734
 735   function ret = f_2 (x, p, y)
 736
 737     y(3, 2) = p(4);
 738     x(9, 3) = p(5);
 739     p = p(:);
 740     mp = cat (2, p([1, 2, 3]), p([3, 1, 2]), p([3, 2, 1]));
 741     ret = x * mp - y;
 742
 743   function ret = f_3 (x, p)
 744
 745     ret = fixed_step_rk4 (x.', [1, 1, 0, 0, 0], 1, ...
 746                           @ (x, t) f_3_xdot (x, t, p));
 747     ret = ret.';
 748
 749   function ret = f_3_noweights (x, p)
 750
 751     ret = fixed_step_rk4 (x.', [.1, .1, 0, 0, 0], .2, ...
 752                           @ (x, t) f_3_xdot_noweights (x, t, p));
 753     ret = ret.';
 754
 755   function ret = f_3_xdot (x, t, p)
 756
 757     ret = zeros (5, 1);
 758     tp = p(2) * x(3) - p(1) * x(1) * x(2);
 759     ret(1) = tp;
 760     ret(2) = tp - p(4) * x(2) * x(3) + p(5) * x(5) - p(6) * x(2) * x(4);
 761     ret(3) = - tp - p(3) * x(3) - p(4) * x(2) * x(3);
 762     ret(4) = p(3) * x(3) + p(5) * x(5) - p(6) * x(2) * x(4);
 763     ret(5) = p(4) * x(2) * x(3) - p(5) * x(5) + p(6) * x(2) * x(4);
 764
 765   function ret = f_3_xdot_noweights (x, t, p)
 766
 767     x(1:2) = x(1:2) / .1;
 768     x(4) = x(4) / 10;
 769     x(5) = x(5) / 100;
 770     ret = f_3_xdot (x, t, p);
 771     ret(1:2) = ret(1:2) * .1;
 772     ret(4) = ret(4) * 10;
 773     ret(5) = ret(5) * 100;
 774
 775   function ret = f_r (x, p, hook)
 776
 777     n = size (hook.mc, 2);
 778     ns = cat (1, 0, cumsum (hook.ns));
 779     xdhook.p = p;
 780     ret = zeros (1, ns(end));
 781     %% temporary variables
 782     dls = p(3) ^ 2;
 783     dmhp = p(5) * dls / p(4);
 784     mhp = dmhp / 2;
 785     %%
 786     for id = 1:n
 787       xdhook.c = hook.mc(:, id);
 788       l = xdhook.c(3);
 789       x0 = mhp - sqrt (max (0, mhp ^ 2 + dls + (p(6) - l) * dmhp));
 790       ids = ns(id) + 1;
 791       ide = ns(id + 1);
 792
 793       tp = odeset ();
 794       %% necessary in Matlab (7.1)
 795       tp.OutputSave = [];
 796       tp.Refine = 0;
 797       %%
 798       tp.RelTol = 1e-7;
 799       tp.AbsTol = 1e-7;
 800       [cx, Xcx] = essential_ode23 (@ (t, X) f_r_xdot (X, t, xdhook), ...
 801                                    x([ids, ide]).', x0, tp);
 802       X = lin_interp (cx.', Xcx.', x(ids:ide).');
 803
 804       X = X.';
 805       [discarded, lr] = ...
 806           f_r_xdot (X, hook.t(ids:ide), xdhook);
 807       ret(ids:ide) = max (0, lr - p(6) - X) * p(5);
 808     end
 809     ret = ret.';
 810
 811   function [ret, l] = f_r_xdot (x, t, hook)
 812
 813     %% keep this working with non-scalar x and t
 814
 815     p = hook.p;
 816     c = hook.c;
 817     idl = t <= c(1);
 818     idg = t >= c(2);
 819     idb = ~ (idl | idg);
 820     l = zeros (size (t));
 821     l(idl) = c(3);
 822     l(idg) = c(4) * t(idg) + c(5);
 823     l(idb) = polyval (c(6:9), t(idb));
 824     dls = max (1e-6, l - p(6) - x);
 825     tf = x / p(3);
 826     ido = tf >= 1;
 827     idx = ~ido;
 828     ret(ido) = 0;
 829     ret(idx) = - ((p(4) + p(1)) * p(2)) ./ ...
 830         ((p(5) * dls(idx)) ./ (1 - tf(idx) .^ 2) + p(1)) + p(2);
 831
 832   function ret = alpha_391 (p, id)
 833
 834     %% for .general.schittkowski_391; id is a numeric index(-vector)
 835     %% into p
 836
 837     p = p(:);
 838     n = size (p, 1);
 839
 840     nid = length (id);
 841     id = reshape (id, 1, nid);
 842
 843     v = sqrt (repmat (p .^ 2, 1, nid) + 1 ./ ((1:n).') * id);
 844
 845     log_v = log (v);
 846
 847     ret = 420 * p(id) + (id(:) - 15) .^ 3 + ...
 848         sum (v .* (sin (log_v) .^ 5 + cos (log_v) .^ 5)).';
 849
 850   function ret = schittkowski_281_dfdp (p)
 851
 852     tp = (sum (((1:10).') .^ 3 .* (p - 1) .^ 2)) ^ (- 2 / 3) / 3;
 853
 854     ret = 2 * ((1:10).') .^ 3 .* (p - 1) * tp;
 855
 856   function state = fixed_step_rk4 (t, x0, step, f)
 857
 858     %% minimalistic fourth order ODE-solver, as said to be a popular one
 859     %% by Wikipedia (to make these optimization tests self-contained;
 860     %% for the same reason 'lookup' and even 'interp1' are not used
 861     %% here)
 862
 863     n = ceil ((t(end) - t(1)) / step) + 1;
 864     m = length (x0);
 865     tstate = zeros (m, n);
 866     tstate(:, 1) = x0;
 867     tt = linspace (t(1), t(1) + step * (n - 1), n);
 868     for id = 1 : n - 1
 869       k1 = f (tstate(:, id), tt(id));
 870       k2 = f (tstate(:, id) + .5 * step * k1, tt(id) + .5 * step);
 871       k3 = f (tstate(:, id) + .5 * step * k2, tt(id) + .5 * step);
 872       k4 = f (tstate(:, id) + step * k3, tt(id + 1));
 873       tstate(:, id + 1) = tstate(:, id) + ...
 874           (step / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
 875     end
 876     state = lin_interp (tt, tstate, t);
 877
 878   function ret = pc2 (p, id)
 879     %% a combination out of 2 possible values for each parameter
 880     r = size (p, 1);
 881     n = 2 ^ r;
 882     if (id < 0 || id > n)
 883       error ('no parameter set for this index');
 884     end
 885     if (id == 0) % return maximum id
 886       ret = n;
 887       return;
 888     end
 889     idx = dec2bin (id - 1, r) == '1';
 890     nidx = ~idx;
 891     ret = zeros (r, 1);
 892     ret(nidx) = p(nidx, 1);
 893     ret(idx) = p(idx, 2);
 894
 895   function [varargout] = essential_ode23 (vfun, vslot, vinit, vodeoptions)
 896
 897     %% This code is taken from the ode23 solver of Thomas Treichl
 898     %% <thomas.treichl@gmx.net>, some flexibility of the
 899     %% interface has been removed. The idea behind this duplication is
 900     %% to have a fixed version of the solver here which runs both in
 901     %% Octave and Matlab.
 902
 903     %% Some of the option treatment has been left out.
 904     if (length (vslot) > 2)
 905       vstepsizefixed = true;
 906     else
 907       vstepsizefixed = false;
 908     end
 909     if (strcmp (vodeoptions.NormControl, 'on'))
 910       vnormcontrol = true;
 911     else
 912       vnormcontrol = false;
 913     end
 914     if (~isempty (vodeoptions.NonNegative))
 915       if (isempty (vodeoptions.Mass))
 916         vhavenonnegative = true;
 917       else
 918         vhavenonnegative = false;
 919       end
 920     else
 921       vhavenonnegative = false;
 922     end
 923     if (isempty (vodeoptions.OutputFcn) && nargout == 0)
 924       vodeoptions.OutputFcn = @odeplot;
 925       vhaveoutputfunction = true;
 926     elseif (isempty (vodeoptions.OutputFcn))
 927       vhaveoutputfunction = false;
 928     else
 929       vhaveoutputfunction = true;
 930     end
 931     if (~isempty (vodeoptions.OutputSel))
 932       vhaveoutputselection = true;
 933     else
 934       vhaveoutputselection = false;
 935     end
 936     if (isempty (vodeoptions.OutputSave))
 937       vodeoptions.OutputSave = 1;
 938     end
 939     if (vodeoptions.Refine > 0)
 940       vhaverefine = true;
 941     else
 942       vhaverefine = false;
 943     end
 944     if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
 945       vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
 946       vodeoptions.InitialStep = vodeoptions.InitialStep / ...
 947           10^vodeoptions.Refine;
 948     end
 949     if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
 950       vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10;
 951     end
 952     if (~isempty (vodeoptions.Events))
 953       vhaveeventfunction = true;
 954     else
 955       vhaveeventfunction = false;
 956     end
 957     if (~isempty (vodeoptions.Mass) && ismatrix (vodeoptions.Mass))
 958       vhavemasshandle = false;
 959       vmass = vodeoptions.Mass;
 960     elseif (isa (vodeoptions.Mass, 'function_handle'))
 961       vhavemasshandle = true;
 962     else
 963       vhavemasshandle = false;
 964     end
 965     if (strcmp (vodeoptions.MStateDependence, 'none'))
 966       vmassdependence = false;
 967     else
 968       vmassdependence = true;
 969     end
 970
 971     %% Starting the initialisation of the core solver ode23
 972     vtimestamp  = vslot(1,1);           %% timestamp = start time
 973     vtimelength = length (vslot);       %% length needed if fixed steps
 974     vtimestop   = vslot(1,vtimelength); %% stop time = last value
 975     vdirection  = sign (vtimestop);     %% Flag for direction to solve
 976
 977     if (~vstepsizefixed)
 978       vstepsize = vodeoptions.InitialStep;
 979       vminstepsize = (vtimestop - vtimestamp) / (1/eps);
 980     else %% If step size is given then use the fixed time steps
 981       vstepsize = vslot(1,2) - vslot(1,1);
 982       vminstepsize = sign (vstepsize) * eps;
 983     end
 984
 985     vretvaltime = vtimestamp; %% first timestamp output
 986     vretvalresult = vinit;    %% first solution output
 987
 988     %% Initialize the OutputFcn
 989     if (vhaveoutputfunction)
 990       if (vhaveoutputselection) vretout = ...
 991             vretvalresult(vodeoptions.OutputSel);
 992       else
 993         vretout = vretvalresult;
 994       end
 995       feval (vodeoptions.OutputFcn, vslot.', ...
 996              vretout.', 'init');
 997     end
 998
 999     %% Initialize the EventFcn
1000     if (vhaveeventfunction)
1001       odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
1002                            vretvalresult.', 'init');
1003     end
1004
1005     vpow = 1/3;            %% 20071016, reported by Luis Randez
1006     va = [  0, 0, 0;       %% The Runge-Kutta-Fehlberg 2(3) coefficients
1007           1/2, 0, 0;       %% Coefficients proved on 20060827
1008           -1, 2, 0];      %% See p.91 in Ascher & Petzold
1009     vb2 = [0; 1; 0];       %% 2nd and 3rd order
1010     vb3 = [1/6; 2/3; 1/6]; %% b-coefficients
1011     vc = sum (va, 2);
1012
1013     %% The solver main loop - stop if the endpoint has been reached
1014     vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1,3);
1015     vcntiter = 0; vunhandledtermination = true; vcntsave = 2;
1016     while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
1017            (vdirection * (vstepsize) >= vdirection * (vminstepsize)))
1018
1019       %% Hit the endpoint of the time slot exactely
1020       if ((vtimestamp + vstepsize) > vdirection * vtimestop)
1021         %% if (((vtimestamp + vstepsize) > vtimestop) || ...
1022         %% (abs(vtimestamp + vstepsize - vtimestop) < eps))
1023         vstepsize = vtimestop - vdirection * vtimestamp;
1024       end
1025
1026       %% Estimate the three results when using this solver
1027       for j = 1:3
1028         vthetime  = vtimestamp + vc(j,1) * vstepsize;
1029         vtheinput = vu.' + vstepsize * vk(:,1:j-1) * va(j,1:j-1).';
1030         if (vhavemasshandle)   %% Handle only the dynamic mass matrix,
1031           if (vmassdependence) %% constant mass matrices have already
1032             vmass = feval ...  %% been set before (if any)
1033                 (vodeoptions.Mass, vthetime, vtheinput);
1034           else                 %% if (vmassdependence == false)
1035             vmass = feval ...  %% then we only have the time argument
1036                 (vodeoptions.Mass, vthetime);
1037           end
1038           vk(:,j) = vmass \ feval ...
1039               (vfun, vthetime, vtheinput);
1040         else
1041           vk(:,j) = feval ...
1042               (vfun, vthetime, vtheinput);
1043         end
1044       end
1045
1046       %% Compute the 2nd and the 3rd order estimation
1047       y2 = vu.' + vstepsize * (vk * vb2);
1048       y3 = vu.' + vstepsize * (vk * vb3);
1049       if (vhavenonnegative)
1050         vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
1051         y2(vodeoptions.NonNegative) = abs (y2(vodeoptions.NonNegative));
1052         y3(vodeoptions.NonNegative) = abs (y3(vodeoptions.NonNegative));
1053       end
1054       vSaveVUForRefine = vu;
1055
1056       %% Calculate the absolute local truncation error and the
1057       %% acceptable error
1058       if (~vstepsizefixed)
1059         if (~vnormcontrol)
1060           vdelta = abs (y3 - y2);
1061           vtau = max (vodeoptions.RelTol * abs (vu.'), ...
1062                       vodeoptions.AbsTol);
1063         else
1064           vdelta = norm (y3 - y2, Inf);
1065           vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), ...
1066                                                 1.0), ...
1067                       vodeoptions.AbsTol);
1068         end
1069       else %% if (vstepsizefixed == true)
1070         vdelta = 1; vtau = 2;
1071       end
1072
1073       %% If the error is acceptable then update the vretval variables
1074       if (all (vdelta <= vtau))
1075         vtimestamp = vtimestamp + vstepsize;
1076         vu = y3.'; %% MC2001: the higher order estimation as 'local
1077         %% extrapolation' Save the solution every vodeoptions.OutputSave
1078         %% steps
1079         if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
1080           vretvaltime(vcntsave,:) = vtimestamp;
1081           vretvalresult(vcntsave,:) = vu;
1082           vcntsave = vcntsave + 1;
1083         end
1084         vcntloop = vcntloop + 1; vcntiter = 0;
1085
1086         %% Call plot only if a valid result has been found, therefore
1087         %% this code fragment has moved here. Stop integration if plot
1088         %% function returns false
1089         if (vhaveoutputfunction)
1090           for vcnt = 0:vodeoptions.Refine %% Approximation between told
1091             %% and t
1092             if (vhaverefine)              %% Do interpolation
1093               vapproxtime = (vcnt + 1) * vstepsize / ...
1094                   (vodeoptions.Refine + 2);
1095               vapproxvals = vSaveVUForRefine.' + vapproxtime * (vk * ...
1096                                                                 vb3);
1097               vapproxtime = (vtimestamp - vstepsize) + vapproxtime;
1098             else
1099               vapproxvals = vu.';
1100               vapproxtime = vtimestamp;
1101             end
1102             if (vhaveoutputselection)
1103               vapproxvals = vapproxvals(vodeoptions.OutputSel);
1104             end
1105             vpltret = feval (vodeoptions.OutputFcn, vapproxtime, ...
1106                              vapproxvals, []);
1107             if vpltret %% Leave refinement loop
1108               break;
1109             end
1110           end
1111           if (vpltret) %% Leave main loop
1112             vunhandledtermination = false;
1113             break;
1114           end
1115         end
1116
1117         %% Call event only if a valid result has been found, therefore
1118         %% this code fragment has moved here. Stop integration if
1119         %% veventbreak is true
1120         if (vhaveeventfunction)
1121           vevent = ...
1122               odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
1123                                    vu(:), []);
1124           if (~isempty (vevent{1}) && vevent{1} == 1)
1125             vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
1126             vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
1127             vunhandledtermination = false; break;
1128           end
1129         end
1130       end %% If the error is acceptable ...
1131
1132       %% Update the step size for the next integration step
1133       if (~vstepsizefixed)
1134         %% 20080425, reported by Marco Caliari vdelta cannot be negative
1135         %% (because of the absolute value that has been introduced) but
1136         %% it could be 0, then replace the zeros with the maximum value
1137         %% of vdelta
1138         vdelta(find (vdelta == 0)) = max (vdelta);
1139         %% It could happen that max (vdelta) == 0 (ie. that the original
1140         %% vdelta was 0), in that case we double the previous vstepsize
1141         vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
1142
1143         if (vdirection == 1)
1144           vstepsize = min (vodeoptions.MaxStep, ...
1145                            min (0.8 * vstepsize * (vtau ./ vdelta) .^ ...
1146                                 vpow));
1147         else
1148           vstepsize = max (vodeoptions.MaxStep, ...
1149                            max (0.8 * vstepsize * (vtau ./ vdelta) .^ ...
1150                                 vpow));
1151         end
1152       else %% if (vstepsizefixed)
1153         if (vcntloop <= vtimelength)
1154           vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
1155         else %% Get out of the main integration loop
1156           break;
1157         end
1158       end
1159
1160       %% Update counters that count the number of iteration cycles
1161       vcntcycles = vcntcycles + 1; %% Needed for cost statistics
1162       vcntiter = vcntiter + 1;     %% Needed to find iteration problems
1163
1164       %% Stop solving because the last 1000 steps no successful valid
1165       %% value has been found
1166       if (vcntiter >= 5000)
1167         error (['Solving has not been successful. The iterative', ...
1168         ' integration loop exited at time t = %f before endpoint at', ...
1169         ' tend = %f was reached. This happened because the iterative', ...
1170         ' integration loop does not find a valid solution at this time', ...
1171         ' stamp. Try to reduce the value of ''InitialStep'' and/or', ...
1172         ' ''MaxStep'' with the command ''odeset''.\n'], vtimestamp, vtimestop);
1173       end
1174
1175     end %% The main loop
1176
1177     %% Check if integration of the ode has been successful
1178     if (vdirection * vtimestamp < vdirection * vtimestop)
1179       if (vunhandledtermination == true)
1180         error ('OdePkg:InvalidArgument', ...
1181        ['Solving has not been successful. The iterative', ...
1182         ' integration loop exited at time t = %f', ...
1183         ' before endpoint at tend = %f was reached. This may', ...
1184         ' happen if the stepsize grows smaller than defined in', ...
1185         ' vminstepsize. Try to reduce the value of ''InitialStep'' and/or', ...
1186         ' ''MaxStep'' with the command ''odeset''.\n'], vtimestamp, vtimestop);
1187       else
1188         warning ('OdePkg:InvalidArgument', ...
1189          ['Solver has been stopped by a call of ''break'' in', ...
1190           ' the main iteration loop at time t = %f before endpoint at', ...
1191           ' tend = %f was reached. This may happen because the @odeplot', ...
1192           ' function returned ''true'' or the @event function returned ''true''.'], ...
1193          vtimestamp, vtimestop);
1194       end
1195     end
1196
1197     %% Postprocessing, do whatever when terminating integration
1198     %% algorithm
1199     if (vhaveoutputfunction) %% Cleanup plotter
1200       feval (vodeoptions.OutputFcn, vtimestamp, ...
1201              vu.', 'done');
1202     end
1203     if (vhaveeventfunction)  %% Cleanup event function handling
1204       odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
1205                            vu.', 'done');
1206     end
1207     %% Save the last step, if not already saved
1208     if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
1209       vretvaltime(vcntsave,:) = vtimestamp;
1210       vretvalresult(vcntsave,:) = vu;
1211     end
1212
1213
1214     varargout{1} = vretvaltime;     %% Time stamps are first output argument
1215     varargout{2} = vretvalresult;   %% Results are second output argument
1216
1217   function yi = lin_interp (x, y, xi)
1218
1219     %% Actually interp1 with 'linear' should behave equally in Octave
1220     %% and Matlab, but having this subset of functionality here is being
1221     %% on the safe side.
1222
1223     n = size (x, 2);
1224     m = size (y, 1);
1225     %% This elegant lookup is from an older version of 'lookup' by Paul
1226     %% Kienzle, and had been suggested by Kai Habel <kai.habel@gmx.de>.
1227     [v, p] = sort ([x, xi]);
1228     idx(p) = cumsum (p <= n);
1229     idx = idx(n + 1 : n + size (xi, 2));
1230     %%
1231     idx(idx == n) = n - 1;
1232     yi = y(:, idx) + ...
1233         repmat (xi - x(idx), m, 1) .* ...
1234         (y(:, idx + 1) - y(:, idx)) ./ ...
1235         repmat (x(idx + 1) - x(idx), m, 1);
1236
1237   function ret = apply_idx_if_given  (ret, idx)
1238
1239     if (nargin > 1)
1240       ret = ret(idx);
1241     end
1242
1243   function fval = scalar_ifelse (cond, tval, fval)
1244
1245     %% needed for some anonymous functions, builtin ifelse only available
1246     %% in Octave > 3.2; we need only the scalar case here
1247
1248     if (cond)
1249       fval = tval;
1250     end
1251
1252 %!demo
1253 %! p_t = optim_problems ().curve.p_1;
1254 %! global verbose;
1255 %! verbose = false;
1256 %! [cy, cp, cvg, iter] = leasqr (p_t.data.x, p_t.data.y, p_t.init_p, p_t.f)
1257 %! disp (p_t.result.p)
1258 %! sumsq (cy - p_t.data.y)