function [Data] = pvwcgpeak(M,blk,opts,PrevData)
%function [Data] = pvwcgpeak(M,blk,opts,PrevData)

% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: 1.1.6.2 $ $Date: 2004/08/25 19:12:02 $

muopt = 'cw';
starttime = clock;
GTH.M = opts.MGoodThreshold;
GTH.A = opts.AGoodThreshold;
BTH.M = opts.MBadThreshold;
BTH.A = opts.ABadThreshold;
ABSTOL = opts.AbsTol;
RELTOL = opts.RelTol;
MAXTIME = opts.MaxTime;
NTIMES = opts.NTimes;
MAXCNT = opts.MaxCnt;
%PFlag = opts.Display;
PFlag = 'off';

bidx = blk2idx(blk);
if bidx.repfull.num>0
   error('No repeated full blocks');
   return
end 
% Number of real parameters
nblk = bidx.sreal.num + bidx.repreal.num;
% [realDidx,imagDidx,realGidx,imagGidx,realD,imagD,realG,imagG] = mu2amidx(szm(1),blk);
[ftdidx,varcnt] = uball(M(:,:,1),bidx,[],[],1);
% [left_ends right_ends LOWER UPPER ACTIVE DGbeta]
ALLCOR = 1:(2*nblk);
LEFTCOR = 1:nblk;
RIGHTCOR = (nblk+1):(2*nblk);
LOWER = 2*nblk + 1;
UPPER = 2*nblk + 2;
ACTIVE = 2*nblk + 3;
% realblk = [ones(bidx.sreal.num,1);bidx.repreal.repeated];
% DGbeta = 1 + 2*sum(realblk.*realblk) + bidx.full.num; % beta
DGbeta = 1 + varcnt; % beta
SCALINGS = (2*nblk+4):(2*nblk+3+DGbeta);
cols = 2*nblk + 3 + DGbeta;
SCALARS = 1:bidx.sreal.num;
REPS = bidx.sreal.num+1:nblk;

if isempty(PrevData)
   row = bidx.rdimm+1:size(M,1);
	col = bidx.cdimm+1:size(M,2);
	npts = size(M,3);
   Sr = ones(1,npts);
   Sc = ones(1,npts);
	for i = 1:npts            
      Sr(i) = 1/max( norm(M(row,1:bidx.cdimm,i)) , 0.01 );
      Sc(i) = 1/max( norm(M(1:bidx.rdimm,col,i)) , 0.01 );
      rnorm = norm(M(row,1:bidx.cdimm,i));
      cnorm = norm(M(1:bidx.rdimm,col,i));
      if rnorm<0.01
         rnorm = 0.01;
      elseif rnorm>100
         rnorm = 100;
      end
      if cnorm<0.01
         cnorm = 0.01;
      elseif cnorm>100
         cnorm = 100;
      end
      Sr(i) = 1/rnorm;
      Sc(i) = 1/cnorm;
	end
else
   Sr = PrevData.Sr;
   Sc = PrevData.Sc;
end
%Uncomment next 2 lines for no scaling
%Sr = ones(1,npts);
%Sc = ones(1,npts);

S = Sr.*Sc;
for i = 1:npts
   M(row,:,i) = Sr(i)*M(row,:,i);
   M(:,col,i) = M(:,col,i)*Sc(i);
end

if isempty(PrevData)
   [PrevData] = localinit(M,nblk,blk,bidx,...
        ALLCOR,LOWER,UPPER,ACTIVE,SCALINGS,DGbeta,cols,ftdidx,S,...
        muopt,opts);
   % PrevData is a scalar structure
end
AllWsreal = PrevData.AllWsreal;
AllWrepreal = PrevData.AllWrepreal;
AllWlvec = PrevData.AllWlvec;
AllWrvec = PrevData.AllWrvec;
AllWrepcomp = PrevData.AllWrepcomp;
BBList = PrevData.BBList;
freqactive = PrevData.freqactive;
freqUmax = PrevData.freqUmax;
freqUmaxcube = PrevData.freqUmaxcube;       
freqLmax = PrevData.freqLmax;
nomcost = PrevData.nomcost;
npts = PrevData.npts;
bestL = PrevData.bestL;
bestU = PrevData.bestU;
SbestL = PrevData.SbestL;
SbestU = PrevData.SbestU;

reasonstostop = [all(freqUmax <= S*GTH.A + GTH.M*nomcost); ...
      any(freqLmax >= S*BTH.A + BTH.M*nomcost); ...
      SbestU-SbestL<ABSTOL; ...
      SbestU-SbestL<RELTOL*SbestU; ...
      etime(clock,starttime)>MAXTIME; ...
      strcmpi(opts.LowerBoundOnly,'On')];

