% bootstraptucker3.m
%
% Note:  make sure that A-mode is resampling mode!
% 
% input: X = data matrix
%        Xprep  preprocessed data
%        A,B,C,G  = (possibly rotated) solution
%		 n,m,p =  order of array
%        r1,r2,r3 = numbers of components
%		 conv = convergence criterion
%        centopt = centering options
%        normopt = centering options
%
% Uses currently available sample solution, and bootstrap 
% via orthogonal rotation towards full solution
% 
% Note: preprocessing must be done in same way as for sample analysis
%		resampling: within A
%		start bootstraps: two starts: rational start, and solution from the sample
%
% uses norm3.m cent3.m procr.m tuck3abkrep.m

ccc=0;
while ccc==0
  n_boot=input(' How many bootstrap samples do you want to draw? (1000 recommended) ');
  if isempty(n_boot);ccc=0;else ccc=1;end;
end;

% setting up supermatrices for collecting bootstrap matrices

BB=zeros(m*r2,n_boot);
CC=zeros(p*r3,n_boot);
GG=zeros(r1*r2*r3,n_boot);
fpfp=zeros(1,n);
BBB=zeros(m*r2,n);
CCC=zeros(p*r3,n);
GGG=zeros(r1*r2*r3,n);
fpfpfp=zeros(1,n);
for i=1:n_boot     %+n
  if i<=n_boot
      fprintf('bootstrap run %g \n',i);
      b=ceil(rand(n,1)*n)';  % bootstrap sample
      Xb=X(b,:);
      n_bootsample=n;
      Astart=A;
  end;
  if i>n_boot
      fprintf('expanded data set run, for BCa %g \n',i);
      Xb=[X;X(i-n_boot,:)];
      n_bootsample=n+1;
      Astart=[A;A(i-n_boot,:)];
  end;
  
  % preprocessing
  Xprepb=Xb;
  cc=centopt;
  if cc==1 | cc==12 | cc==13
    Xprepb=cent3(Xprepb,n_bootsample,m,p,1);
  end;
  if cc==2 | cc==12 | cc==23
   Xprepb=cent3(Xprepb,n_bootsample,m,p,2);
  end;
  if cc==3 | cc==13 | cc==23
   Xprepb=cent3(Xprepb,n_bootsample,m,p,3);
  end;
  cc=normopt;
  if cc==1 
   Xprepb=norm3(Xprepb,n_bootsample,m,p,1);
  end;
  if cc==2 
   Xprepb=norm3(Xprepb,n_bootsample,m,p,2);
  end;
  if cc==3
   Xprepb=norm3(Xprepb,n_bootsample,m,p,3);
  end;
  
  % Tucker3, use two runs: rational start and sample solution start 
  start=2;  % start from sample solution 
  [Ab,Bb,Cb,Gb,fb,iter,fpb,Lab,Lbb,Lcb]=tuck3abkrep(Xprepb,n_bootsample,m,p,r1,r2,r3,start,conv,Astart,B,C,G);
  start=0;  % rational start
  [Abb,Bbb,Cbb,Gbb,fbb,iter,fpbb,Lab,Lbb,Lcb]=tuck3abkrep(Xprepb,n_bootsample,m,p,r1,r2,r3,start,conv,Astart,B,C,G);
  if fbb<fb,fb=fbb;Ab=Abb;Bb=Bbb;Cb=Cbb;Gb=Gbb;fpb=fpbb;end;

  if optimalmatch==0
  %matching via orthogonal rotation towards full solution
  % minimize || Bb T - B ||, || Cb U - C ||, || S' Gb - G ||.
  % using svd B'Bb = PDQ', then T=QP', etc.
  [P,D,Q]=svd(B'*Bb);
  T=Q*P';
  [P,D,Q]=svd(C'*Cb);
  U=Q*P';
  Bb=Bb*T;
  Cb=Cb*U;
  Gb=Gb*kron(U,T);
  [P,D,Q]=svd(G*Gb');
  S=Q*P';
  Gb=S'*Gb;
  end;

  if optimalmatch==1
  %matching via optimal transformation towards full solution
  % minimize || Bb T - B ||, || Cb U - C ||, || S Gb - G ||.
  % using linear regressions
  T=Bb\B;
  U=Cb\C;
  Bb=Bb*T;
  Cb=Cb*U;
  Gb=Gb*kron(inv(U'),inv(T'));
  S=G/Gb;
  Gb=S*Gb;
 % Ab=Ab/S;
 % fprintf('fit check %8.6f \n',abs(ssq(Ab*Gb*kron(Cb,Bb)'-Xprepb)-fb));
  end;


  % collect bootstrap solutions
  if i<=n_boot
    BB(:,i)=Bb(:);
    CC(:,i)=Cb(:);
    GG(:,i)=Gb(:);
    fpfp(i)=fpb;
  end;
  % collect expanded data solutions
  if i>n_boot
    BBB(:,i-n_boot)=Bb(:);
    CCC(:,i-n_boot)=Cb(:);
    GGG(:,i-n_boot)=Gb(:);
    fpfpfp(i-n_boot)=fpb;
  end;

end;

se_B=zeros(m,r2);
se_C=zeros(p,r3);
se_G=zeros(r1,r2*r3);
se_B(:)=std(BB',1);
se_C(:)=std(CC',1);
se_G(:)=std(GG',1);
se_fp=std(fpfp',1);

m_B=zeros(m,r2);
m_C=zeros(p,r3);
m_G=zeros(r1,r2*r3);
m_B(:)=mean(BB',1);
m_C(:)=mean(CC',1);
m_G(:)=mean(GG',1);
m_fp=mean(fpfp',1);

lo_B=zeros(m,r2);
lo_C=zeros(p,r3);
lo_G=zeros(r1,r2*r3);
up_B=zeros(m,r2);
up_C=zeros(p,r3);
up_G=zeros(r1,r2*r3);

Bint=[];Cint=[];Gint=[];
[lo,up]=percentile95(BB');
lo_B(:)=lo;
up_B(:)=up;
Bint(:,1:2:2*r2)=lo_B;
Bint(:,2:2:2*r2)=up_B;
[lo,up]=percentile95(CC');
lo_C(:)=lo;
up_C(:)=up;
Cint(:,1:2:2*r3)=lo_C;
Cint(:,2:2:2*r3)=up_C;
[lo,up]=percentile95(GG');
lo_G(:)=lo;
up_G(:)=up;
Gint(:,1:2:2*r2*r3)=lo_G;
Gint(:,2:2:2*r2*r3)=up_G;
[lo,up]=percentile95(fpfp');
lo_fp=lo;
up_fp=up;