function [alpha]=minepsvl(MI,SG,alpha,dalpha,tmax,vdelta)
% MINEPSVL finds optimal alpha (Generalized Anderson's task).
% [alpha]=minepsvl(MI,SG,alpha,dalpha,tmax,vdelta)
%
% 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.
%
%  This method discretizes function value of the objective
%  function.
%
% See also MINEPS, MINEPSRT, 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.

% default setting
if nargin < 5,
   tmax = inf;
   delta=1e-6;
elseif nargin < 6,
   delta=0;
end

% get dimension N and the # of distribution
K = size(MI,2);
N = size(MI,1);

% computes Fdown, Fup
minrs=(alpha'*MI)./sqrt( reshape(alpha'*SG,N,K)'*alpha )';
Fdown=min(minrs);
Fup=max(minrs);

% 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

[U]=fteps(Fdown,K,s,ss,ds,sga,sgd,sgad);
if isempty(U),
   % alpha can`t be better
   return;
end

T=U;
while (Fup-Fdown) > vdelta & tmax > 0,
   tmax=tmax-1;

   % compute Fmid
   Fmid=1/2*(Fup+Fdown);

   % find t from U, min r(alpha*(1-t)+dalpha*t,mi(j),sigma(j)) >= Fmid, for every j
   [U]=fteps(Fmid,K,s,ss,ds,sga,sgd,sgad);

   if isempty(U)==1,
      Fup=Fmid;
   else
      T=U;
      Fdown=Fmid;
   end

end % while ...

% computes alpha
t=1/2*(T(1)+T(2));
alpha=alpha*(1-t)+dalpha*t;

return

%%%%%%%%%%%%%%%%%%%%%%%%
function [U]=fteps(c,K,s,ss,ds,sga,sgd,sgad)
% checks if F(t) >= c, t is from U
%

Tup=1;
c2=c^2;
for j=1:K,
   A(j)=ds(j)^2 - sga(j)*c2 + 2*sgad(j)*c2 - sgd(j)*c2;
   B(j)=2*s(j)*ds(j) + 2*sga(j)*c2 - 2*sgad(j)*c2;
   C(j)=s(j)^2 - sga(j)*c2;

   %%%%%%%%%%%
   if ds(j) < 0,
      Tup=min([Tup, -s(j)/ds(j)]);
   end
end

U=[0 Tup];
for j=1:K,

   if A(j) == 0,
      % linear
      if B(j) ~= 0,
         t=-C(j)/B(j);
         if B(j) < 0, U=setand(U,[-inf t]); U=[]; else U=setand(U,[t inf]); end
      elseif C(j) < 0,
         U=[];
      end
   else
      % quadratic
      D=B(j)^2 - 4*A(j)*C(j);

      if D < 0,
         if A(j) < 0, U=[]; end
      elseif D == 0,
         t=-B(j)/(2*A(j));
         U=setand(U,[t t]);
      else % D > 0
         t1=(-B(j)-sqrt(D))/(2*A(j));
         t2=(-B(j)+sqrt(D))/(2*A(j));
         if t1 > t2, [t1,t2]=swap(t1,t2); end

         if A(j) < 0,
            U=setand(U,[t1 t2]);
         else
            U=setand(U,[-inf,t1,t2,inf]);
         end
      end % if , else
   end

   if isempty(U), return; end
end

return

function [T]=setand(T1,T2)
%
%

T=[];
for i=1:size(T1,2)/2,
   for j=1:size(T2,2)/2,
      U=and2([T1((i-1)*2+1),T1(i*2)],[T2((j-1)*2+1),T2(j*2)]);
      if ~isempty(U),
         T=[T,U];
      end
   end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [T]=and2(T1,T2)
% T= T1 U T2
%

T1 = sort(T1);
T2 = sort(T2);

if T1(1) > T2(1), [T1,T2]=swap(T1,T2); end

if T1(2) < T2(1),
   T=[];
elseif T2(2) < T1(2),
   T=T2;
else
   T=[T2(1) T1(2)];
end

return

function [A,B]=swap(A,B)
%
ASWAP=A;
A=B;
B=ASWAP;
return