ShorEllipsoid(prob)

function r = ShorEllipsoid(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 - modificated elipsoid method by Ukrainean academician Naum Z. Shor
%    Last modification: 15.1.2007     /Dmitrey Kroshko, icq 275976670(inviz)/
%    this solver is for small-scale problems only (nVars = 1...10)
%    elseware ralg() is preferable
%    requires r0 knowing
%    see ooexample*.m for ShorEllipsoid usage (similar to ralg)

%todo: prevent matrix B to be singular

r = oor('authors: Dmitrey Kroshko, [e-mail], icq 275976670', ...
    'alg:Naum Z. Shor modificated elipsoid method, adaptive space transformation',...
    'lastChanger: Dmitrey', ...
    'moreinfo:Ukrainian Science Academy, Instityte of cybernetics (www.icyb.kiev.ua), optimization department');


if isempty(prob.df); prob.df = @oo_df; end% forcify numerical (sub)gradient anywhere
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 prob.nc || prob.nh || ~isempty(prob.A) || ~isempty(prob.Aeq)
    prob.solver = @ShorEllipsoid;%for more safety - maybe, it was called from other solver, global or so
    r = nonlinconstr(prob); return
end

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

maxitn = prob.MaxIter;
x0 = prob.x0;

n = prob.n;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Shor Elipsoid engine                      %
% 11-Nov-2006 20:01:09                      %
% //by Dmitrey Kroshko                      %
% //icq 275976670(inviz)                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r.xf=x0;

if prob.useScaling;prob.err('ShorEllipsoid + scaling isn''t implemented yet');end
    
if prob.MaxIter<inf
    radius = zeros(prob.MaxIter,1);
else
    radius = zeros(3500,1);
end
if ~isfield(prob,'r0') || isempty(prob.r0); 
    prob.err('ShorEllipsoid: You must supply prob.r0 such that norm(x0-x_opt)<= r0')
else
    radius(1) = prob.r0;
end


B = diag(ones(n,1));
x=x0;


f = prob.f0;
g = prob.df(x, prob);%check for prob.df0

if any(isnan(g))
    prob.err('gradient(x0) has coordinate that equals to NaN')
    r.istop = -3;
    return
end

r.ff = f;
r.df = g;

beta = sqrt(1+1/n^2) - 1/n;

ExtractRoutineParamsFromProb;%script-file

for k=1:maxitn%not (while 1) because k value is used in alg
    r.f = f; r.x = x; r.df = g;
    if f<r.ff
        r.xf=x;
        r.ff=f; 
    end
    r = prob.iterfcn(prob, r);
    if r.istop; return; end

    dzeta = B' * g;
    dzeta = dzeta / sqrt(sum(dzeta.^2));
    hk = radius(k) / n * beta;
    x = x - hk * (B * dzeta);
    f = prob.f(x, prob);
    radius(k+1) = radius(k) * sqrt(1+1/n^2);
    B = B + (B * dzeta) * (dzeta' .* (1-1/beta));
    g = prob.df(x, prob);
end