X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Foptim-1.2.0%2Fde_min.m;fp=octave_packages%2Foptim-1.2.0%2Fde_min.m;h=93b4c26e86991c45b2a5bbc8f79684ab9d1f915f;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/optim-1.2.0/de_min.m b/octave_packages/optim-1.2.0/de_min.m new file mode 100644 index 0000000..93b4c26 --- /dev/null +++ b/octave_packages/optim-1.2.0/de_min.m @@ -0,0 +1,451 @@ +## Copyright (C) 1996, 1997 R. Storn +## Copyright (C) 2009-2010 Christian Fischer +## +## 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 . + +## de_min: global optimisation using differential evolution +## +## Usage: [x, obj_value, nfeval, convergence] = de_min(fcn, control) +## +## minimization of a user-supplied function with respect to x(1:D), +## using the differential evolution (DE) method based on an algorithm +## by Rainer Storn (http://www.icsi.berkeley.edu/~storn/code.html) +## See: http://www.softcomputing.net/tevc2009_1.pdf +## +## +## Arguments: +## --------------- +## fcn string : Name of function. Must return a real value +## control vector : (Optional) Control variables, described below +## or struct +## +## Returned values: +## ---------------- +## x vector : parameter vector of best solution +## obj_value scalar : objective function value of best solution +## nfeval scalar : number of function evaluations +## convergence : 1 = best below value to reach (VTR) +## 0 = population has reached defined quality (tol) +## -1 = some values are close to constraints/boundaries +## -2 = max number of iterations reached (maxiter) +## -3 = max number of functions evaluations reached (maxnfe) +## +## Control variable: (optional) may be named arguments (i.e. "name",value +## ---------------- pairs), a struct, or a vector, where +## NaN's are ignored. +## +## XVmin : vector of lower bounds of initial population +## *** note: by default these are no constraints *** +## XVmax : vector of upper bounds of initial population +## constr : 1 -> enforce the bounds not just for the initial population +## const : data vector (remains fixed during the minimization) +## NP : number of population members +## F : difference factor from interval [0, 2] +## CR : crossover probability constant from interval [0, 1] +## strategy : 1 --> DE/best/1/exp 7 --> DE/best/1/bin +## 2 --> DE/rand/1/exp 8 --> DE/rand/1/bin +## 3 --> DE/target-to-best/1/exp 9 --> DE/target-to-best/1/bin +## 4 --> DE/best/2/exp 10--> DE/best/2/bin +## 5 --> DE/rand/2/exp 11--> DE/rand/2/bin +## 6 --> DEGL/SAW/exp else DEGL/SAW/bin +## refresh : intermediate output will be produced after "refresh" +## iterations. No intermediate output will be produced +## if refresh is < 1 +## VTR : Stopping criterion: "Value To Reach" +## de_min will stop when obj_value <= VTR. +## Use this if you know which value you expect. +## tol : Stopping criterion: "tolerance" +## stops if (best-worst)/max(1,worst) < tol +## This stops basically if the whole population is "good". +## maxnfe : maximum number of function evaluations +## maxiter : maximum number of iterations (generations) +## +## The algorithm seems to work well only if [XVmin,XVmax] covers the +## region where the global minimum is expected. +## DE is also somewhat sensitive to the choice of the +## difference factor F. A good initial guess is to choose F from +## interval [0.5, 1], e.g. 0.8. +## CR, the crossover probability constant from interval [0, 1] +## helps to maintain the diversity of the population and is +## rather uncritical but affects strongly the convergence speed. +## If the parameters are correlated, high values of CR work better. +## The reverse is true for no correlation. +## Experiments suggest that /bin likes to have a slightly +## larger CR than /exp. +## The number of population members NP is also not very critical. A +## good initial guess is 10*D. Depending on the difficulty of the +## problem NP can be lower than 10*D or must be higher than 10*D +## to achieve convergence. +## +## Default Values: +## --------------- +## XVmin = [-2]; +## XVmax = [ 2]; +## constr= 0; +## const = []; +## NP = 10 *D +## F = 0.8; +## CR = 0.9; +## strategy = 12; +## refresh = 0; +## VTR = -Inf; +## tol = 1.e-3; +## maxnfe = 1e6; +## maxiter = 1000; +## +## +## Example to find the minimum of the Rosenbrock saddle: +## ---------------------------------------------------- +## Define f as: +## function result = f(x); +## result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2; +## end +## Then type: +## +## ctl.XVmin = [-2 -2]; +## ctl.XVmax = [ 2 2]; +## [x, obj_value, nfeval, convergence] = de_min (@f, ctl); +## +## Keywords: global-optimisation optimisation minimisation + +function [bestmem, bestval, nfeval, convergence] = de_min(fcn, varargin) + +## Default options +XVmin = [-2 ]; +XVmax = [ 2 ]; +constr= 0; +const = []; +NP = 0; # NP will be set later +F = 0.8; +CR = 0.9; +strategy = 12; +refresh = 0; +VTR = -Inf; +tol = 1.e-3; +maxnfe = 1e6; +maxiter = 1000; + +## ------------ Check input variables (ctl) -------------------------------- +if nargin >= 2, # Read control arguments + va_arg_cnt = 1; + if nargin > 2, + ctl = struct (varargin{:}); + else + ctl = varargin{va_arg_cnt++}; + end + if isnumeric (ctl) + if length (ctl)>=1 && !isnan (ctl(1)), XVmin = ctl(1); end + if length (ctl)>=2 && !isnan (ctl(2)), XVmax = ctl(2); end + if length (ctl)>=3 && !isnan (ctl(3)), constr = ctl(3); end + if length (ctl)>=4 && !isnan (ctl(4)), const = ctl(4); end + if length (ctl)>=5 && !isnan (ctl(5)), NP = ctl(5); end + if length (ctl)>=6 && !isnan (ctl(6)), F = ctl(6); end + if length (ctl)>=7 && !isnan (ctl(7)), CR = ctl(7); end + if length (ctl)>=8 && !isnan (ctl(8)), strategy = ctl(8); end + if length (ctl)>=9 && !isnan (ctl(9)), refresh = ctl(9); end + if length (ctl)>=10&& !isnan (ctl(10)), VTR = ctl(10); end + if length (ctl)>=11&& !isnan (ctl(11)), tol = ctl(11); end + if length (ctl)>=12&& !isnan (ctl(12)), maxnfe = ctl(12); end + if length (ctl)>=13&& !isnan (ctl(13)), maxiter = ctl(13); end + else + if isfield (ctl,"XVmin") && !isnan (ctl.XVmin), XVmin=ctl.XVmin; end + if isfield (ctl,"XVmax") && !isnan (ctl.XVmax), XVmax=ctl.XVmax; end + if isfield (ctl,"constr")&& !isnan (ctl.constr), constr=ctl.constr; end + if isfield (ctl,"const") && !isnan (ctl.const), const=ctl.const; end + if isfield (ctl, "NP" ) && ! isnan (ctl.NP ), NP = ctl.NP ; end + if isfield (ctl, "F" ) && ! isnan (ctl.F ), F = ctl.F ; end + if isfield (ctl, "CR" ) && ! isnan (ctl.CR ), CR = ctl.CR ; end + if isfield (ctl, "strategy") && ! isnan (ctl.strategy), + strategy = ctl.strategy ; end + if isfield (ctl, "refresh") && ! isnan (ctl.refresh), + refresh = ctl.refresh ; end + if isfield (ctl, "VTR") && ! isnan (ctl.VTR ), VTR = ctl.VTR ; end + if isfield (ctl, "tol") && ! isnan (ctl.tol ), tol = ctl.tol ; end + if isfield (ctl, "maxnfe") && ! isnan (ctl.maxnfe) + maxnfe = ctl.maxnfe; + end + if isfield (ctl, "maxiter") && ! isnan (ctl.maxiter) + maxiter = ctl.maxiter; + end + end +end + +## set dimension D and population size NP +D = length (XVmin); +if (NP == 0); + NP = 10 * D; +end + +## -------- do a few sanity checks -------------------------------- +if (length (XVmin) != length (XVmax)) + error("Length of upper and lower bounds does not match.") +end +if (NP < 5) + error("Population size NP must be bigger than 5.") +end +if ((F <= 0) || (F > 2)) + error("Difference Factor F out of range (0,2].") +end +if ((CR < 0) || (CR > 1)) + error("CR value out of range [0,1].") +end +if (maxiter <= 0) + error("maxiter must be positive.") +end +if (maxnfe <= 0) + error("maxnfe must be positive.") +end +refresh = floor(abs(refresh)); + +## ----- Initialize population and some arrays -------------------------- + +pop = zeros(NP,D); # initialize pop + +## pop is a matrix of size NPxD. It will be initialized with +## random values between the min and max values of the parameters +for i = 1:NP + pop(i,:) = XVmin + rand (1,D) .* (XVmax - XVmin); +end + +## initialise the weighting factors between 0.0 and 1.0 +w = rand (NP,1); +wi = w; + +popold = zeros (size (pop)); # toggle population +val = zeros (1, NP); # create and reset the "cost array" +bestmem = zeros (1, D); # best population member ever +bestmemit = zeros (1 ,D); # best population member in iteration +nfeval = 0; # number of function evaluations + +## ------ Evaluate the best member after initialization ------------------ + +ibest = 1; # start with first population member +val(1) = feval (fcn, [pop(ibest,:) const]); +bestval = val(1); # best objective function value so far +bestw = w(1); # weighting of best design so far +for i = 2:NP # check the remaining members + val(i) = feval (fcn, [pop(i,:) const]); + if (val(i) < bestval) # if member is better + ibest = i; # save its location + bestval = val(i); + bestw = w(i); + end +end +nfeval = nfeval + NP; +bestmemit = pop(ibest,:); # best member of current iteration +bestvalit = bestval; # best value of current iteration + +bestmem = bestmemit; # best member ever + +## ------ DE - Minimization --------------------------------------- +## popold is the population which has to compete. It is static +## through one iteration. pop is the newly emerging population. + +bm_n= zeros (NP, D); # initialize bestmember matrix in neighbourh. +lpm1= zeros (NP, D); # initialize local population matrix 1 +lpm1= zeros (NP, D); # initialize local population matrix 2 +rot = 0:1:NP-1; # rotating index array (size NP) +rotd= 0:1:D-1; # rotating index array (size D) + +iter = 1; +while ((iter < maxiter) && (nfeval < maxnfe) && (bestval > VTR) && ... + ((abs (max (val) - bestval) / max (1, abs (max (val))) > tol))) + popold = pop; # save the old population + wold = w; # save the old weighting factors + + ind = randperm (4); # index pointer array + + a1 = randperm (NP); # shuffle locations of vectors + rt = rem (rot + ind(1), NP); # rotate indices by ind(1) positions + a2 = a1(rt+1); # rotate vector locations + rt = rem (rot + ind(2), NP); + a3 = a2(rt+1); + rt = rem (rot +ind(3), NP); + a4 = a3(rt+1); + rt = rem (rot + ind(4), NP); + a5 = a4(rt+1); + + pm1 = popold(a1,:); # shuffled population 1 + pm2 = popold(a2,:); # shuffled population 2 + pm3 = popold(a3,:); # shuffled population 3 + w1 = wold(a4); # shuffled weightings 1 + w2 = wold(a5); # shuffled weightings 2 + + bm = repmat (bestmemit, NP, 1); # population filled with the best member + # of the last iteration + bw = repmat (bestw, NP, 1); # the same for the weighting of the best + + mui = rand (NP, D) < CR; # mask for intermediate population + # all random numbers < CR are 1, 0 otherwise + + if (strategy > 6) + st = strategy - 6; # binomial crossover + else + st = strategy; # exponential crossover + mui = sort (mui'); # transpose, collect 1's in each column + for i = 1:NP + n = floor (rand * D); + if (n > 0) + rtd = rem (rotd + n, D); + mui(:,i) = mui(rtd+1,i); #rotate column i by n + endif + endfor + mui = mui'; # transpose back + endif + mpo = mui < 0.5; # inverse mask to mui + + if (st == 1) # DE/best/1 + ui = bm + F*(pm1 - pm2); # differential variation + elseif (st == 2) # DE/rand/1 + ui = pm3 + F*(pm1 - pm2); # differential variation + elseif (st == 3) # DE/target-to-best/1 + ui = popold + F*(bm-popold) + F*(pm1 - pm2); + elseif (st == 4) # DE/best/2 + pm4 = popold(a4,:); # shuffled population 4 + pm5 = popold(a5,:); # shuffled population 5 + ui = bm + F*(pm1 - pm2 + pm3 - pm4); # differential variation + elseif (st == 5) # DE/rand/2 + pm4 = popold(a4,:); # shuffled population 4 + pm5 = popold(a5,:); # shuffled population 5 + ui = pm5 + F*(pm1 - pm2 + pm3 - pm4); # differential variation + else # DEGL/SAW + ## The DEGL/SAW method is more complicated. + ## We have to generate a neighbourhood first. + ## The neighbourhood size is 10% of NP and at least 1. + nr = max (1, ceil ((0.1*NP -1)/2)); # neighbourhood radius + ## FIXME: I don't know how to vectorise this. - if possible + for i = 1:NP + neigh_ind = i-nr:i+nr; # index range of neighbourhood + neigh_ind = neigh_ind + ((neigh_ind <= 0)-(neigh_ind > NP))*NP; + # do wrap around + [x, ix] = min (val(neigh_ind)); # find the local best and its index + bm_n(i,:) = popold(neigh_ind(ix),:); # copy the data from the local best + neigh_ind(nr+1) = []; # remove "i" + pq = neigh_ind(randperm (length (neigh_ind))); + # permutation of the remaining ind. + lpm1(i,:) = popold(pq(1),:); # create the local pop member matrix + lpm2(i,:) = popold(pq(2),:); # for the random point p,q + endfor + ## calculate the new weihting factors + wi = wold + F*(bw - wold) + F*(w1 - w2); # use DE/target-to-best/1/nocross + # for optimisation of weightings + ## fix bounds for weightings + o = ones (NP, 1); + wi = sort ([0.05*o, wi, 0.95*o],2)(:,2); # sort and take the second column + ## fill weighting matrix + wm = repmat (wi, 1, D); + li = popold + F*(bm_n- popold) + F*(lpm1 - lpm2); + gi = popold + F*(bm - popold) + F*(pm1 - pm2); + ui = wm.*gi + (1-wm).*li; # combine global and local part + endif + ## crossover + ui = popold.*mpo + ui.*mui; + + ## enforce initial bounds/constraints if specified + if (constr == 1) + for i = 1:NP + ui(i,:) = max (ui(i,:), XVmin); + ui(i,:) = min (ui(i,:), XVmax); + end + end + + ## ----- Select which vectors are allowed to enter the new population ------ + for i = 1:NP + tempval = feval (fcn, [ui(i,:) const]); # check cost of competitor + if (tempval <= val(i)) # if competitor is better + pop(i,:) = ui(i,:); # replace old vector with new one + val(i) = tempval; # save value in "cost array" + w(i) = wi(i); # save the weighting factor + + ## we update bestval only in case of success to save time + if (tempval <= bestval) # if competitor better than the best one ever + bestval = tempval; # new best value + bestmem = ui(i,:); # new best parameter vector ever + bestw = wi(i); # save best weighting + end + end + endfor #---end for i = 1:NP + + nfeval = nfeval + NP; # increase number of function evaluations + + bestmemit = bestmem; # freeze the best member of this iteration for the + # coming iteration. This is needed for some of the + # strategies. + + ## ---- Output section ---------------------------------------------------- + + if (refresh > 0) + if (rem (iter, refresh) == 0) + printf ('Iteration: %d, Best: %8.4e, Worst: %8.4e\n', ... + iter, bestval, max(val)); + for n = 1:D + printf ('x(%d) = %e\n', n, bestmem(n)); + end + end + end + + iter = iter + 1; +endwhile #---end while ((iter < maxiter) ... + +## check that all variables are well within bounds/constraints +boundsOK = 1; +for i = 1:NP + range = XVmax - XVmin; + if (ui(i,:) < XVmin + 0.01*range) + boundsOK = 0; + end + if (ui(i,:) > XVmax - 0.01*range) + boundsOK = 0; + end +end + +## create the convergence result +if (bestval <= VTR) + convergence = 1; +elseif (abs (max (val) - bestval) / max (1, abs (max (val))) <= tol) + convergence = 0; +elseif (boundsOK == 0) + convergence = -1; +elseif (iter >= maxiter) + convergence = -2; +elseif (nfeval >= maxnfe) + convergence = -3; +end + +endfunction + +%!function result = f(x); +%! result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2; +%!test +%! tol = 1.0e-4; +%! ctl.tol = 0.0; +%! ctl.VTR = 1.0e-6; +%! ctl.XVmin = [-2 -2]; +%! ctl.XVmax = [ 2 2]; +%! rand("state", 11) +%! [x, obj_value, nfeval, convergence] = de_min (@f, ctl); +%! assert (convergence == 1); +%! assert (f(x) == obj_value); +%! assert (obj_value < ctl.VTR); + +%!demo +%! ## define a simple example function +%! f = @(x) peaks(x(1), x(2)); +%! ## plot the function to see where the minimum might be +%! peaks() +%! ## first we set the region where we expect the minimum +%! ctl.XVmin = [-3 -3]; +%! ctl.XVmax = [ 3 3]; +%! ## and solve it with de_min +%! [x, obj_value, nfeval, convergence] = de_min (f, ctl)