function [C,a,F] = Unmix (M, k)
% Code for tumor type unmixing based on the algorith of Ehrlich and Full (1987)
% Find k vertices bracketing M minimizing volume
% M is the expression matrix (rows=genes and columns=samples)
% k is a desired number of components
% V is the inferred components (rows=genes and columns=samples)
% a is the area of the simplex (used as a measure of fit quality)
% F is the inferred mixture fractions (rows=samples and columns=components)
% Code is copyright Russell Schwartz, Carnegie Mellon University 2009

% First perform PCA to get dimension down to k-1
[V,S]=princomp(M','econ');
P=S(:,1:k-1);
points=size(P,1);

% Initialize for the simplex search
[H,minVol]=convhulln(P);
minArea=0;
init=0;
i1p=0;
i2p=0;

% Conduct (points) trials of simplex search
for round=1:points
  % Initialize candidate vertex search
  cands=zeros(k,k-1);
  plist=zeros(k,1);
  pchosen=zeros(points,1);

  % Compute pairwise distances between points
  pairWeights=zeros((points+1)*(points+1),1);
  pairSum=0;
  for i1=1:(points-1)
   for i2=(i1+1):points
    pairWeights(i1*points+i2)=norm(P(i1,:)-P(i2,:))^k;
    pairSum=pairSum+pairWeights(i1*points+i2);
   end
  end

  % Sample a starting pair of points weighted by distance
  pairChoice=rand*pairSum;
  done=0;
  for i1=1:(points-1)
   if (done==1)
    break;
   end
   for i2=(i1+1):points
    if (pairChoice<pairWeights(i1*points+i2))
     cands(1,:)=P(i1,:);
     cands(2,:)=P(i2,:);
     plist(1)=i1;
     plist(2)=i2;
     pchosen(i1)=1;
     pchosen(i2)=1;
     i1p=i1;
     i2p=i2;
     done=1;
     break;
    else
     pairChoice=pairChoice-pairWeights(i1*points+i2);
    end
   end
  end

  % Pick points 3 through k weighted by volume
  for i=3:k

   % Find volume induced by each possible next point
   areaList=zeros(points,1);
   areaSum=0;
   for i1=1:points
    if (pchosen(i1)==0)
     cands(i,:)=P(i1,:);

     % Reduce dimension to i-1
     % Create orthogonalized matrix
     VM=cands(1:i,:);
     for i2=i:-1:1
      VM(i2,:)=VM(i2,:)-VM(1,:);
     end
     UM=VM(2:i,:);
     for i2=1:(i-1)
      if (norm(UM(i2,:))<1e-6)
       UM(i2,:)=zeros(size(UM(i2,:)));
      else
       UM(i2,:)=UM(i2,:)/norm(UM(i2,:));
      end
      for i3=(i2+1):(i-1)
       UM(i3,:)=UM(i3,:)-UM(i3,:)*UM(i2,:)'*UM(i2,:);
      end
     end
     ZM=(UM*VM')';

     % Measure volume
     [Ht,area]=convhulln(ZM,{'Pp','QJ'});

     % Check weight
     areaList(i1)=area^k;
     areaSum=areaSum+area^k;
    end
   end

   % Sample a next point
   areaChoice=rand*areaSum;
   for i1=1:points
    if (pchosen(i1)==0)
     if (areaChoice<areaList(i1))
      pchosen(i1)=1;
      plist(i)=i1;
      areaChoice=2*areaSum;
      cands(i,:)=P(i1,:);
      break;
     else
      areaChoice=areaChoice-areaList(i1);
     end
    end
   end
  end

  % Candidate set is now chosen
  % Next, move points outward to eliminate violations
  err=1;
  while (true)
   err=0;

   best=-inf;
   besti=0;

   % For each subset of k-1 points
   for i=1:k

    % Find a normal to the hyperplane defined by the k-1 points
    if (i==1)
     hplane=cands(2:k,:);
    elseif (i==k)
     hplane=cands(1:(k-1),:);
    else
     hplane=[cands(1:(i-1),:); cands(i+1:k,:)];
    end
    normA=ones(k-1,k-1);
    for j=2:k-1
     normA(j,:)=hplane(j,:)-hplane(1,:);
    end % for
    normb=zeros(k-1,1);
    normb(1,1)=1;
    normx=normA\normb;
    % check if this failed
    if (norm(normb-normA*normx)>1e-3)
     continue;
    end
    normx=normx/sqrt(normx'*normx);

    % Find direction from the plane to the kth point
    extra=cands(i,:);
    v=extra-hplane(1,:);
    sign=v*normx/abs(v*normx);
   
    % Find the farthest point outside the plane (if any)
    maxdist=-inf;
    argdist=0;
    violated=0;
    for j=1:points
     % Find the point
     pt=P(j,:);
     % Take a vector from the hyperplane
     vt=pt-hplane(1,:);
     % Project onto the normal
     dist=vt*normx;
     if (dist*sign<-1e-3)
      violated=violated+1;
     end
     if (-dist*sign>maxdist)
      maxdist=-dist*sign;
      argdist=dist;
     end
    end

    % See if we have a new farthest point
    if (violated>0 && violated>best)
     best=violated;
     besti=i;
     bestnorm=argdist*normx;
    end
    if (abs(maxdist)>err)
     err=abs(maxdist);
    end

   end % for

   % If no more violations then quit otherwise adjust and continue
   if (best<1)
    break;
   else
    % Translate the plane in the direction of the best normal to fix it
    for j=1:k
     if (j==besti)
      continue;
     end % if
     cands(j,:)=cands(j,:)+bestnorm';
    end % for
   end

  end % while

 % Test volume of simplex - keep if it's the best so far
 [H,area]=convhulln(cands,{'Pp','QJ'});
 if ((area<minArea && area>minVol) || (init==0))
  minArea=area;
  minK=cands;
  init=1;
 end % if

end % for

K=minK;
a=minArea;

% Explain each point as a convex combination of the inferred points
L=zeros(points*k,points*k);
b=zeros(points*k,1);
% Set up constraints sum_j f_{ij} c_{jk} = a_{ik}
for i=1:points
 for kp=1:(k-1)
  for j=1:k
   L((i-1)*(k-1)+kp,(i-1)*k+j)=K(j,kp);
  end
  b((i-1)*(k-1)+kp,1)=P(i,kp);
 end
end
% Set up constraints sum_j f_{ij} = 1
for i=1:points
 for j=1:k
  L(points*(k-1)+i,(i-1)*k+j)=1;
 end
 b(points*(k-1)+i,1)=1;
end
% Solve
ftmp=L\b;
F=zeros(points,k);
% Write point distributions to a matrix
for i=1:points
 for j=1:k
  F(i,j)=ftmp((i-1)*k+j,1);
 end
end

% Translate components from PCA space to gene expression space
C=V(:,1:(k-1))*K';
for i=1:k
 C(:,i)=C(:,i)+mean(M,2);
end

