Code covered by the BSD License  

Highlights from
OpenOpt

image thumbnail

OpenOpt

by

 

25 Nov 2006 (Updated )

nonSmoothSolve (similar to fsolve), non-smooth & noisy local + some global solvers; works in Octave

ooRunProbSolver(prob, solver)
function r = ooRunProbSolver(prob, solver)
global globstat graphicsInfo


if ~isempty(globstat)
%     disp('parallel task mode isn''t implemented yet')
%     disp('clearing globstat...')
%     disp('restarting...')
    globstat = [];
end

globstat.isFinished =  0;

globstat.nFunEvals = 0;
globstat.nGradEvals = 0;
globstat.nHessEvals = 0;

globstat.nCEvals = 0;
globstat.nDCEvals = 0;

globstat.nHEvals = 0;
globstat.nDHEvals = 0;

globstat.xPrevF = [];
globstat.xPrevC = [];
globstat.xPrevH = [];

globstat.nIter = 0;
globstat.ff = inf;
globstat.xf = nan*ones(prob.n, 1);
if prob.graphics.drawingInOneWindow && ~isempty(graphicsInfo); globstat.graphics = graphicsInfo;end
if ~isfield(globstat, 'graphics'); globstat.graphics=[]; end
fn = fieldnames(prob.graphics);
for i=1:length(fn)
    globstat.graphics = setfield(globstat.graphics, fn{i}, getfield(prob.graphics, fn{i})); %#ok<SFLD,GFLD>
end
prob0 = prob;

if ~iscell(prob.f)
    if ischar(prob.f); prob.primal.fName = prob.f;
    else prob.primal.fName = func2str(prob.f); end
else
    prob.primal.fName = 'UnnamedFunc';
end

if nargin>1
    switch class(solver)
        case 'function_handle'
            prob.solver = solver;
            prob.solverName = func2str(solver);
        case 'char'
            prob.solverName = solver;
            prob.solver = str2func(solver);
        otherwise
            prob.err('incorrect solver is chosen, must be func name (string) or func handler')
    end
end

isNonSmoothSolve = isfield(prob, 'TolFunSolve') || isequal(lower(prob.solverName), 'nonsmoothsolve');
prob.primal.isUC = isUC(prob);%is unconstraint

fn = {'f' 'df' 'd2f' 'c' 'dc' 'd2c' 'h' 'dh' 'd2h'};%todo: d2L, dL?
for i = 1:length(fn)
    func = fn{i};
    if isfield(prob, func) && ~isempty(getfield(prob, func)) %#ok<GFLD>
        funcVal = getfield(prob,func); %#ok<GFLD>
        switch class(funcVal)
            case {'function_handle' 'function handle' 'cell'}% 2nd - for Octave
%                 prob.Fun{i} = funcVal;
                prob.primal = setfield(prob.primal, func, funcVal); %#ok<SFLD>
%                 funcNargins(end+1) = nargin(funcVal); %#ok<AGROW>
            case 'char'
%                 prob.primal = setfield(prob.primal, func, str2func(funcVal)); %#ok<SFLD>
                %don't replace both str2func by one - it yields error in current
                %MATLAB ver!!
%                 prob.Fun{i} = str2func(funcVal);
                prob.primal = setfield(prob.primal, func, str2func(funcVal)); %#ok<SFLD>
%                 funcNargins(end+1) = nargin(str2func(funcVal)); %#ok<AGROW>
            otherwise
                prob.err(['unknown or incorrect class "' class(funcVal) '"of objective function descriptor'])
        end
        
        %todo: do something for Octave to be able get nargin(FunHandlers)
%         ff = prob.Fun{i};
%         if iscell(ff); ff = ff{1};end
%         if isequal(prob.env, 'matlab')
%             nArgIn = nargin(ff);
%         elseif isequal(prob.env, 'octave')
%             nArgIn = nargin(func2str(ff));%func2str because Octave can't handle nargin(handlers) for now
%         else
%             prob.err('unknown environment')
%         end
%         switch nArgIn
%             case 1
%                 %for Octave compability
% %                 eval(['prob.primal = setfield(prob.primal, func , @(x, prob) prob.Fun{' int2str(i) ' }(x));']); %#ok<SFLD>
%             case 2
%                 %for Octave compability
% %                 eval(['prob.primal = setfield(prob.primal,  func , @(x, prob) prob.Fun{' int2str(i) ' }(x, prob));']); %#ok<SFLD>
%             otherwise
%                 prob.err('please, use fun(x) or fun(x,prob) syntax only; pass your parameters via prob.user field')
%         end
        %todo: replace upper by lower when Octave will be compatible
%         switch nargin(prob.Fun{i})
%             case 1
%                 prob.primal = setfield(prob.primal, func, @(x, prob) prob.Fun{i}(x)); %#ok<SFLD>
%             case 2
%                 prob.primal = setfield(prob.primal, func, @(x, prob) prob.Fun{i}(x,prob)); %#ok<SFLD>
%             otherwise
%                 prob.err('please, use fun(x) or fun(x,prob) syntaxis only; pass your parameters via prob.user field')
%         end
        eval(['prob.' func ' = @oo_' func ';']);
    else
        eval(['prob.' func ' = [];']);
        eval(['prob.primal.' func ' = [];']);
    end
end
% if any(funcNargins>2 | (~funcNargins))
%     prob.err('please, use fun(x) or fun(x,prob) syntaxis only; pass your parameters via prob.user field')
% end
% if any(funcNargins<0)
%     prob.warn()
% end



if isinf(prob.MaxIter); globstat.iterFvals = zeros(1e4, 1); else
    globstat.iterFvals = zeros(prob.MaxIter,1); end

if prob.advanced.traceIterX
    if isinf(prob.MaxIter); globstat.iterX = zeros(prob.n,1e4); else
         globstat.iterX = zeros(prob.n, prob.MaxIter); end
end

if prob.advanced.traceF
    if isinf(prob.MaxFunEvals); globstat.F = zeros(1e4, 1); else
         globstat.F = zeros(prob.MaxFunEvals,1); end
end
if prob.advanced.traceX
    if isinf(prob.MaxIter); globstat.X = zeros(prob.n, 1e4); else
         globstat.X = zeros(prob.n, prob.MaxFunEvals); end
end



if prob.iterPrint>=0
    prob.info(['problem ' prob.name ': solver ' prob.solverName ' starting...'])
end
globstat.tmp.dispBestVal = prob.primal.isUC && ~isequal(lower(prob.solverName), 'nonsmoothsolve');
if  isfield(prob.advanced, 'dispBestVal') && prob.advanced.dispBestVal 
    globstat.tmp.dispBestVal = 1;
end



% this field is designed for Base Module, NOT for solvers
% objective func substitution will be used authomatically

prob.b = prob.b(:);
prob.beq = prob.beq(:);
prob.lb = prob.lb(:);
prob.ub = prob.ub(:);


for fn = {'x0' 'ScaleFactor' 'TypicalX' 'lb' 'ub'}
    if ~isfield(prob, fn{1}); continue; end
    f = getfield(prob, fn{1}); %#ok<GFLD>
    if size(f, 2)>1 % this condition handles empty f as well
        prob.warn([fn{1} ' must be N x 1 vector, reshaping will be used, but errors can be obtained'])
        prob = setfield(prob, fn{1}, f(:)); %#ok<SFLD>
    end
end

for fn = {'A' 'Aeq' 'b' 'beq' 'lb' 'ub'}
    prob.primal = setfield(prob.primal, fn{1}, getfield(prob, fn{1})); %#ok<GFLD,SFLD>
end

if prob.useScaling
    if ~isfield(prob, 'ScaleFactor')
        prob.info('setting scale factor')
        if isfield(prob, 'TypicalX')
            prob.ScaleFactor = prob.TypicalX;
        else
            prob.ScaleFactor = prob.x0;
        end
        ind = find(~prob.ScaleFactor);
        if length(ind) == prob.n
            prob.ScaleFactor = ones(prob.n,1);
        else
            indm = find(abs(prob.ScaleFactor) == min(abs(prob.ScaleFactor)));
            prob.ScaleFactor(ind) = prob.ScaleFactor(indm(1));
        end
    else
        prob.ScaleFactor = prob.ScaleFactor(:);
    end
    prob.ScaleFactor = prob.ScaleFactor(:);
    prob.x0 = prob.x0 ./ prob.ScaleFactor;
    if ~isempty(prob.lb); prob.lb = prob.lb ./ prob.ScaleFactor; end
    if ~isempty(prob.ub); prob.ub = prob.ub ./ prob.ScaleFactor; end    
    if ~isempty(prob.A) 
        for j=1:size(prob.A, 1)
            prob.A(j,:) = prob.A(j,:) .* prob.ScaleFactor.';
        end
    end    
    if ~isempty(prob.Aeq) 
        for j=1:size(Aeq, 1)
            prob.Aeq(j,:) = prob.Aeq(j,:) .* prob.ScaleFactor;
        end
    end
%     prob.decode = @(x) x .* prob.ScaleFactor;
else
%     prob.decode = @(x) x;
end

