function [Data] = localinit(M,nblk,blk,bidx,...
   ALLCOR,LOWER,UPPER,ACTIVE,SCALINGS,DGbeta,cols,ftdidx,S,...
   muopt,opts);
% Copyright 2003-2004 The MathWorks, Inc.

predim = 5;%50;
CNTMAX = 4;
dims = ndims(M);
szm = size(M);
if dims==2
   npts = 1;
elseif dims==3
   npts = szm(3);
else
   error('M must be 2 or 3-d');
end
% create BBLists, and compute first set of bounds
freqactive = (1:npts);  % NO0's and 1's corresponding to which frequencies are active
freqUmax = repmat(inf,[1 npts]);    % current list of maxupper bound at each frequency
freqUmaxcube = ones(1,npts); % index of cube where FREQUMAX occurs        
freqLmax = zeros(1,npts);    % current list of maxLOWER bound at each frequency
nomcost = zeros(1,npts); 

bestL = 0;
NTIMES = opts.NTimes;
MAXCNT = opts.MaxCnt;
MBad = opts.MBadThreshold;
SABad = S*opts.ABadThreshold;
%PFlag = strcmpi(opts.Display,'On');
PFlag = 'off';
BBList = cell(npts,3); % each holds a 2-d DOUBLE, predim(50) x SpecificNumberForInfoAndScalings

% The DOUBLE represents the CUBES at a given frequency.  Some are ACTIVE, some are not.
for i=1:npts
    BBList{i,2} = 1;
%    BBList{i,3} = 0; % Parent
%    BBList{i,2} = 1; % Eliminated from considerations
end
rowone = ones(1,nblk);

i = 1;
go = 1;

while i<=npts & go==1
   BBList{i,1} = zeros(predim,cols); 
   BBList{i,3} = cell(predim,1); 
   Mcube = rcnrs(M(:,:,i),bidx,-ones(1,nblk),ones(1,nblk));
   [lowerbnd,pertsreal,pertrepreal,pertLvec,pertRvec,pertrepcomp,nomcost(i)] = ...
        wclowc(Mcube,bidx,NTIMES,MAXCNT,MBad,SABad(i));
   AllWsreal{i} = [pertsreal(:).'];
   AllWrepreal{i} = [pertrepreal(:).'];
   AllWlvec{i} = pertLvec;
   AllWrvec{i} = pertRvec;
   AllWrepcomp{i} = pertrepcomp;
   freqLmax(i) = lowerbnd;
   
   if lowerbnd>=MBad*nomcost(i) + SABad(i);
      go = 0;
      for k=i+1:npts
         AllWsreal{k} = zeros(size(pertsreal));
         AllWrepreal{k} = zeros(size(pertrepreal));
         for kk=1:length(pertLvec)
            AllWlvec{k}{kk} = zeros(size(pertLvec{kk}));
            AllWrvec{k}{kk} = zeros(size(pertRvec{kk}));
         end
         AllWrepcomp{k} = zeros(size(pertrepcomp));
      end
   elseif strcmpi(opts.LowerBoundOnly,'Off')      
       xfeas = [];
       cnt = 1;
       beta = lowerbnd;
       lbstr = [', LB: ' num2str(beta)];
       if bidx.sreal.num+bidx.repreal.num>0
          while isempty(xfeas) && cnt<=CNTMAX
             if cnt==1
                low = lowerbnd;
             else
                low = beta;
             end
             cnt = cnt + 1;
             beta = (1.25*beta+0.0001);  
             [xfeas,ubout] = usmuwcp2(M(:,:,i),blk,bidx,low,beta,ftdidx,muopt);
         end
          if isempty(xfeas)
%              if strcmpi(PFlag,'On')                                                        
%                 disp([int2str(i) '/' int2str(npts) ', NO cheap feasable' lbstr]);
%              end
             upperbnd = inf;
             dgbvec = zeros(1,DGbeta);                             
          else
%              if strcmpi(PFlag,'On')
%                 ubstr = [', UB: try' num2str(beta) ', actual ' num2str(ubout)];
%                 disp([int2str(i) '/' int2str(npts) ':' int2str(cnt) ', obtain cheap feas' lbstr ubstr]);
%              end
             upperbnd = ubout;
             dgbvec = xfeas(:)';
          end
          BBList{i,1}(1,[ALLCOR LOWER UPPER ACTIVE SCALINGS]) = ...
              [-rowone rowone lowerbnd upperbnd 1 dgbvec];
          %BBlist{i,3}{1} = {D D G G};
       else
          [upper1,di1,dbs1,gzr1,gzl1,psdami1,dgb1] = uball(M(:,:,i),bidx);
          upperbnd = upper1;
%           if strcmpi(PFlag,'On')
%              ubstr = [', UB: ' num2str(upperbnd)];
%              disp([int2str(i) '/' int2str(npts) ':' int2str(cnt) ', obtain cheap feas' lbstr ubstr]);
%           end
       end
       freqUmax(i) = upperbnd;
       i = i + 1;
   elseif strcmpi(opts.LowerBoundOnly,'On')   
      upperbnd = inf;
      dgbvec = zeros(1,DGbeta);                             
      BBList{i,1}(1,[ALLCOR LOWER UPPER ACTIVE SCALINGS]) = ...
              [-rowone rowone lowerbnd upperbnd 1 dgbvec];
      freqUmax(i) = upperbnd;
      i = i + 1;
   end
end
   

Data.AllWsreal = AllWsreal;
Data.AllWrepreal = AllWrepreal;
Data.AllWlvec = AllWlvec;
Data.AllWrvec = AllWrvec;
Data.AllWrepcomp = AllWrepcomp;
Data.BBList = BBList;
Data.freqactive = freqactive;  % unchanged since initialization
Data.freqUmax = freqUmax;
Data.freqUmaxcube = freqUmaxcube;       
Data.freqLmax = freqLmax;
Data.nomcost = nomcost;  %partially filled if loop didn't complete
Data.npts = npts;
Data.bestL = max(freqLmax);
Data.bestU = max(freqUmax);

Data.SbestL = max(freqLmax./S);
Data.SbestU = max(freqUmax./S);