nT = 6;
U = [-1, -1, 0, 0, 0, 0; 0, 0, -1, -1, 0, 0; 0, 0, 0, 0, -1, -1];
u = [-1; -2; -4];
v = [0; 0; 0];

% Initialize the pdf to zero
pdf = zeros(nT,1);

% Initialize the indices to keep in the problem to everything
DecidedIndex   = find(pdf < 0);
UnDecidedIndex = find(pdf == 0);

% Start the loop to construct the pdf
while(length(DecidedIndex) < nT)

   % Try to solve this puppy
   LogMesg('    -> solving for piece of pdf');
   [m,n] = size(U(:,UnDecidedIndex));
   Soln = CplexLink([zeros(length(UnDecidedIndex),1); 1 ], ...
      [U(:,UnDecidedIndex), v ; eye(n), -ones(n,1) ], ...
      [u - U*pdf; zeros(n,1)], [], [], ...
      zeros(length(UnDecidedIndex)+1,1), [], [], 'dsimplex');
Soln.status.value
pause
Soln.dual
Soln.statustext
Soln.variables'
Soln.objective
   % Need to embedd this soln into the overall solution to keep indices right
   TmpSol = pdf;
   TmpSol(UnDecidedIndex) = Soln.variables(1:length(Soln.variables)-1);
   TmpDual = zeros(6,1);
   TmpDual(UnDecidedIndex) = Soln.dual(m+1:length(Soln.dual));
TmpDual
pause
   % Find where the max is and is not attained
   PrimalDeltaDecidedIndex = find( ...
      (TmpSol >= Soln.variables(length(Soln.variables)) - 10^(-5)) & ...
      (TmpSol <= Soln.variables(length(Soln.variables)) + 10^(-5)) );
   DualDeltaDecidedIndex = find( (TmpDual > 10^(-5)) | (TmpDual < -10^(-5)) );
   DeltaDecidedIndex = ...
      intersect(PrimalDeltaDecidedIndex, DualDeltaDecidedIndex);
   DecidedIndex = union(DecidedIndex, DeltaDecidedIndex);
   UnDecidedIndex = setdiff(UnDecidedIndex, DeltaDecidedIndex);
PrimalDeltaDecidedIndex
DualDeltaDecidedIndex
DeltaDecidedIndex
DecidedIndex
UnDecidedIndex
pause

   % Update the pdf
   pdf(DeltaDecidedIndex) = Soln.objective;
   
end;

save /tmp/pdfsolve pdf