extdim = 50;
frac = 0.5;
while all(reasonstostop==0) & bidx.sreal.num+bidx.repreal.num>0
	% Clean up: go through each frequency, and delete cubes with upper
	%   bounds that are lower than the current best LOWER bound
	for ii=1:length(freqactive)
       i = freqactive(ii);
       useless = find(BBList{i,1}(:,ACTIVE)==1 & BBList{i,1}(:,UPPER)<S(i)*SbestL);    
       % deletecube (cube j in frequency i)
       BBList{i,1}(useless,ACTIVE) = 0;
       % delete frequency
       if all(BBList{i,1}(:,ACTIVE)==0)
          BBList{i,2} = 0;
          freqactive(ii) = 0;
       end
	end
	freqactive(find(freqactive==0)) = [];
	% Find offending cube
	[SbestU,idx] = max(freqUmax./S);
	usefreq = idx(1);
	cubetosplit = freqUmaxcube(usefreq);

   % Reallocate BBList if necessary
   freespace = find(BBList{usefreq,1}(:,ACTIVE)==0);
   used = size(BBList{usefreq,1},1) - length(freespace);
   if length(freespace)==0
      freespace = size(BBList{usefreq,1},1)+(1:extdim);
      BBList{usefreq,1} = [BBList{usefreq,1};zeros(extdim,cols)];
   end

   % Array index is: usefreq
   % Cube to split is: cubetosplit 
   % Split cube, and duplicate into two cubes
   cubej = BBList{usefreq,1}(cubetosplit,:);
   [dum,dimidx] = max(cubej(RIGHTCOR)-cubej(LEFTCOR));
   k = dimidx(1);  % coordinate with longest side-length
   left = cubej(LEFTCOR(k));
   right = cubej(RIGHTCOR(k));
   mid = left + frac*(right-left);
   cubejleft = cubej;
   cubejright = cubej;
   cubejleft(RIGHTCOR(k)) = mid;
   cubejright(LEFTCOR(k)) = mid;

   % Recenter and normalize
   McubeL = rcnrs(M(:,:,usefreq),bidx,cubejleft(LEFTCOR),cubejleft(RIGHTCOR));
   McubeR = rcnrs(M(:,:,usefreq),bidx,cubejright(LEFTCOR),cubejright(RIGHTCOR));

   % Call lower bounds
   % Implicitly, we have [-1 , 1] assumed in the LOWERBOUND_CALC
   MXGAIN = BTH.M*nomcost(usefreq) + S(usefreq)*BTH.A;
   [lowerL,baddeltaLsreal,baddeltaLrepreal,...
     baddeltaLlvec,baddeltaLrvec,baddeltaLrepcomp,il] = ...
     wclowc(McubeL,bidx,NTIMES,MAXCNT,0,MXGAIN);
   [lowerR,baddeltaRsreal,baddeltaRrepreal,...
     baddeltaRlvec,baddeltaRrvec,baddeltaRrepcomp,ir] = ...
     wclowc(McubeR,bidx,NTIMES,MAXCNT,0,MXGAIN);
   
   activeLower = cubej(LOWER);
   activeUpper = cubej(UPPER);
   
   % New lower bounds (hopefully at least one is larger, though not necessarily)
   cubejleft(1,LOWER) = lowerL;
   cubejright(1,LOWER) = lowerR;
        
   % Call upper bounds, using subdivided cube's scalings.
   % Implicitly, we have [-1 , 1] assumed in the UPPERBOUND_CALC
   if isinf(cubej(UPPER))
      [upperL,di,dbs,gzr,gzl,psdami,dgbL] = uball(McubeL,bidx);
      %BBList{cubetosplit,3} = {di dbs gzr gzl};
      [upperR,di,dbs,gzr,gzl,psdami,dgbR] = uball(McubeR,bidx);
      %BBList{freespace(1),3} = {di dbs gzr gzl};
   else
      [upperL,di,dbs,gzr,gzl,psdami,dgbL] = uball(McubeL,bidx,cubej(SCALINGS),1);
      %BBList{cubetosplit,3} = {di dbs gzr gzl};
      [upperR,di,dbs,gzr,gzl,psdami,dgbR] = uball(McubeR,bidx,cubej(SCALINGS),1);
      %BBList{freespace(1),3} = {di dbs gzr gzl};
   end
   if ~isinf(upperL)
      cubejleft(SCALINGS) = dgbL(:)';
      cubejleft(1,UPPER) = upperL;
   end
   if ~isinf(upperR)
      cubejright(SCALINGS) = dgbR(:)';
      cubejright(1,UPPER) = upperR;
   end

   % New upper bounds (hopefully both are smaller, algorithm may
   %   not yield this, but it usually will.  We probably should not accept
   %   numbers when they get worse - just use the old ones -- but for now
   cubejleft(1,UPPER) = upperL;
   cubejright(1,UPPER) = upperR;
   
   Wgain = max(freqLmax./S);
   
   % Check if either is worst, unnormalize if so.
   if (lowerL/S(usefreq)>=Wgain) & (lowerL>=lowerR)
      offset = (cubejleft(RIGHTCOR)+cubejleft(LEFTCOR))/2;
      slope = (cubejleft(RIGHTCOR)-cubejleft(LEFTCOR))/2;
      AllWsreal{usefreq} = offset(SCALARS) + ([baddeltaLsreal(:).']).*slope(SCALARS);
      AllWrepreal{usefreq} = offset(REPS) + ([baddeltaLrepreal(:).']).*slope(REPS);
      AllWlvec{usefreq} = baddeltaLlvec;
      AllWrvec{usefreq} = baddeltaLrvec;
      AllWrepcomp{usefreq} = baddeltaLrepcomp;
   elseif (lowerR/S(usefreq)>=Wgain) & (lowerR>=lowerL)
      offset = (cubejright(RIGHTCOR)+cubejright(LEFTCOR))/2;
      slope = (cubejright(RIGHTCOR)-cubejright(LEFTCOR))/2;
      AllWsreal{usefreq} = offset(SCALARS) + ([baddeltaRsreal(:).']).*slope(SCALARS);
      AllWrepreal{usefreq} = offset(REPS) + ([baddeltaRrepreal(:).']).*slope(REPS);
      AllWlvec{usefreq} = baddeltaRlvec;
      AllWrvec{usefreq} = baddeltaRrvec;
      AllWrepcomp{usefreq} = baddeltaRrepcomp;
   elseif lowerR/S(usefreq)<Wgain & lowerL/S(usefreq)<Wgain
      %disp('Splitting did not help lower bound');
   end
   % All info is ready in new cubes, so put this into BBList.  We overwrite
   %   the one we split, as well as add to the BBList at FREESPACE(1)        
   BBList{usefreq,1}(cubetosplit,:) = cubejleft;
   BBList{usefreq,1}(freespace(1),:) = cubejright;

   % Update freq{LU}max
	lookat = find(BBList{usefreq,1}(:,ACTIVE)==1);
	[mxl,idxl] = max(BBList{usefreq,1}(lookat,LOWER));
	[mxu,idxu] = max(BBList{usefreq,1}(lookat,UPPER));
   % mxl,usefreq,size(freqLmax)
	freqLmax(usefreq) = mxl;
	freqUmax(usefreq) = mxu;
   
   freqUmaxcube(usefreq) = lookat(idxu(1));            
   bestL = max(freqLmax);
   bestU = max(freqUmax);
   SbestL = max(freqLmax./S);
   SbestU = max(freqUmax./S);

   if strcmpi(PFlag,'On')
      fprintf('Frequency %d, NewUpper %g, NewLower %g\n',usefreq,SbestU,SbestL);
      fprintf('   Initial Lower %g, LowLeft %g, LowRight %g\n',...
         activeLower,lowerL,lowerR);
      fprintf('   Initial Upper %g, UppLeft %g, UppRight %g\n',...
         activeUpper,upperL,upperR);
   end

      % Clean up: go through each frequency, and delete cubes with upper
      %   bounds that are lower than the current best LOWER bound
      voluse = 0;
      mxvol = 0;
      for ii=1:length(freqactive)
         i = freqactive(ii);
         useless = find(BBList{i,1}(:,ACTIVE)==1 & BBList{i,1}(:,UPPER)/S(i)<Wgain);    
         % deletecube (cube j in frequency i)
         BBList{i,1}(useless,ACTIVE) = 0;
         % delete frequency
         if all(BBList{i,1}(:,ACTIVE)==0)
            BBList{i,2} = 0;
            freqactive(ii) = 0;
         else
            actfin = find(BBList{i,1}(:,ACTIVE)==1);
            sidelens = BBList{i,1}(actfin,RIGHTCOR) - BBList{i,1}(actfin,LEFTCOR);
            tmpvol = sum(prod(sidelens'));
            if tmpvol > mxvol
               mxvol = tmpvol;
            end
            voluse = voluse + tmpvol;
         end
      end
      %volrat = voluse/initvol;
      %mxvolrat = mxvol/(initvol/npts);
      freqactive(find(freqactive==0)) = [];

      reasonstostop = [all(freqUmax <= S*GTH.A + GTH.M*nomcost); ...
         any(freqLmax >= S*BTH.A + BTH.M*nomcost); ...
         SbestU-SbestL<ABSTOL; ...
         SbestU-SbestL<RELTOL*SbestU; ...
         etime(clock,starttime)>MAXTIME];
      
end
Data.AllWsreal = AllWsreal;
Data.AllWrepreal = AllWrepreal;
Data.AllWlvec = AllWlvec;
Data.AllWrvec = AllWrvec;
Data.AllWrepcomp = AllWrepcomp;
Data.BBList = BBList;
Data.freqactive = freqactive;
Data.freqUmax = freqUmax;
Data.SfreqUmax = freqUmax./S;
Data.freqUmaxcube = freqUmaxcube;       
Data.freqLmax = freqLmax;
Data.SfreqLmax = freqLmax./S;
Data.nomcost = nomcost;
Data.npts = npts;
Data.bestL = max(freqLmax);
Data.bestU = max(freqUmax);
Data.SbestL = max(freqLmax./S);
Data.SbestU = max(freqUmax./S);
Data.Sr = Sr;
Data.Sc = Sc;