function [pdf] = BalancedPDF(Prob)

%
% 12/22/04 - Make a Balanced pdf
%

% grab some sizes
[mT,nT] = size(Prob.TumorDoseMatrix);
[mC,nC] = size(Prob.CriticalDoseMatrix);
[mN,nN] = size(Prob.RegularDoseMatrix);

% This uses the optimal values of alpha, beta, and gamma to
% iteratively minimize the infinity norm of x.
% Depends on AnalysisType - 1 = absolute, 2 = average

fprintf('    -> Making data for Balanced pdf construction \n');
if Prob.AnalysisType == 1
   U = sparse( [ Prob.TumorDoseMatrix;
                -Prob.TumorDoseMatrix;
                 Prob.CriticalDoseMatrix;
                 Prob.RegularDoseMatrix ]);
   v = sparse( [ zeros(mT + mT + mC + mN, 1) ] );
   % Add a bit to the solution space for numerical stability
   fuzz = 10^(-4);
   u = sparse( [ Prob.TUB + fuzz;
                -Prob.TLB + Prob.alpha;
                 Prob.CUB + Prob.beta + fuzz;
                 Prob.RUB + Prob.gamma + fuzz] );
else
   fprintf('RAD -> ERROR: BalancedPDF only works with AnalysisType 1 (absolute analysis) \n')
end

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

% Initialize a tmp soln
OldSoln = Prob.x;
OldUnDecidedIndex = find(TmpPDF == 0);

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

% Start the loop to construct the pdf
StopIterations = 0;
while(StopIterations == 0)

   % Try to solve this puppy
   [m,n] = size(U(:,UnDecidedIndex));
   Soln = CplexLink([zeros(length(UnDecidedIndex),1); 1], ...
      [U(:,UnDecidedIndex), v ; eye(n), -ones(n,1)], ...
      [u - U*TmpPDF; zeros(n,1)], [], [], zeros(length(UnDecidedIndex)+1,1), ...
      [], [], 'dsimplex');

   if (Soln.status == 1)

   % Need to embedd this soln into the overall solution to keep indices right
   TmpSol = TmpPDF;
   TmpSol(UnDecidedIndex) = Soln.variables(1:length(Soln.variables)-1);
   TmpDual = zeros(nT,1);
   TmpDual(UnDecidedIndex) = Soln.dual(m+1:length(Soln.dual));

   % Update Old Soln in case next iteration dies
   OldSoln(UnDecidedIndex) = Soln.variables(1:length(Soln.variables)-1);
   OldUnDecidedIndex = UnDecidedIndex;

   % 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);

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

   % Lets check to see if we need to continue
   [gap, I] = min( [Prob.TUB - Prob.TumorDoseMatrix * TmpPDF; ...
                    Prob.CUB + Prob.beta - Prob.CriticalDoseMatrix * TmpPDF; ...
                    Prob.RUB + Prob.gamma - Prob.RegularDoseMatrix * TmpPDF ] );
   fprintf('\n')
%   if (gap / max(Prob.TUB)) < -0.02
%      fprintf('RAD -> We are terminating with a relative gap of %f \n', ...
%         (100*gap / max(Prob.TUB)) )
%      TmpPDF(UnDecidedIndex) = Soln.variables(1:length(Soln.variables)-1);
%      [NumARows,NumACols] = size(Prob.DoseMatrix);
%      pdf = zeros(NumACols,1);
%      pdf(Prob.TumorColIndex) = TmpPDF;
%      StopIterations = 1;
   if length(DeltaDecidedIndex) == 0
      fprintf('RAD -> BalancedPDf terminating, can not fine more variables to decided \n')
      TmpPDF(UnDecidedIndex) = Soln.variables(1:length(Soln.variables)-1);
      [NumARows,NumACols] = size(Prob.DoseMatrix);
      pdf = zeros(NumACols,1);
      pdf(Prob.TumorColIndex) = TmpPDF;
      StopIterations = 1;
   elseif length(DecidedIndex) == nT
      fprintf('RAD -> BalancedPDF terminating, all variables decided \n')
      pdf = zeros(NumACols,1);
      pdf(Prob.TumorColIndex) = TmpPDF;
      StopIterations = 1;
   else
      fprintf('RAD -> BalancedPDF continuing with a relative gap of %f \n', ...
         (100*gap / max(Prob.TUB)) )
      fprintf('    -> %d variable(s) fixed at %f, only %d left \n', ...
         length(DeltaDecidedIndex), Soln.objective, length(UnDecidedIndex) )
   end

   else

      fprintf('\nRAD -> WARNING in BlancedPDF:  Soln.status was %d. \n',Soln.status)
      fprintf('    -> WARNING in BalancedPDF:  Solver message was ')
      disp(Soln.statustext)
      fprintf('    -> Im returning the original solution as the pdf \n')
      TmpPDF(OldUnDecidedIndex) = OldSoln(OldUnDecidedIndex);
      [NumARows,NumACols] = size(Prob.DoseMatrix);
      pdf = zeros(NumACols,1);
      pdf(Prob.TumorColIndex) = TmpPDF;
      StopIterations = 1;
      pdf = Prob.ThePlan;
   end

end

return
