%%%%%%%%%% adesc_minpolytope.m %%%%%%%%%%
function [Ptmin] = crmz_minpolytope(P);
% Finds the nearest point Ptmin in the convex hull of 
% the set of points P1,...,PM given by the column of 
% the N*M matrix P.
% Implementation of a method by P. Wolfe presented in 
% "Mathematical Programming", vol. 11, 1976, pp. 128-149. 
%
S = []; w = []; no_stp = 1;
[N,M] = size(P);
if (M==1)
  no_stp = 0;
  X = P;
end;
%Z1 = 10^(-12);
Z1 = 10^(-10);
Z2 = 10^(-10);
Z3 = 10^(-10);
P_norm = sum(abs(P).^2);
[pmn,J] = min(P_norm);
S = [S;J];
w = [w;1];
while (no_stp)
  X = P(:,S)*w;
  [pmn,J] = min(X'*P);
  PS_norm = max(sum(abs(P(:,S)).^2));
  PJ_norm = P(:,J)'*P(:,J);
  no_stp = (X'*P(:,J) <= (X'*X - Z1*max([PJ_norm;PS_norm])));
 if (no_stp)
    I = find(~(S-J));
    if (length(I) > 0)
      no_stp = 0;
      disp('temporary disaster: see Note 3')
    end;
 end;
 if (no_stp)
  S = [S;J];
  w = [w;0];
  flg = 1;
  while (flg)
    e = ones(size(S));
    A = ones(length(e),length(e)) + P(:,S)'*P(:,S);
    u = A \ e;
    v = u./(e'*u);
    I = find(v <= Z2);
    if ~(length(I))
      w = v;
      flg = 0;
    else
      I = find((w-v) > Z3);
      t = w(I)./(w(I)-v(I));
      theta = min([1;1-min(t)]);
      w = theta*w + (1-theta)*v;
      I = find(w <= Z2);
      if (length(I) > 0)
        w(I) = zeros(length(I),1);
        w(I(1)) = [];
        S(I(1)) = [];
      end;
    end; 
  end;
 end;
end;
Ptmin = X;
%%%%%%%%%% End adesc_minpolytope.m %%%%%%%%%%
