%
% Parse - A. Holder and J. Reese, 6/4/03
%
% ver. 1.0
%
% This function takes the dose matrix A, the prescription, and the 
% uniformity level and parses the information to form a problem instance.
%

function [Prob] = Parse(A, prescription, UniLevel, ColIndex);

%
% This defines our problem structure.  A problem instance is defined 
% by the dosematrix (and its submatrices, as well as the bounds of 
% the problem, and the solution. We also indicate the solver used 
% and some other useful numbers.
%

Prob = struct('DoseMatrix', [], 'TumorDoseMatrix', [], ...
   'CriticalDoseMatrix', [], 'RegularDoseMatrix', [], ...
   'TumorRowIndex', [], 'CriticalRowIndex', [], ...
   'RegularRowIndex', [], 'TumorColIndex', [], ...
   'TGoal', [], 'TUB', [], 'TLB', [], 'CUB', [], 'RUB', [], ...
   'AnalysisType', [], 'SolverType', [], 'omega', [], ...
   'tol', [], 'x', [], 'fval', [], 'alpha', [], 'beta', [], ...
   'gamma', [], 'ThePlan', [], 'RadLevel', [], 'Angles', [], ...
   'NumAngles', [], 'NumPencils', [], 'BalancedSoln', [], ...
   'BalancedAvg', [], 'BalancedMax', [], 'PrimalAvg', [], ...
   'PrimalMax', [], 'DualAvg', [], 'DualMax', [], ...
   'InteriorAvg', [], 'InteriorMax', [], 'SetCover1', [], ...
   'SetCover2', [], 'SetCover3', [], 'Scoring', [], ...
   'EntropyAvgBalanced', []);

%
% Set the angle geometry
%
LogMesg('    -> Assuming an angle every 5 degrees')
Prob.NumAngles  = 72;
Prob.Angles = [1:72]';
LogMesg('    -> NumPencils is variable array')
Prob.NumPencils = ColIndex;

%
% Grab some indices
% 0 - Air; 1 - normal; 2 - target; 3+ critical
%
Prob.RegularRowIndex = find(prescription(:, : , :, 2) == 1);
Prob.TumorRowIndex = find(prescription(:, : , :, 2) == 2);
Prob.CriticalRowIndex = find(prescription(:, : , :, 2) >= 3);

%
% Partition the dose deposition matrix A into those rows corresponding
% to tumor pixels, critical structure pixels, and regular tissue pixels.
%

[s1, s2] = size(A);
Prob.DoseMatrix         = A;
Prob.TumorDoseMatrix    = sparse(A(Prob.TumorRowIndex,:));
if length(Prob.CriticalRowIndex) >= 1
  Prob.CriticalDoseMatrix = sparse(A(Prob.CriticalRowIndex,:));
  Prob.CUB   = prescription(Prob.CriticalRowIndex);
elseif length(Prob.CriticalRowIndex) == 0
  Prob.CriticalDoseMatrix = sparse(zeros(1,s2));
  Prob.CUB   = [0];
else
  LogMesg('RAD -> Confused about Prob.CricitalRowIndex')
end

Prob.RegularDoseMatrix  = sparse(A(Prob.RegularRowIndex,:));

%
% Grab some sizes.
%

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

%
% Remove the columns corresponding to the sub-beams that do not
% strike the tumor.
%

CheckColSum             = ones(1,tr)*Prob.TumorDoseMatrix;
Prob.TumorColIndex      = find(CheckColSum >= 0);
Prob.TumorDoseMatrix    = sparse(Prob.TumorDoseMatrix(:,Prob.TumorColIndex));
Prob.CriticalDoseMatrix = sparse(Prob.CriticalDoseMatrix(:,Prob.TumorColIndex));
Prob.RegularDoseMatrix  = sparse(Prob.RegularDoseMatrix(:,Prob.TumorColIndex));

%
% Extract the upper and lower bounds from the prescription matrix
%
% Prob.tol - the amount that we are going to increase the tumor lower
%            bound

Prob.tol   = 10^(-4);
p          = 3 / 100;
Prob.TGoal = prescription(Prob.TumorRowIndex);
Prob.TUB   = Prob.TGoal * (1 + p);
Prob.TLB   = (Prob.TGoal * (1 - p)) + Prob.tol;
Prob.RUB   = prescription(Prob.RegularRowIndex);

% Hard code the anlysis type
LogMesg('    -> Analysis Type is set to Absolute')
Prob.AnalysisType = 1;

%
% Calculate the weight omega that guarantees the tumor uniformity.
%
% The bound presented in the paper leads to numerical instability
% with the interior point method always thinking that it is infeasible.
% This value seems to work well.
%

kappa      = (norm(Prob.TLB, inf));
Prob.omega = kappa / Prob.tol;

% Set the Solver Type
LogMesg('    -> SolverType set to Interior, may not matter')
Prob.SolverType = 0;

