Dawn@4: function [x,fval,exitflag,output]=fminsearchbnd(fun,x0,LB,UB,options,varargin)
Dawn@4: % FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
Dawn@4: % usage: x=FMINSEARCHBND(fun,x0)
Dawn@4: % usage: x=FMINSEARCHBND(fun,x0,LB)
Dawn@4: % usage: x=FMINSEARCHBND(fun,x0,LB,UB)
Dawn@4: % usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
Dawn@4: % usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
Dawn@4: % usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
Dawn@4: % 
Dawn@4: % arguments:
Dawn@4: %  fun, x0, options - see the help for FMINSEARCH
Dawn@4: %
Dawn@4: %  LB - lower bound vector or array, must be the same size as x0
Dawn@4: %
Dawn@4: %       If no lower bounds exist for one of the variables, then
Dawn@4: %       supply -inf for that variable.
Dawn@4: %
Dawn@4: %       If no lower bounds at all, then LB may be left empty.
Dawn@4: %
Dawn@4: %       Variables may be fixed in value by setting the corresponding
Dawn@4: %       lower and upper bounds to exactly the same value.
Dawn@4: %
Dawn@4: %  UB - upper bound vector or array, must be the same size as x0
Dawn@4: %
Dawn@4: %       If no upper bounds exist for one of the variables, then
Dawn@4: %       supply +inf for that variable.
Dawn@4: %
Dawn@4: %       If no upper bounds at all, then UB may be left empty.
Dawn@4: %
Dawn@4: %       Variables may be fixed in value by setting the corresponding
Dawn@4: %       lower and upper bounds to exactly the same value.
Dawn@4: %
Dawn@4: % Notes:
Dawn@4: %
Dawn@4: %  If options is supplied, then TolX will apply to the transformed
Dawn@4: %  variables. All other FMINSEARCH parameters should be unaffected.
Dawn@4: %
Dawn@4: %  Variables which are constrained by both a lower and an upper
Dawn@4: %  bound will use a sin transformation. Those constrained by
Dawn@4: %  only a lower or an upper bound will use a quadratic
Dawn@4: %  transformation, and unconstrained variables will be left alone.
Dawn@4: %
Dawn@4: %  Variables may be fixed by setting their respective bounds equal.
Dawn@4: %  In this case, the problem will be reduced in size for FMINSEARCH.
Dawn@4: %
Dawn@4: %  The bounds are inclusive inequalities, which admit the
Dawn@4: %  boundary values themselves, but will not permit ANY function
Dawn@4: %  evaluations outside the bounds. These constraints are strictly
Dawn@4: %  followed.
Dawn@4: %
Dawn@4: %  If your problem has an EXCLUSIVE (strict) constraint which will
Dawn@4: %  not admit evaluation at the bound itself, then you must provide
Dawn@4: %  a slightly offset bound. An example of this is a function which
Dawn@4: %  contains the log of one of its parameters. If you constrain the
Dawn@4: %  variable to have a lower bound of zero, then FMINSEARCHBND may
Dawn@4: %  try to evaluate the function exactly at zero.
Dawn@4: %
Dawn@4: %
Dawn@4: % Example usage:
Dawn@4: % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
Dawn@4: %
Dawn@4: % fminsearch(rosen,[3 3])     % unconstrained
Dawn@4: % ans =
Dawn@4: %    1.0000    1.0000
Dawn@4: %
Dawn@4: % fminsearchbnd(rosen,[3 3],[2 2],[])     % constrained
Dawn@4: % ans =
Dawn@4: %    2.0000    4.0000
Dawn@4: %
Dawn@4: % See test_main.m for other examples of use.
Dawn@4: %
Dawn@4: %
Dawn@4: % See also: fminsearch, fminspleas
Dawn@4: %
Dawn@4: %
Dawn@4: % Author: John D'Errico
Dawn@4: % E-mail: woodchips@rochester.rr.com
Dawn@4: % Release: 4
Dawn@4: % Release date: 7/23/06
Dawn@4: 
Dawn@4: % size checks
Dawn@4: xsize = size(x0);
Dawn@4: x0 = x0(:);
Dawn@4: n=length(x0);
Dawn@4: 
Dawn@4: if (nargin<3) || isempty(LB)
Dawn@4:   LB = repmat(-inf,n,1);
Dawn@4: else
Dawn@4:   LB = LB(:);
Dawn@4: end
Dawn@4: if (nargin<4) || isempty(UB)
Dawn@4:   UB = repmat(inf,n,1);
Dawn@4: else
Dawn@4:   UB = UB(:);
Dawn@4: end
Dawn@4: 
Dawn@4: if (n~=length(LB)) || (n~=length(UB))
Dawn@4:   error 'x0 is incompatible in size with either LB or UB.'
Dawn@4: end
Dawn@4: 
Dawn@4: % set default options if necessary
Dawn@4: if (nargin<5) || isempty(options)
Dawn@4:   options = optimset('fminsearch');
Dawn@4: end
Dawn@4: 
Dawn@4: % stuff into a struct to pass around
Dawn@4: params.args = varargin;
Dawn@4: params.LB = LB;
Dawn@4: params.UB = UB;
Dawn@4: params.fun = fun;
Dawn@4: params.n = n;
Dawn@4: params.OutputFcn = [];
Dawn@4: 
Dawn@4: % 0 --> unconstrained variable
Dawn@4: % 1 --> lower bound only
Dawn@4: % 2 --> upper bound only
Dawn@4: % 3 --> dual finite bounds
Dawn@4: % 4 --> fixed variable
Dawn@4: params.BoundClass = zeros(n,1);
Dawn@4: for i=1:n
Dawn@4:   k = isfinite(LB(i)) + 2*isfinite(UB(i));
Dawn@4:   params.BoundClass(i) = k;
Dawn@4:   if (k==3) && (LB(i)==UB(i))
Dawn@4:     params.BoundClass(i) = 4;
Dawn@4:   end
Dawn@4: end
Dawn@4: 
Dawn@4: % transform starting values into their unconstrained
Dawn@4: % surrogates. Check for infeasible starting guesses.
Dawn@4: x0u = x0;
Dawn@4: k=1;
Dawn@4: for i = 1:n
Dawn@4:   switch params.BoundClass(i)
Dawn@4:     case 1
Dawn@4:       % lower bound only
Dawn@4:       if x0(i)<=LB(i)
Dawn@4:         % infeasible starting value. Use bound.
Dawn@4:         x0u(k) = 0;
Dawn@4:       else
Dawn@4:         x0u(k) = sqrt(x0(i) - LB(i));
Dawn@4:       end
Dawn@4:       
Dawn@4:       % increment k
Dawn@4:       k=k+1;
Dawn@4:     case 2
Dawn@4:       % upper bound only
Dawn@4:       if x0(i)>=UB(i)
Dawn@4:         % infeasible starting value. use bound.
Dawn@4:         x0u(k) = 0;
Dawn@4:       else
Dawn@4:         x0u(k) = sqrt(UB(i) - x0(i));
Dawn@4:       end
Dawn@4:       
Dawn@4:       % increment k
Dawn@4:       k=k+1;
Dawn@4:     case 3
Dawn@4:       % lower and upper bounds
Dawn@4:       if x0(i)<=LB(i)
Dawn@4:         % infeasible starting value
Dawn@4:         x0u(k) = -pi/2;
Dawn@4:       elseif x0(i)>=UB(i)
Dawn@4:         % infeasible starting value
Dawn@4:         x0u(k) = pi/2;
Dawn@4:       else
Dawn@4:         x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
Dawn@4:         % shift by 2*pi to avoid problems at zero in fminsearch
Dawn@4:         % otherwise, the initial simplex is vanishingly small
Dawn@4:         x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
Dawn@4:       end
Dawn@4:       
Dawn@4:       % increment k
Dawn@4:       k=k+1;
Dawn@4:     case 0
Dawn@4:       % unconstrained variable. x0u(i) is set.
Dawn@4:       x0u(k) = x0(i);
Dawn@4:       
Dawn@4:       % increment k
Dawn@4:       k=k+1;
Dawn@4:     case 4
Dawn@4:       % fixed variable. drop it before fminsearch sees it.
Dawn@4:       % k is not incremented for this variable.
Dawn@4:   end
Dawn@4:   
Dawn@4: end
Dawn@4: % if any of the unknowns were fixed, then we need to shorten
Dawn@4: % x0u now.
Dawn@4: if k<=n
Dawn@4:   x0u(k:n) = [];
Dawn@4: end
Dawn@4: 
Dawn@4: % were all the variables fixed?
Dawn@4: if isempty(x0u)
Dawn@4:   % All variables were fixed. quit immediately, setting the
Dawn@4:   % appropriate parameters, then return.
Dawn@4:   
Dawn@4:   % undo the variable transformations into the original space
Dawn@4:   x = xtransform(x0u,params);
Dawn@4:   
Dawn@4:   % final reshape
Dawn@4:   x = reshape(x,xsize);
Dawn@4:   
Dawn@4:   % stuff fval with the final value
Dawn@4:   fval = feval(params.fun,x,params.args{:});
Dawn@4:   
Dawn@4:   % fminsearchbnd was not called
Dawn@4:   exitflag = 0;
Dawn@4:   
Dawn@4:   output.iterations = 0;
Dawn@4:   output.funcount = 1;
Dawn@4:   output.algorithm = 'fminsearch';
Dawn@4:   output.message = 'All variables were held fixed by the applied bounds';
Dawn@4:   
Dawn@4:   % return with no call at all to fminsearch
Dawn@4:   return
Dawn@4: end
Dawn@4: 
Dawn@4: % Check for an outputfcn. If there is any, then substitute my
Dawn@4: % own wrapper function.
Dawn@4: if ~isempty(options.OutputFcn)
Dawn@4:   params.OutputFcn = options.OutputFcn;
Dawn@4:   options.OutputFcn = @outfun_wrapper;
Dawn@4: end
Dawn@4: 
Dawn@4: % now we can call fminsearch, but with our own
Dawn@4: % intra-objective function.
Dawn@4: [xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);
Dawn@4: 
Dawn@4: % undo the variable transformations into the original space
Dawn@4: x = xtransform(xu,params);
Dawn@4: 
Dawn@4: % final reshape
Dawn@4: x = reshape(x,xsize);
Dawn@4: 
Dawn@4: % Use a nested function as the OutputFcn wrapper
Dawn@4:   function stop = outfun_wrapper(x,varargin);
Dawn@4:     % we need to transform x first
Dawn@4:     xtrans = xtransform(x,params);
Dawn@4:     
Dawn@4:     % then call the user supplied OutputFcn
Dawn@4:     stop = params.OutputFcn(xtrans,varargin{1:(end-1)});
Dawn@4:     
Dawn@4:   end
Dawn@4: 
Dawn@4: end % mainline end
Dawn@4: 
Dawn@4: % ======================================
Dawn@4: % ========= begin subfunctions =========
Dawn@4: % ======================================
Dawn@4: function fval = intrafun(x,params)
Dawn@4: % transform variables, then call original function
Dawn@4: 
Dawn@4: % transform
Dawn@4: xtrans = xtransform(x,params);
Dawn@4: 
Dawn@4: % and call fun
Dawn@4: fval = feval(params.fun,xtrans,params.args{:});
Dawn@4: 
Dawn@4: end % sub function intrafun end
Dawn@4: 
Dawn@4: % ======================================
Dawn@4: function xtrans = xtransform(x,params)
Dawn@4: % converts unconstrained variables into their original domains
Dawn@4: 
Dawn@4: xtrans = zeros(1,params.n);
Dawn@4: % k allows some variables to be fixed, thus dropped from the
Dawn@4: % optimization.
Dawn@4: k=1;
Dawn@4: for i = 1:params.n
Dawn@4:   switch params.BoundClass(i)
Dawn@4:     case 1
Dawn@4:       % lower bound only
Dawn@4:       xtrans(i) = params.LB(i) + x(k).^2;
Dawn@4:       
Dawn@4:       k=k+1;
Dawn@4:     case 2
Dawn@4:       % upper bound only
Dawn@4:       xtrans(i) = params.UB(i) - x(k).^2;
Dawn@4:       
Dawn@4:       k=k+1;
Dawn@4:     case 3
Dawn@4:       % lower and upper bounds
Dawn@4:       xtrans(i) = (sin(x(k))+1)/2;
Dawn@4:       xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
Dawn@4:       % just in case of any floating point problems
Dawn@4:       xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));
Dawn@4:       
Dawn@4:       k=k+1;
Dawn@4:     case 4
Dawn@4:       % fixed variable, bounds are equal, set it at either bound
Dawn@4:       xtrans(i) = params.LB(i);
Dawn@4:     case 0
Dawn@4:       % unconstrained variable.
Dawn@4:       xtrans(i) = x(k);
Dawn@4:       
Dawn@4:       k=k+1;
Dawn@4:   end
Dawn@4: end
Dawn@4: 
Dawn@4: end % sub function xtransform end
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: