ralg(prob)

function r = ralg(prob)
%    Ukrainian Science Academy
%    Instityte of cybernetics (www.icyb.kiev.ua), optimization department
%    Contacts: 10-38-044-526-21-68,  e-mail: [e-mail]
%    Engine - r-algorithm of Ukrainean academician Naum Z. Shor
%   INPUT:
%     prob - structure, describing the problem
%     if prob.ralg has fields 'alp', 'h0', 'nh', 'q1', 'q2'
%     then they will override the default vals
%    Working arrays:
%      b(n,n) -- backward spase transformation matrix
%      g1(n),g2(n) -- gradients
%
%    if you have choise : implement nonlinear constraint as equality 
%    or inequality one, latter is more prefered.
%
%    Last modification: 28.1.2007     /Dmitrey Kroshko/
%   see ooexample2, ooexample5 & other for correct ralg usage
%   for small nVars=1...10 ShorEllipsoid can be better
%   OpenOpt also has nonSmoothSolve() based on ralg (default) - 
%   MATLAB/Octave fsolve analog for nonsmooth funcs



if isempty(prob.df); prob.df = @oo_df; end% forcify numerical (sub)gradient anyway
if ~isempty(prob.c) && isempty(prob.dc); prob.dc = @oo_dc; end% forcify numerical (sub)gradient anyway
if ~isempty(prob.h) && isempty(prob.dh); prob.dh = @oo_dh; end% forcify numerical (sub)gradient anyway

% if ~isUC(prob)%todo: if any equality constraints
if prob.nc || prob.nh || ~isempty(prob.A) || ~isempty(prob.Aeq) || ~all(isinf(prob.lb)) || ~all(isinf(prob.ub))
    prob.solver = @ralg;%for more safety - maybe, it was called from other solver, global or so
    r = nonlinconstr(prob); return
end


r = oor('authors: Dmitrey Kroshko, [e-mail]', ...
    'alg:Naum Z. Shor r-algorithm with AST & some modifications',...
    'lastChanger: Dmitrey', ...
    'ralg_info:Ukrainian Science Academy, Instityte of cybernetics (www.icyb.kiev.ua). Sometimes hand-turning ralg parameters can enhance convergence');

%setting default parameters
[alp h0 nh q1 q2 hsmin] = ...
    deal(2.0, 1.0, 3, 1.0, 1.1, max(0.00025, 100*prob.TolX));%defaults
% if prob.classF>0 && prob.classC>0% for functions with at least 1st derivatives
%     q1 = 0.95;
% end
n = prob.n;

needRej = @(b) norm(b(:), 2)/n^2  > 1e-20;

%checking are any solver parameters is structure prob.(name of solver) (here - prob.ralg) 
ExtractRoutineParamsFromProb;%script-file
%if present - they will replace defaults 

%for economy of memory
if ~exist('b', 'var'); b = diag(ones(n,1)); end
% sometimes b can be obtained from ExtractRoutineParamsFromProb()
if ~exist('hs', 'var'); hs = h0; end

w = 1.0/alp-1.0;


% if prob.classF<0 || prob.classC<0
%     prob.warn('ralg can''t properly handle discontinious problems!')
% end

% if prob.needGlobal
%     prob.warn('ralg can''t handle global extremum problems, local extremum will be found')
% end

x0 = prob.x0;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Shor r-alg engine                         %
% translated from fortran with some changes %
% 19-Sep-2006 11:12:25                      %
% //by Dmitrey Kroshko                      %
% //icq 275976670(inviz)                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = x0;
xPrev = x0;
r.xf = x0;
r.x = x0;

f = prob.f0;
r.f = f;
r.ff = f;

g = prob.df(x, prob);

r.df = g;
r = prob.iterfcn(prob, r);
if r.istop; 
    %r.advanced.ralg.hs = hs; 
    return
end

for itn = 1:prob.MaxIter+8
    ls1 = 0;
    g2 = b.' * g;
    if ~any(g2); g2(1) = 1e-15; end
    g2 = g2 ./ norm(g2);
    g1 = b * g2;

    x_0 = x;
    for ls=1:prob.MaxLineSearch
        ls1 = ls1 + 1;
        x = x - hs .* g1;

        if ls1 > nh
            hs = hs .* q2;
            ls1 = 0;
        end
        
        f = prob.f(x, prob);
        g2 = prob.df(x, prob);
        

        
        if f < r.ff
            r.ff = f;
            r.xf = x;
        end
        d = sum(g1.*g2);

        if d<=0
%             if ~any(g1); g1(1)=1e-11; end
%             if ~any(g2); g2(1)=1e-11; end
%             phi = acos(d/norm(g1)/norm(g2));
            break
        end % angle > pi/2
    end
%     disp(ls)
    
    %debug
%     if prob.advanced.debug
%         g3 = prob.df(x, prob);
%         f3 = prob.f(x, prob);
%         disp('checking numerical gradient')
%         df2 = zeros(prob.n,1);
%         prob.diffInt = 1e-6;
%         for i=1:prob.n
%             x(i)=x(i)+prob.diffInt;
%             df2(i) = (prob.f(x,prob) - f3)/prob.diffInt;
%             x(i)=x(i)-prob.diffInt;
%         end
%         
%         NN = norm(g3-df2,1);
%         if NN>1e-1 || any(prob.c(x,prob)>0)
%             disp([g3 df2 g3-df2]);
%             disp(['norm1 = ' num2str(NN)])
%             keyboard
%         end
%     end
    
    %debug end
    
    if ls == prob.MaxLineSearch; r.istop = -5; return; end


    if itn == 1 ... % 1st iter only!
            && ls==1
         while 1
%              disp(hs)
             if hs * 0.75 < hsmin; break; end
             hs = hs * 0.75;
             x_prev = x;
             f_prev = f;
             hs_prev = hs;
             x = x_0 - hs .* g1;
             f = prob.f(x, prob);
             g2 = prob.df(x, prob);
             if f < r.ff
                 r.ff = f;
                 r.xf = x;
             end
             d = sum(g1.*g2);
             if d>0 || f > f_prev
                 x = x_prev;
                 f = f_prev;
                 hs = hs_prev;
                 break
             end % angle > pi/2
         end
    end
    
%     if ls==1; hs=hs.*q1; end
    %CHANGES!
    q1 = 0.9;
    if ls == 1; hs = max(hsmin, hs .* q1); end
    %changes end


    g1 = g2 - g;
    g = b.' * g1;
    if ~any(g); g(1) = 1e-15; end
    ng = norm(g);
    if ng < prob.TolGrad; r.istop = 2; end
    g = g / ng;
%     disp(norm(b(:), 2) / n^2)
    if needRej(b)
        b = b + (b * g) * (g.' .* w);
    else
        b = diag(ones(n,1));
        hs = max(norm(xPrev - x), hsmin);
    end

    g = g2;
    r.df = g;%CHECK ME! - MODIFIED!!
    r.f = f;
    r.x = x;
%     if r.ff > r.f
%         r.ff = r.f;
%         r.xf = r.x;
%     end
    r = prob.iterfcn(prob, r);
    if r.istop
        r.advanced.ralg.hs = max(norm(xPrev - x), hsmin);
        r.df = g2;
        return
    end
    xPrev = x;
%     prob = updateLM(r.x, prob);
end