function [alpha]=minepsrt(MI,SG,alpha,dalpha)
% MINEPSrt finds optimal alpha (Generalized Anderson's task).
%  [alpha]=minepsrt(MI,SG,alpha,dalpha)
%
% MINEPSRT is an auxiliary function used in algorithms GANDERS
%  and GANDERS2. For more details refer to book SH10.
%
%  This function maximizes unimodal objective function which
%  acts in algorithm solving Generalized Anderson's task.
%
% See also MINEPS, MINEPSVL, GANDERS, GANDERS2, book SH10.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 11.5.2000
% Modifications
% 24. 6.00 V. Hlavac, comments polished.

% Attempt to improve epsilon calculation. It is needed to calculate
% roots of a polynomial of the degree 4.

K = size(MI,2);
N = size(MI,1);

% compute constants
for j = 1:K,
   s(j)= alpha'*MI(:,j);
   ss(j) = dalpha'*MI(:,j);
   ds(j) = ss(j) - s(j);
   sga(j) = alpha'*SG(:,(j-1)*N+1:j*N)*alpha;
   sgd(j) = dalpha'*SG(:,(j-1)*N+1:j*N)*dalpha;
   sgad(j) = dalpha'*SG(:,(j-1)*N+1:j*N)*alpha;
end


MAX_ROOTS=5;
allroots=[];      % roots
extr=[];          % extremes
for i=1:K-1,
   for j=i+1:K,

      %--------------
      a1=s(i); a2=s(j);
      b1=ds(i); b2=ds(j);
      c1=(sga(i)-2*sgad(i)+sgd(i)); c2=(sga(j)-2*sgad(j)+sgd(j));
      d1=(2*sgad(i)-2*sga(i)); d2=(2*sgad(j)-2*sga(j));
      e1=sga(i); e2=sga(j);

      A=(b1^2*c2 - b2^2*c1);
      B=(2*a1*b1*c2 + b1^2*d2 - 2*a2*b2*c1 - b2^2*d1 );
      C=(a1^2*c2 + 2*a1*b1*d2 + b1^2*e2 - a2^2*c1 - 2*a2*b2*d1 - b2^2*e1 );
      D=(a1^2*d2 + 2*a1*b1*e2 - a2^2*d1 - 2*a2*b2*e1 );
      E=(a1^2*e2 - a2^2*e1 );
      %------------
      % finds roots of the polynomial
      rts=roots([ A, B, C, D, E ])';

      % selects the real roots
      for k=1:size(rts,2),
         if isreal(rts(k)) & rts(k) > 0 & rts(k) < 1,
            allroots=[allroots,rts(k)];
         end
      end
   end

   % extreme
   t=(2*b1*e1-a1*d1)/(2*a1*c1-b1*d1);
   if t > 0 & t < 1,
      extr=[extr,t];
   end
end


% computes f(t)=min r(t)
posmin=[0,sort([extr,allroots]),1];
f=[];
for i=1:size(posmin,2),
   t=posmin(i);
   f=[f,min( (s+t*ds)./sqrt( (1-t)^2*sga + 2*t*(1-t)*sgad + t^2*sgd ) )];
end

% finds max min r
[maxf,inx]=max(f);
t=posmin(inx);

% compute new alpha
alpha=alpha*(1-t)+dalpha*t;