if strcmp(prob.goalType, 'maximum')
%     prob.f0 = -prob.f0;% empty prob.f0 is handled correctly
    if ~isinf(prob.fEnough)% handling default value -inf
        prob.fEnough = - prob.fEnough;
    end
end

if isfield(prob, 'x0')
    f0_faled = 0;
%     try
        prob.f0 = prob.f(prob.x0, prob);%todo:handle exeptions or nan
%     catch
%         f0_faled=1;
%     end
    if isnan(prob.f0); f0_faled=1; end
    if f0_faled
        r.authors = '';
        r.alg = '';
        r.lastChanger = '';
        r.istop = -2;
        r.msg = msgstop(r.istop);
        r.ff = nan;
        r.xf = nan*ones(prob.n,1);
        return
    end
end



prob.nh = 0;
prob.nc = 0;
prob.primal.nc = 0;
prob.primal.nh = 0;
% inequality constraints
if isfield(prob, 'c') && ~isempty(prob.c)
    if ~isfield(prob, 'c0') || isempty(prob.c0)
        prob.c0 = prob.c(prob.x0, prob);
    end
    prob.nc = length(prob.c0);
    prob.primal.nc = prob.nc;
end


% equality constraints
if isfield(prob, 'h') && ~isempty(prob.h) 
    if ~isfield(prob, 'h0') || isempty(prob.h0)
        prob.h0 = prob.h(prob.x0, prob);
    end
    prob.nh = length(prob.h0);
    prob.primal.nh = prob.nh;
end

ooCheck(prob);% todo: skip if the flag is 0


prob.advanced.finiteLB = find(isfinite(prob.lb));
prob.advanced.finiteUB = find(isfinite(prob.ub));

globstat.advanced.drawTime=0;
globstat.advanced.drawCPUTime=0;
globstat.iterTime = 0;
globstat.iterCPUTime = 0;
prob.advanced.cputimeStart = cputime;
prob.advanced.timeStart = now;
r = prob.solver(prob);%<- ENGINE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r.CPUTimeElapsed = cputime-prob.advanced.cputimeStart - globstat.advanced.drawCPUTime;
r.TimeElapsed = 24*60*60*(now-prob.advanced.timeStart) - globstat.advanced.drawTime;
r.xf = r.xf(:);
globstat.isFinished = 1;
r.ff = prob.f(r.xf, prob);% for more safety
r.f = r.ff;% for more safety

violations = prob.advanced.getViolations(r.xf, prob);%need r.xf without scaling
if isempty(violations); violations=0; end
r.maxConstrViolation = norm(violations, inf);
r.isFeasible = r.maxConstrViolation<=prob.TolCon;

prob.iterfcn(prob, r);

if prob.useScaling
    prob.useScaling = 0;
    r.xf = r.xf .* prob.ScaleFactor;
end

if ~isfield(r, 'istop') && ~isfield(globstat, 'istop')
    prob.warn('no .istop field of output structure, setting to 0');
    r.istop = 0;
elseif isfield(globstat, 'istop')
    r.istop = globstat.istop;
end

if ~isfield(r, 'xf')
    if isfield(globstat, 'xf'); r.xf = globstat.xf;
    else
        prob.warn(['no .xf field in solver ' prob.solverName ' output structure; replacing by NaNs'])
        r.xf = nan*ones(prob.n, 1);
    end
end
if ~isfield(r, 'ff')
    if isfield(globstat, 'ff'); r.ff = globstat.ff;
    else
        prob.warn(['no .ff field in solver ' prob.solverName ' output structure; replacing by NaN'])
        r.ff = nan;
    end
end

if isNonSmoothSolve
    globstat.nFunEvals = globstat.nHEvals;
    if isfield(prob, 'TolFunSolve')
        prob.TolCon = prob.TolFunSolve;
        r.prob.TolCon = prob.TolFunSolve;
        r.prob.TolFun = prob.TolFunSolve;
    end
    if ~isempty(prob.primal.h); FUN = prob.primal.h;
    else FUN = prob.primal.f; end
    if iscell(FUN)
        if length(prob.primal.f)>1
            r.ff = zeros(length(FUN),1);
            for i=1:length(FUN)
                r.ff(i) = FUN{i}(r.xf, prob.varargin{:});
            end
        else
            r.ff = FUN{1}(r.xf, prob.varargin{:});
        end
    else
        r.ff = FUN(r.xf, prob.varargin{:});
    end
    
    r.maxConstrViolation = norm(r.ff, inf);
end

% if ~r.isFeasible; 
r = ooCheckSolution(prob, r);
% end

globstat.iterFvals(globstat.nIter+1:end)=[];

if prob.advanced.traceIterX
    globstat.iterX(prob.n, globstat.nIter+1:end)=[];
