function MIPSolve(Prob, prescription);

%load /tmp/OncologyExample

%
% Grab some sizes.
%

[tr,tc] = size(Prob.TumorDoseMatrix);
[gr,gc] = size(Prob.CriticalDoseMatrix);
[rr,rc] = size(Prob.RegularDoseMatrix);

%
% Some numbers to add angle constraints to MIP:
% number of beams, bixels per beam, and matrix size
% 

H = Prob.NumAngles;
P = Prob.NumPencils;
M = tr+gr+rr;
N = tc;

%
% Set BigM for beam switch constraints and number of beams
% The number of beams should be part of the problem input
%

BigM = 1000;

DiffAngles = [5, 7, 10, 15];

for q = 1:length(DiffAngles)

NumBeams = DiffAngles(q);

%
% Construct constraint matrix for angle selection
% For each angle an identity matrix for the bixels of that angle
% and -bigM coefficient for angle variable
%

IntConst = zeros(P*H+1,(P+1)*H);
for i = 1:H
    NewBlock = [zeros(P,(i-1)*P), eye(P), zeros(P,(H-i)*P), zeros(P,i-1), ones(P,1)*(-BigM), zeros(P,H-i)];
IntConst((i-1)*P+1:i*P,1:(P+1)*H) = NewBlock;
end


%
% Row to restrict number of angles: all ones
%

IntConst(P*H+1,1:(P+1)*H) = [zeros(1,P*H), ones(1,H)];

%
% Remove columns corresponding to sub-beams not hitting the tumor
%

NeededCols = union(Prob.TumorColIndex,[H*P+1:H*(1+P)]);
IntConst = IntConst(:,NeededCols);

%
% Remove rows where sub-beam variable has been deleted, i.e. row sum > -BigM
%
 
[mi,ni]     = size(IntConst);
CheckRowSum = IntConst*ones(ni,1);
IntRows     = find(CheckRowSum > -BigM);
IntConst    = IntConst(IntRows,:);
[ir,ic] = size(IntRows);

%
% Create A matrix, rhs b, lower and upper bounds and variable type for solver
%

if Prob.AnalysisType == 0
 
  U  = sparse([ Prob.TumorDoseMatrix,    zeros(tr, H), zeros(tr, tr), zeros(tr, gr), zeros(tr, rr);
               -Prob.TumorDoseMatrix,    zeros(tr, H),      -eye(tr), zeros(tr, gr), zeros(tr, rr);
                Prob.CriticalDoseMatrix, zeros(gr, H), zeros(gr, tr),      -eye(gr), zeros(gr, rr);
                Prob.RegularDoseMatrix,  zeros(rr, H), zeros(rr, tr), zeros(rr, gr),      -eye(rr);
                IntConst,                              zeros(ir, tr), zeros(ir, gr), zeros(ir, rr)]);
  u  = sparse([Prob.TUB; -Prob.TLB; Prob.CUB; Prob.RUB; zeros(ir-1,1); NumBeams]);
  lb = sparse([zeros(tc, 1); zeros(H,1); zeros(tr, 1); -Prob.CUB; zeros(rr,1)]);
  ub = sparse([inf * ones(tc, 1); ones(H, 1); Prob.TLB; inf * ones(gr, 1); inf * ones(rr, 1)]);
  c  = sparse([zeros(tc,1); zeros(H, 1); Prob.omega * ones(tr, 1); ones(gr, 1); ones(rr, 1)]);
  VarType = [char(ones(1, tc)*'C'),char(ones(1, H)*'B'),char(ones(1,tr+gr+rr)*'C')];
 
elseif Prob.AnalysisType == 1
 
  U  = sparse([ Prob.TumorDoseMatrix,    zeros(tr, H), zeros(tr, 1), zeros(tr, 1), zeros(tr, 1);
               -Prob.TumorDoseMatrix,    zeros(tr, H), -ones(tr, 1), zeros(tr, 1), zeros(tr, 1);
                Prob.CriticalDoseMatrix, zeros(gr, H), zeros(gr, 1),  -ones(gr,1), zeros(gr, 1);
                Prob.RegularDoseMatrix,  zeros(rr, H), zeros(rr, 1), zeros(rr, 1), -ones(rr, 1);
 		IntConst,                              zeros(ir, 1), zeros(ir, 1), zeros(ir, 1)]);
  u  = sparse([Prob.TUB; -Prob.TLB; Prob.CUB; Prob.RUB; zeros(ir-1,1); NumBeams]);
  lb = sparse([zeros(tc, 1); zeros(H,1); 0; -min(Prob.CUB); 0]);
  ub = sparse([inf * ones(tc, 1); ones(H,1); min(Prob.TLB); inf; inf]);
  c  = sparse([zeros(tc,1); zeros(H,1); Prob.omega; 1; 1]);
  VarType = [char(ones(1,tc)*'C'),char(ones(1,H)*'B'),'CCC'];

else
 
 LogMesg('Analysis Type is Not Selected')
  return
 
end;

%
% Set Solver to CPLEX MIP
%

SolverType = 'mip';
 
%
% Grab the the clock
%
 
StartTime = clock;
 
%
% Off to the Solver
%
 
Soln = CplexLink(c, U, u, [], [], lb, ub, VarType, SolverType);

%
% Display optimal beam angles and solution value
%

Prob.x = Soln.variables;
angles = find(Prob.x(tc+1:tc+H)>0.9999)*(360/Prob.NumAngles);
[nb,nix]=size(angles);
fval = Soln.objective;

LogMesg(sprintf('Optimal value: %f3 ', fval));
LogMesg('Optimal angles: ')
     for i=1:nb 
         LogMesg(sprintf('%d ',angles(i)));
     end;
 
%
% Calculuate the time required by the solver.
%
 
TotalTime = etime(clock,StartTime);
 
%
% Partition the solution up into the different types of
% decision variables.
%
 
if Prob.AnalysisType == 0
                                                                                                                    
  Prob.alpha = Prob.x(tc+H+1:tc+tr+H);
  Prob.beta  = Prob.x(tc+tr+H+1:tc+tr+gr+H);
  Prob.gamma = Prob.x(tc+tr+gr+H+1:tc+tr+gr+rr+H);
  Prob.x     = Prob.x(1:tc);
 
 
 elseif Prob.AnalysisType == 1
 
  Prob.alpha = Prob.x(tc+H+1);
  Prob.beta  = Prob.x(tc+H+2);
  Prob.gamma = Prob.x(tc+H+3);
  Prob.x     = Prob.x(1:tc);
 
 end;
 
%
% Embedd the smaller plan, x, into a larger one and
% calculate the radiation levels over the image.
%
 
  [NumARows,NumACols] = size(Prob.DoseMatrix);
  Prob.ThePlan = zeros(NumACols,1);
  Prob.ThePlan(Prob.TumorColIndex) = Prob.x;
  Prob.RadLevel = Prob.DoseMatrix*Prob.ThePlan;
 
Figure(Prob);

end
