Mercurial > mxe-octave
view src/of-optim-2-deprecated.patch @ 4942:805494784837
of-optim: Cherry-picked upstream fix (bug #55325).
* src/of-optim-2-deprecated.patch: new file
* dist-files.mk: add ref to of-optim-2-deprecated.patch
author | Markus Mützel <markus.muetzel@gmx.de> |
---|---|
date | Sun, 20 Jan 2019 17:59:50 +0100 |
parents | |
children |
line wrap: on
line source
# HG changeset patch # User Olaf Till <i7tiol@t-online.de> # Date 1547218871 -3600 # Fri Jan 11 16:01:11 2019 +0100 # Node ID 5fcbf411d5c2a4b076e3ababeb6126125243e49d # Parent 158b61118c7557f549656bf76ddc144c70bd3aa5 remove deprecated samin.cc, still available as backend to nonlin_min diff -r 158b61118c75 -r 5fcbf411d5c2 src/Makefile.in --- a/src/Makefile.in Fri Jan 11 15:29:45 2019 +0100 +++ b/src/Makefile.in Fri Jan 11 16:01:11 2019 +0100 @@ -46,7 +46,7 @@ TEXIFILE := $(addsuffix .texi,$(basename $(INFOFILE))) TXIFILE := $(addsuffix .txi,$(basename $(INFOFILE))) HTMLDIR := ../doc/html/ -OCTFILES := __bfgsmin.oct numgradient.oct numhessian.oct samin.oct __disna_optim__.oct +OCTFILES := __bfgsmin.oct numgradient.oct numhessian.oct __disna_optim__.oct OCTSOURCEFILES := $(addsuffix .cc,$(basename $(OCTFILES))) DSFILES := $(addsuffix .docstrings,$(OCTSOURCEFILES)) diff -r 158b61118c75 -r 5fcbf411d5c2 src/samin.cc --- a/src/samin.cc Fri Jan 11 15:29:45 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,549 +0,0 @@ -// Copyright (C) 2004, 2006 Michael Creel <michael.creel@uab.es> -// -// This program is free software; you can redistribute it and/or modify it under -// the terms of the GNU General Public License as published by the Free Software -// Foundation; either version 3 of the License, or (at your option) any later -// version. -// -// This program is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -// details. -// -// You should have received a copy of the GNU General Public License along with -// this program; if not, see <http://www.gnu.org/licenses/>. - -// References: -// -// The code follows the article -// Goffe, William L. (1996) "SIMANN: A Global Optimization Algorithm -// using Simulated Annealing " Studies in Nonlinear Dynamics & Econometrics -// Oct96, Vol. 1 Issue 3. -// -// The code uses the same names for control variables, -// for the most part. A notable difference is that the initial -// temperature is found automatically to ensure that the active -// bounds when the temperature begins to reduce cover the entire -// parameter space (defined as a n-dimensional -// rectangle that is the Cartesian product of the -// (lb_i, ub_i), i = 1,2,..n -// -// Also of note: -// Corana et. al., (1987) "Minimizing Multimodal Functions of Continuous -// Variables with the "Simulated Annealing" Algorithm", -// ACM Transactions on Mathematical Software, V. 13, N. 3. -// -// Goffe, et. al. (1994) "Global Optimization of Statistical Functions -// with Simulated Annealing", Journal of Econometrics, -// V. 60, N. 1/2. - -#include <oct.h> -#include <octave/parse.h> -#include <octave/Cell.h> -#include <octave/lo-mappers.h> -#include <octave/oct-rand.h> -#include <float.h> - -#include "error-helpers.h" - -// define argument checks -static bool any_bad_argument(const octave_value_list& args) -{ - - // objective function name is a string? - if (!args(0).is_string()) - { - _p_error("samin: first argument must be string holding objective function name"); - return true; - } - - // are function arguments contained in a cell? - if (!args(1).OV_ISCELL ()) - { - _p_error("samin: second argument must cell array of function arguments"); - return true; - } - - // is control a cell? - Cell control; - bool err; - SET_ERR (control =args(2).cell_value(), err); - if (err) - { - _p_error("samin: third argument must cell array of algorithm controls"); - return true; - } - - // does control have proper number of elements? - if (!(control.numel () == 11)) - { - _p_error("samin: third argument must be a cell array with 11 elements"); - return true; - } - - // now check type of each element of control - if (!(control(0).is_real_matrix()) && !(control(0).is_real_scalar())) - { - _p_error("samin: 1st element of controls must be LB: a vector of lower bounds"); - return true; - } - - if ((control(0).is_real_matrix()) && (control(0).columns() != 1)) - { - _p_error("samin: 1st element of controls must be LB: a vector of lower bounds"); - return true; - } - - if (!(control(1).is_real_matrix()) && !(control(1).is_real_scalar())) - { - _p_error("samin: 1st element of controls must be UB: a vector of lower bounds"); - return true; - } - - if ((control(1).is_real_matrix()) && (control(1).columns() != 1)) - { - _p_error("samin: 2nd element of controls must be UB: a vector of lower bounds"); - return true; - } - - int tmp; - SET_ERR (tmp = control(2).int_value(), err); - if (err || tmp < 1) - { - _p_error("samin: 3rd element of controls must be NT: positive integer\n\ -loops per temperature reduction"); - return true; - } - - SET_ERR (tmp = control(3).int_value(), err); - if (err || tmp < 1) - { - _p_error("samin: 4th element of controls must be NS: positive integer\n\ -loops per stepsize adjustment"); - return true; - } - - double tmp2; - SET_ERR (tmp2 = control(4).double_value(), err); - if (err || tmp < 0) - { - _p_error("samin: 5th element of controls must be RT:\n\ -temperature reduction factor, RT > 0"); - return true; - } - - SET_ERR (tmp2 = control(5).double_value(), err); - if (err || tmp < 0) - { - _p_error("samin: 6th element of controls must be integer MAXEVALS > 0 "); - return true; - } - - SET_ERR (tmp = control(6).int_value(), err); - if (err || tmp < 0) - { - _p_error("samin: 7th element of controls must be NEPS: positive integer\n\ -number of final obj. values that must be within EPS of eachother "); - return true; - } - - SET_ERR (tmp2 = control(7).double_value(), err); - if (err || tmp2 < 0) - { - _p_error("samin: 8th element of controls must must be FUNCTOL (> 0)\n\ -used to compare the last NEPS obj values for convergence test"); - return true; - } - - SET_ERR (tmp2 = control(8).double_value(), err); - if (err || tmp2 < 0) - { - _p_error("samin: 9th element of controls must must be PARAMTOL (> 0)\n\ -used to compare the last NEPS obj values for convergence test"); - return true; - } - - SET_ERR (tmp = control(9).int_value(), err); - if (err || tmp < 0 || tmp > 2) - { - _p_error("samin: 9th element of controls must be VERBOSITY (0, 1, or 2)"); - return true; - } - - SET_ERR (tmp = control(10).int_value(), err); - if (err || tmp < 0) - { - _p_error("samin: 10th element of controls must be MINARG (integer)\n\ - position of argument to minimize wrt"); - return true; - } - - // make sure that minarg points to an existing element - if ((tmp > args(1).length ()) || (tmp < 1)) - { - _p_error("bfgsmin: 4th argument must be a positive integer that indicates \n\ -which of the elements of the second argument is the one minimization is over"); - return true; - } - - return false; -} - -//-------------- The annealing algorithm -------------- -DEFUN_DLD(samin, args, , "samin: simulated annealing minimization of a function. See samin_example.m\n\ -\n\ -samin will be removed from a future version of the optim package.\n\ -Equivalent functionality is now in the samin backend of nonlin_min.\n\ -\n\ -usage: [x, obj, convergence, details] = samin(\"f\", {args}, {control})\n\ -\n\ -Arguments:\n\ -* \"f\": function name (string)\n\ -* {args}: a cell array that holds all arguments of the function,\n\ -* {control}: a cell array with 11 elements\n\ - * LB - vector of lower bounds\n\ - * UB - vector of upper bounds\n\ - * nt - integer: # of iterations between temperature reductions\n\ - * ns - integer: # of iterations between bounds adjustments\n\ - * rt - (0 < rt <1): temperature reduction factor\n\ - * maxevals - integer: limit on function evaluations\n\ - * neps - integer: number of values final result is compared to\n\ - * functol - (> 0): the required tolerance level for function value\n\ - comparisons\n\ - * paramtol - (> 0): the required tolerance level for parameters\n\ - * verbosity - scalar: 0, 1, or 2.\n\ - * 0 = no screen output\n\ - * 1 = only final results to screen\n\ - * 2 = summary every temperature change\n\ - * minarg - integer: which of function args is minimization over?\n\ -\n\ -Returns:\n\ -* x: the minimizer\n\ -* obj: the value of f() at x\n\ -* convergence:\n\ - 0 if no convergence within maxevals function evaluations\n\ - 1 if normal convergence to a point interior to the parameter space\n\ - 2 if convergence to point very near bounds of parameter space\n\ - (suggest re-running with looser bounds)\n\ -* details: a px3 matrix. p is the number of times improvements were found.\n\ - The columns record information at the time an improvement was found\n\ - * first: cumulative number of function evaluations\n\ - * second: temperature\n\ - * third: function value\n\ -\n\ -Example: see samin_example\n\ -") -{ - static bool warned = false; - if (! warned) - { - warned = true; - warning_with_id ("Octave:deprecated-function", - "samin will be removed from a future version of the optim package, equivalent functionality is now in the samin backend of nonlin_min"); - } - - int nargin = args.length(); - if (!(nargin == 3)) - { - error("samin: you must supply 3 arguments"); - return octave_value_list(); - } - - // check the arguments - if (any_bad_argument (args)) - { - error ("error in samin"); - return octave_value_list(); - } - - std::string obj_fn (args(0).string_value()); - Cell f_args_cell = args(1).cell_value (); // args to obj fn come in as a cell to allow arbitrary number - Cell control (args(2).cell_value()); - - octave_value_list f_args; - octave_value_list f_return; // holder for feval returns - - int m, i, j, k, h, n, nacc, func_evals; - int nup, nrej, nnew, ndown, lnobds; - int converge, test, coverage_ok; - - // user provided controls - const ColumnVector lb (control(0).column_vector_value()); - const ColumnVector ub (control(1).column_vector_value()); - const int nt (control(2).int_value()); - const int ns (control(3).int_value()); - const double rt (control(4).double_value()); - const int maxevals (control(5).int_value()); - const int neps (control(6).int_value()); - const double functol (control(7).double_value()); - const double paramtol (control(8).double_value()); - const int verbosity (control(9).int_value()); - const int minarg (control(10).int_value()); - - // type checking for minimization parameter done here, since we don't know minarg - // until now - if (!(f_args_cell(minarg - 1).is_real_matrix() || (f_args_cell(minarg - 1).is_real_scalar()))) - { - error("samin: minimization must be with respect to a column vector"); - return octave_value_list(); - } - if ((f_args_cell(minarg - 1).is_real_matrix()) && (f_args_cell(minarg - 1).columns() != 1)) - { - error("samin: minimization must be with respect to a column vector"); - return octave_value_list(); - } - - double f, fp, p, fopt, rand_draw, ratio, t; - - Matrix details(1,3); // record function evaluations, temperatures and function values - RowVector info(3); - - // copy cell contents over to octave_value_list to use feval() - k = f_args_cell.numel (); - f_args(k); // resize only once - for (i = 0; i<k; i++) f_args(i) = f_args_cell(i); - - ColumnVector x = f_args(minarg - 1).column_vector_value(); - ColumnVector bounds = ub - lb; - n = x.rows(); - ColumnVector xopt = x; - ColumnVector xp(n); - ColumnVector nacp(n); - - // Set initial values - nacc = 0; // total accepted trials - t = 1000.0; // temperature - will initially rise or fall to cover parameter space. Then it will fall - converge = 0; // convergence indicator 0 (failure), 1 (normal success), or 2 (convergence but near bounds) - coverage_ok = 0; // has parameter space been covered? When turns to 1, temperature starts to fall - // most recent values, to compare to when checking convergend - ColumnVector fstar(neps,1); - fstar.fill(DBL_MAX); - octave_rand::distribution("uniform"); // we'll be using draws from U(0,1) - - // check for out-of-bounds starting values - for(i = 0; i < n; i++) - { - if(( x(i) > ub(i)) || (x(i) < lb(i))) - { - error("samin: initial parameter %d out of bounds", i+1); - return octave_value_list(); - } - } - - // Initial obj_value - f_return = OCTAVE__FEVAL (obj_fn, f_args); - f = f_return(0).double_value(); - fopt = f; // give it something to compare to - func_evals = 0; // total function evaluations (limited by maxeval) - details(0,0) = func_evals; - details(0,1) = t; - details(0,2) = fopt; - - // main loop, first increase temperature until parameter space covered, then reduce until convergence - while(converge==0) - { - // statistics to report at each temp change, set back to zero - nup = 0; - nrej = 0; - nnew = 0; - ndown = 0; - lnobds = 0; - - // repeat nt times then adjust temperature - for(m = 0;m < nt;m++) - { - // repeat ns times, then adjust bounds - for(j = 0;j < ns;j++) - { - // generate new point by taking last and adding a random value - // to each of elements, in turn - for(h = 0;h < n;h++) - { - // new Sept 2011, if bounds are same, skip the search for that vbl. Allows restrictions without complicated programming - if (lb(h) != ub(h)) - { - xp = x; - rand_draw = octave_rand::scalar(); - xp(h) = x(h) + (2.0 * rand_draw - 1.0) * bounds(h); - if ((xp(h) < lb(h)) || (xp(h) > ub(h))) - { - rand_draw = octave_rand::scalar(); // change 07-Nov-2007: avoid correlation with hitting bounds - xp(h) = lb(h) + (ub(h) - lb(h)) * rand_draw; - lnobds = lnobds + 1; - } - // Evaluate function at new point - f_args(minarg - 1) = xp; - f_return = OCTAVE__FEVAL (obj_fn, f_args); - fp = f_return(0).double_value(); - func_evals = func_evals + 1; - // Accept the new point if the function value decreases - if(fp <= f) - { - x = xp; - f = fp; - nacc = nacc + 1; // total number of acceptances - nacp(h) = nacp(h) + 1; // acceptances for this parameter - nup = nup + 1; - // If lower than any other point, record as new optimum - if(fp < fopt) - { - xopt = xp; - fopt = fp; - nnew = nnew + 1; - info(0) = func_evals; - info(1) = t; - info(2) = fp; - details = details.stack(info); - } - } - // If the point is higher, use the Metropolis criteria to decide on - // acceptance or rejection. - else - { - p = exp(-(fp - f) / t); - rand_draw = octave_rand::scalar(); - if(rand_draw < p) - { - x = xp; - f = fp; - nacc = nacc + 1; - nacp(h) = nacp(h) + 1; - ndown = ndown + 1; - } - else nrej = nrej + 1; - } - } - // If maxevals exceeded, terminate the algorithm - if (func_evals >= maxevals) - { - if (verbosity >= 1) - { - printf("\n================================================\n"); - printf("SAMIN results\n"); - printf("NO CONVERGENCE: MAXEVALS exceeded\n"); - printf("================================================\n"); - printf("Convergence tolerances: Func. tol. %e Param. tol. %e\n", functol, paramtol); - printf("Obj. fn. value %f\n\n", fopt); - printf(" parameter search width\n"); - for(i = 0; i < n; i++) printf("%20f%20f\n", xopt(i), bounds(i)); - } - f_return(3) = details; - f_return(2) = 0; - f_return(1) = fopt; - f_return(0) = xopt; - return octave_value_list(f_return); - } - } - } - // Adjust bounds so that approximately half of all evaluations are accepted - test = 0; - for(i = 0;i < n;i++) - { - if (lb(i) != ub(i)) - { - ratio = nacp(i) / ns; - if(ratio > 0.6) bounds(i) = bounds(i) * (1.0 + 2.0 * (ratio - 0.6) / 0.4); - else if (ratio < .4) - bounds(i) = bounds(i) / (1.0 + 2.0 * ((0.4 - ratio) / 0.4)); - // keep within initial bounds - if(bounds(i) >= (ub(i) - lb(i))) - { - bounds(i) = ub(i) - lb(i); - test = test + 1; - } - } - else - test = test + 1; // make sure coverage check passes for the fixed parameters - } - nacp.fill(0.0); - // check if we cover parameter space. if we have yet to do so - if (!coverage_ok) coverage_ok = (test == n); - } - // intermediate output, if desired - if(verbosity == 2) - { - printf("\nsamin: intermediate results before next temperature change\n"); - printf("\ntemperature %e", t); - printf("\ncurrent best function value %f", fopt); - printf("\ntotal evaluations so far %d", func_evals); - printf("\ntotal moves since last temperature reduction %d", nup + ndown + nrej); - printf("\ndownhill %d", nup); - printf("\naccepted uphill %d", ndown); - printf("\nrejected uphill %d", nrej); - printf("\nout of bounds trials %d", lnobds); - printf("\nnew minima this temperature %d", nnew); - printf("\n\n parameter search width\n"); - for(i = 0; i < n; i++) printf("%20f%20f\n", xopt(i), bounds(i)); - printf("\n"); - } - // Check for convergence, if we have covered the parameter space - if (coverage_ok) - { - // last value close enough to last neps values? - fstar(0) = f; - test = 0; - for (i = 1; i < neps; i++) test = test + fabs(f - fstar(i)) > functol; - test = (test > 0); // if different from zero, function conv. has failed - // last value close enough to overall best? - if (((fopt - f) <= functol) && (!test)) - { - // check for bound narrow enough for parameter convergence - for (i = 0;i < n;i++) - { - if (bounds(i) > paramtol) - { - converge = 0; // no conv. if bounds too wide - break; - } - else - converge = 1; - } - } - // check if optimal point is near boundary of parameter space, and change convergence message if so - if (converge && lnobds > 0) converge = 2; - // Like to see the final results? - if (converge > 0) - { - if (verbosity >= 1) - { - printf("\n================================================\n"); - printf("SAMIN results\n\n"); - if (converge == 1) printf("==> Normal convergence <==\n\n"); - if (converge == 2) - { - printf("==> WARNING <==: Last point satisfies convergence criteria,\n"); - printf("but is near boundary of parameter space.\n"); - printf("%d out of %d evaluations were out-of-bounds in the last round.\n", lnobds, (nup+ndown+nrej)); - printf("Expand bounds and re-run, unless this is a constrained minimization.\n\n"); - } - printf("Convergence tolerances:\nFunction: %e\nParameters: %e\n", functol, paramtol); - printf("\nObjective function value at minimum: %f\n\n", fopt); - printf(" parameter search width\n"); - for(i = 0; i < n; i++) printf("%20f%20f\n", xopt(i), bounds(i)); - printf("================================================\n"); - } - f_return(3) = details; - f_return(2) = converge; - f_return(1) = fopt; - f_return(0) = xopt; - return f_return; // this breaks out, if we get here - } - // Reduce temperature, record current function value in the - // list of last "neps" values, and loop again - t = rt * t; - for(i = neps-1; i > 0; i--) fstar(i) = fstar(i-1); - f = fopt; - x = xopt; - } - else - { - // coverage not ok - increase temperature quickly to expand search area, to cover parameter space - t = t*t; - for(i = neps-1; i > 0; i--) fstar(i) = fstar(i-1); - f = fopt; - x = xopt; - } - } - // silence compiler warning - return octave_value_list (); -}