end

if prob.advanced.traceF
      globstat.F(globstat.nFunEvals+1:end)=[];
end
if prob.advanced.traceX
     globstat.X(prob.n, globstat.nFunEvals+1:end)=[];
end




if prob.debug
    prob.useScaling = 0;
    if norm(r.ff - prob.f(r.xf,prob))>1e-8
        prob.warn('check me!')
    end
end

if strcmp(prob.goalType, 'maximum')
    r.ff = -r.ff;
end

graphicsInfo = globstat.graphics;


% fn = setdiff(fieldnames(globstat), fieldnames(r)); 
% don't work in Octave yet
% check me in future
fn = fieldnames(globstat);
for i = 1:length(fn)
    if ~isfield(r, fn{i});
        r = setfield(r, fn{i}, getfield(globstat, fn{i})); %#ok<SFLD,STFLD,GFLD>
    end
end

globstat = [];

r.advanced.iterFvals = r.iterFvals;
r.advanced.iterTime = r.iterTime(:);
r.advanced.iterCPUTime = r.iterCPUTime(:);


for fi = {'inner' 'iterCPUTime' 'iterTime' 'iterFvals' 'legend' 'ylimprev' 'figureHandler' 'fPrev' 'lastX' 'x' 'f' 'isFinished' 'graphics' 'xPrevH' 'hPrevH' 'xPrevC' 'cPrevC' 'xPrevF' 'fPrevF' 'tmp'}
    if isfield(r, fi{1}); r = rmfield(r,fi{1}); end
end

if isempty(prob.primal.df); r = rmfield(r,'nGradEvals'); end
if isempty(prob.primal.d2f); r = rmfield(r,'nHessEvals'); end
if isempty(prob.primal.c); r = rmfield(r,'nCEvals'); end
if isempty(prob.primal.h); r = rmfield(r,'nHEvals'); end
if isempty(prob.primal.dc); r = rmfield(r,'nDCEvals'); end
if isempty(prob.primal.dh); r = rmfield(r,'nDHEvals'); end

if ~isfield(r,'istop'); r.istop = 0; end
r.msg = msgstop(r.istop);
r.prob = prob0;
if ~isfield(r, 'solver'); r.solver = prob.solverName; end

r.randInfo = ooRandInfo(prob);

%to prettify answer 
%I can't use orderfields because it requires exact number of known fields
for fn = {'nIter' 'CPUTimeElapsed' 'TimeElapsed' 'solver' 'istop' 'msg' 'isFeasible' 'xf' 'ff' 'randInfo'}
    tmp = getfield(r, fn{1}); %#ok<GFLD>
    r = rmfield(r, fn{1});
    r = setfield(r, fn{1}, tmp); %#ok<SFLD>
end

if prob.iterPrint>=0
    prob.info(['problem ' prob.name ': solver ' prob.solverName ' finished'])
    if r.istop>=0 || r.istop <= -1; sf = int2str(r.istop);
    else sf = num2str(r.istop); end
    prob.info(['stop flag: ' sf '   (' r.msg ')'])
    s = [];
    if ~isUC(prob); 
        if r.isFeasible; s = 'Solution is Feasible'; else s = 'Solution is Infeasible'; end
    end
    if isNonSmoothSolve
        r.maxConstrViolation = max(r.ff);
        prob.info(['max residual = ' num2str(r.maxConstrViolation)]);
    else
        if ~isempty(s); prob.info([s  ' (' 'max(abs(constr)) = ' num2str(r.maxConstrViolation) ')']); end
        prob.info(['Optim ObjFunVal = ' num2str(r.ff)]);
    end
        prob.info(['rand info: ' r.randInfo]);
end

% if prob.iterPrint>=0
%     prob.info(['problem ' prob.name ': solver ' prob.solverName ' finished'])
%     if r.istop>=0 || r.istop <= -1; sf = int2str(r.istop);
%     else sf = num2str(r.istop); end
%     prob.info(['stop flag: ' sf '   (' r.msg ')'])
%     s = [];
%     if ~isUC(prob); 
%         if r.isFeasible; s = 'Solution is Feasible'; else s = 'Solution is Infeasible'; end
%     end
%     if isNonSmoothSolve
%         prob.info(['max residual = ' num2str(r.maxConstrViolation)]);
%     else
%         if ~isempty(s); prob.info([s  ' (' 'max(abs(constr)) = ' num2str(r.maxConstrViolation) ')']); end
%         prob.info(['Optim ObjFunVal = ' num2str(r.ff)]);
%     end
%         prob.info(['rand info: ' r.randInfo]);
% end

Contact us