function VQ(Prob, N, PDFType, prescription, SaveType);

%
% Make the Dir. for VQ stuff
%
if ~isdir('../images/VQ')
   !mkdir ../images/VQ
end

for n = 1:length(N)

for s = 1:length(fieldnames(PDFType));

% grab system clock
t0 = clock;

%
% Put the PDf in a local vec
%
angles = eval(['Prob.', eval(sprintf('PDFType.t%d',s)) ] );

if max(angles) > .01

% Calculate the accumulative dose and find the break
% points of each interval

AccAng(1) = angles(1);
NumBrkPts = N(n);
TotDose = angles*ones(Prob.NumAngles,1);
AngleIncr = 360/Prob.NumAngles;
BrkPtIntvl = TotDose/NumBrkPts;
NxtBrkPt = BrkPtIntvl;
BrkPtsX = zeros(NumBrkPts+1,1);
BrkPtsY = zeros(NumBrkPts+1,1);
ActualBrkPtsY = zeros(NumBrkPts+1,1);
BrkPtsX(1) = 0;
BrkPtsY(1) = 0;
ActualBrkPtsY(1) = 0;
k = 2;
for i = 2:Prob.NumAngles
   AccAng(i) = AccAng(i-1) + angles(i);
   if AccAng(i) > NxtBrkPt
      BrkPtsX(k) = i;
      ActualBrkPtsY(k) = AccAng(i);
      BrkPtsY(k) = NxtBrkPt;
      k = k + 1;
      NxtBrkPt = NxtBrkPt + BrkPtIntvl;
   end
end
BrkPtsX(NumBrkPts+1) = Prob.NumAngles;
BrkPtsY(NumBrkPts+1) = TotDose;

% Find the Centriods

% Initialize the Centriod vectors.  IndxCentriods is the vector of centriods
% in NumAngle `space'.  So if there are 24 angles, equally spaced every
% 15 degrees, IndxCentriods Centroids over [1, 2, ..., 24].  Centroids is
% the vector of real angle Centriods -i.e. over the cont. interval [0, 360).
IndxCentroids = zeros(NumBrkPts,1);
Centroids = IndxCentroids;
IndxBrkPtsX = Centroids;
Pos = 1;
for i = 2:length(BrkPtsX)
   Step = ((BrkPtsX(i) - BrkPtsX(i-1)) ) - 1;
   indx = Pos : Pos + Step;
   Pos = Pos + Step + 1;
   IndxBrkPtsX(i-1) = Pos - 1;
   % These are the centriods in the NumAnlges index
   if ~isempty(indx)
      IndxCentroids(i-1) = ...
         ((indx)*AccAng(indx)') / ...
         (ones(1,length(AccAng(indx)))*AccAng(indx)');
   else
      IndxCentroids(i-1) = IndxCentroids(i-2);
   end
end

%
% Clean things up by making sure we have integers
%
IndxCentroids = round(IndxCentroids);
Centroids = (IndxCentroids - 1) * AngleIncr;

%
%Make some nice figures
%
% This figure shows the dose over the anlge index, and labels
% the centroids.
figure 
plot(1:Prob.NumAngles, angles)
hold on
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kh')
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kx')
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kp')
hold off
saveas(gcf, ...
   ['../images/VQ/',eval([sprintf('PDFType.t%d', s)]),'-Ctrds'], 'epsc')

% This figure shows the accumulative dose over the angle index and
% labels the centriods.
if min([10^4, ceil(TotDose)]) == 10^4
   fprintf('    -> TotDose is %d, which exceeds 10^4 in VQ PDFType.t%d \n', ...
            ceil(TotDose), s)
   PlotEnd = 10^4;
else
   PlotEnd = ceil(TotDose);
end
figure
plot(1:Prob.NumAngles, AccAng)
hold on
for i = 1:length(BrkPtsY)
   plot(1:Prob.NumAngles, ActualBrkPtsY(i) * ones(1, Prob.NumAngles))
   plot( BrkPtsX(i) * ones(1, PlotEnd) , 1 : PlotEnd )
end
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kh')
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kx')
plot(IndxCentroids, zeros(length(IndxCentroids),1), 'kp')
saveas(gcf, ...
   ['../images/VQ/',eval([sprintf('PDFType.t%d', s)]),'-CumDose-Ctrds'], 'epsc')

TmpProb = RestrictedSolver(Prob, prescription, IndxCentroids);

Figure(TmpProb, prescription, ...
   sprintf('%s_NumBeams%d',eval(sprintf('PDFType.t%d',s)), N(n)),'VQ', SaveType)

% Write a line to the latex file
LaTex('AddRow', 'Vector Quantizer', eval(sprintf('PDFType.t%d',s)), IndxCentroids, etime(clock,t0)/60 );

else

fprintf('    -> PDF is zero, skipping VQ \n')

end

end

end
