% Tucker3InclBoot.m
% program to interactively carry out all steps in a tucker3 analysis
%
% Input: X = full data set as matrix with frontal slabs next to each other');
%			no missings allowed!
% during program, the matrix X is replaced by Xprep (i.e., the preprocessed version of X)
disp('  ');
disp('  ');
disp('  ');
disp(' WELCOME to the interactive TUCKER3 analysis program');
disp('  ');
disp(' (written by Henk A.L. Kiers, April 2001) ');
disp('  ');
disp('  ');
disp(' The program works fully interactively.');
disp(' Simply answer the questions, or press the Enter key for the default option.');
disp('  ');
disp(' NOTES: ');
disp(' 1. The program does not handle missing data: '); 
disp('    you should make sure that all missing values are replaced by sensible values');
disp(' 2. The data should be available in the matrix X, which is set up as follows:');
disp('    - rows correspond to A-mode entities');
disp('    - columns correspond to combinations of B- and C-mode entities,');
disp('       with B-mode entities nested within C-mode entities');
disp('   For example, a 2 x 2 x 3 array is given as follows:');
disp('     x_1_1_1  x_1_2_1  x_1_1_2  x_1_2_2  x_1_1_3  x_1_2_3');
disp('     x_2_1_1  x_2_2_1  x_2_1_2  x_2_2_2  x_2_1_3  x_2_2_3');
disp(' 3. The program may change matrix X; keep a copy of X !');
disp(' 4. Tucker3.m works in the workspace, and makes (and hence overwrites) many matrices;');
disp('     you are advised to BACKUP the contents of your workspace before analysis with Tucker3.m ');
disp(' 5. You can specify labels for A- B- and C- mode entities (in laba, labb, labc);');
disp('     otherwise standard labels will be made by the program, as laba, labb, labc.');
disp(' 6. All matrix output is tab-delimited, and can hence easily be converted to tables, after copying to Word');
disp(' 7. Faster analysis (for advanced users) can be carried out using splithalft3.m');
disp(' Note: If you are interested in doing a bootstrap analysis, it is assumed that the A-mode is the mode from which you want to resampling');
disp(' If this is not the case, first reorder your data');
disp('  ');
disp(' START of program ');
disp('  ');

conv=1e-6;
cc=0;
while cc==0
  n=input(' Specify the number of A-mode entities ');
  m=input(' Specify the number of B-mode entities ');
  p=input(' Specify the number of C-mode entities ');
  if isempty(n)|isempty(m)|isempty(p),cc=0;disp(' Please specify ALL sizes! ');
  else,cc=1;end;
end;

if exist('X')==0,disp(' Error: Data not available on X');return;end;
Xprep=X;
[rr,cc]=size(Xprep);
if rr~=n | cc~=m*p, disp(' Error: array size does not agree with input data matrix X');
   return;
end;
if exist('laba')==0,laba=num2str([1:n]');aa=length(laba(1,:));laba(:,2:aa+1)=laba;laba(:,1)='A';end
if exist('labb')==0,labb=num2str([1:m]');aa=length(labb(1,:));labb(:,2:aa+1)=labb;labb(:,1)='B';end
if exist('labc')==0,labc=num2str([1:p]');aa=length(labc(1,:));labc(:,2:aa+1)=labc;labc(:,1)='C';end
if length(laba(:,1))~=n,disp(' Error: vector of labels for A (laba) has incorrect size');return;end;
if length(labb(:,1))~=m,disp(' Error: vector of labels for B (labb) has incorrect size');return;end;
if length(labc(:,1))~=p,disp(' Error: vector of labels for C (labc) has incorrect size');return;end;



disp('  ');
cc=input(' To see anova results, specify "1":  ');
if isempty(cc);cc=0;end;
if cc==1
   disp('  ');
   threewayanova(Xprep,n,m,p);
   disp('  ');
end;

% Centering
disp('  ');
disp(' How do you want to CENTER your array ?');
disp(' 1= across A-mode ');
disp(' 2= across B-mode ');
disp(' 3= across C-mode ');
disp(' 12= across A-mode and across B-mode ');
disp(' 13= across A-mode and across C-mode ');
disp(' 23= across B-mode and across C-mode ');
cc=input(' Specify your choice:  ');
if isempty(cc);cc=0;end;
centopt=cc;
disp('  ');
if cc==1 | cc==12 | cc==13
   Xprep=cent3(Xprep,n,m,p,1);
   disp(' Data have been centered across A-mode');
end;
if cc==2 | cc==12 | cc==23
   Xprep=cent3(Xprep,n,m,p,2);
   disp(' Data have been centered across B-mode');
end;
if cc==3 | cc==13 | cc==23
   Xprep=cent3(Xprep,n,m,p,3);
   disp(' Data have been centered across C-mode');
end;

% Normalizing 
disp('  ');
disp(' How do you want to NORMALIZE your array ?');
disp(' 1= within A-mode ');
disp(' 2= within B-mode ');
disp(' 3= within C-mode ');
cc=input(' Specify your choice:  ');
if isempty(cc);cc=0;end;
normopt=cc;
disp('  ');
if cc==1 
   Xprep=norm3(Xprep,n,m,p,1);
   disp(' Data have been normalized within A-mode');
end;
if cc==2 
   Xprep=norm3(Xprep,n,m,p,2);
   disp(' Data have been normalized within B-mode');
end;
if cc==3
   Xprep=norm3(Xprep,n,m,p,3);
   disp(' Data have been normalized within C-mode');
end;

disp('  ');
disp(' NOTE: ');
disp('   The preprocessed data are now available in Xprep');
disp('   Xprep can be used for analyses outside the Tucker3.m program');
disp('  ');

disp('  ');
cc=input(' To see PCA of MEAN, specify "1":  ');
if isempty(cc);cc=0;end;
while cc==1
   [Y,ev,A1,B1,C1,A2,B2,C2]=pcamean(Xprep,n,m,p);
   disp('  ');
   disp(' To study/plot the output you can now temporarily leave the interactive program');
   disp('   to get back, type "return"');
   ccc=input(' If you want to temporarily leave the program, specify "1": '); 
   if isempty(ccc);ccc=0;end;
   if ccc==1,
      disp('  ');
      disp(' You leave the interactive program and can manipulate/study current matrices. ');
      disp(' You may use writescr.m, or use any scatter-plotting routine you like')
      disp(' (e.g., scatter, etc.,or "plotje.m").');
      disp(' If you want to make plots, you may use labels in laba, labb, labc.');
      disp('  ');
      disp(' To return to interactive program, type "return"');
      keyboard;
   end;
   disp('  ');
   ccc=input(' If you want another PCA of mean analysis, specify "1":  ');
   if isempty(ccc);ccc=0;end;
   if ccc~=1,cc=2;end;
end;

disp(' ');
disp(' You are now about to do a Tucker3 analysis.');
disp(' You have to specify the dimensionalities to use.');
disp(' To search these, you are advised to first run PCASUP analyses.');
disp(' (PCASUP analyses are PCA''s of supermatrices with slices of the 3way array next to each other;')
disp('  thus 3 supermatrices are analyed by PCA, and for each we get a component matrix, and eigenvalues;');
disp('  the eigenvalues may give an indication as to the required number of components for each mode.)');
disp('  The results can next be used to find useful dimensionalities by the generalized scree test ');
disp('  (see Timmerman & Kiers, 2000, Kiers & der Kinderen, 2003)');
disp(' ');
cc=input(' If you want to do the PCASUPs for Choosing your dimensionality , specify "1":   ');
if isempty(cc);cc=0;end;
disp('  ');
if cc==1
   disp(' For the generalized scree test it is needed to indicate the maximum number of dimensions for each mode you want to study');
   maxa=input(' Up to how many A-mode components do you want to use?  ');
   maxb=input(' Up to how many B-mode components do you want to use?  ');
   maxc=input(' Up to how many C-mode components do you want to use?  ');
   out=tuckrunsApproxFit(Xprep,n,m,p,maxa,maxb,maxc);
   disp('  ');
   disp('           numbers of components  fit    total number');
   disp('            A       B       C     perc.  of components');
   disp('  ');
   writescr(out,'6.2');
   disp('  ');
   disp(' This info remains available in matrix out ');
   disp('  ');
   disp(' Suggestion based on Convex Hull procedure (Ceulemans & Kiers, 2006) ' );
   disp(' -- programmed by Urbano Lorenzo Seva, 2007 --' );
   DimSelector(out);
   disp(' ');
   disp(' ');
   disp(' Figure 1 now gives fit values of all solutions plotted against P+Q+R');
   disp(' Figure 2 gives fit values of all solutions plotted against number of free parameters');
   tucker3dimensionalityplot;
   disp(' ');
end;


ccc=input(' If you want to temporarily leave the program, specify "1": '); 
if isempty(ccc);ccc=0;end;
if ccc==1,
   disp('  ');
   disp(' You leave the interactive program and can manipulate/study current matrices ');
   disp(' To return to interactive program, type "return"');
   keyboard;
end;

disp('  ');
disp(' You can now do a TUCKER3 ANALYSIS');
disp('  ');
cc=0;
while cc==0
  r1=input(' How many A-mode components do you want to use?  ');
  r2=input(' How many B-mode components do you want to use?  ');
  r3=input(' How many C-mode components do you want to use?  ');
  if isempty(r1)|isempty(r2)|isempty(r3),cc=0;disp(' Please specify ALL numbers of components! ');
  else,cc=1,end;
end;
disp(' ');
disp(' By default, only a rationally started analysis run will be carried out.');
disp(' To decrease the chance of missing the optimal solution, you may use');
disp(' additional, randomly started runs.');
addanal=input(' If you want additional runs, specify how many (e.g., 4):  ');
if isempty(addanal),addanal=0;end;
disp('  ');
[A,B,C,H,f,iter,fp,La,Lb,Lc]=tuck3abk(Xprep,n,m,p,r1,r2,r3,0,conv);
func=ones(1+addanal,1);
func=fp;
for run=1:addanal
   disp('  ');
   fprintf(' Run no. %g \n',run+1);
   disp('  ');
   [Ar,Br,Cr,Hr,fr,iterr,fpr,Lar,Lbr,Lcr]=tuck3abk(Xprep,n,m,p,r1,r2,r3,1,conv);
   func(run+1)=fpr;
   if fpr>1.0001*fp % if fit more than .01% better is found, replace solution
      A=Ar;B=Br;C=Cr;H=Hr;f=fr;iter=iterr;fp=fpr;La=Lar;Lb=Lbr;Lc=Lcr;
   end;
   disp('  ');
end;
disp('  ');
if addanal>=1
   disp('  ');
   disp(' Fit % values from all runs:');
   disp(' ');
   writescr(func','6.2');
   disp('  ');
end;
fprintf(' Tucker3 analysis with %gx%gx%g components gave a fit of %4.2f %% \n',r1,r2,r3,fp); 
disp('  ');
disp('  ');
disp('  ');

disp('  ');
ccc=input(' If you want to temporarily leave the program, specify "1": '); 
if isempty(ccc);ccc=0;end;
if ccc==1,
   disp('  ');
   disp(' You leave the interactive program and can manipulate/study current matrices ');
   disp(' To return to interactive program, type "return"');
   keyboard;
end;
disp('  ');


% Prepare plotting matrices
disp('  ');
disp(' PLOTTING ');
disp('  ');
disp(' No special plotting routines are incorporated, but you can use');
disp(' any scatter-plotting routine you like (e.g., scatter, etc.,or "plotje.m")');
disp('  ');
disp(' If you want to make plots of entities of one mode, you can use:');
disp(' matrices Aplot, Bplot, or Cplot, and use labels in laba, labb, labc.');
Aplot=A*diag(diag(La).^.5);
Bplot=B*diag(diag(Lb).^.5);
Cplot=C*diag(diag(Lc).^.5);
disp('  ');
disp(' If you want to make plots of entities of two modes, ');
disp(' with the third projected in it as axes, you can use the coordinates in:');
disp(' BAplot with those in C,or ACplot with those in B, or CBplot with those in A');
disp(' You will still have to make proper labels for these combined mode plots ');
Aplot=A*diag(diag(La).^.5);
Bplot=B*diag(diag(Lb).^.5);
Cplot=C*diag(diag(Lc).^.5);
Ha=H;Hb=permnew(Ha,r1,r2,r3);Hc=permnew(Hb,r2,r3,r1);
CBplot=kron(C,B)*Ha';
ACplot=kron(A,C)*Hb';
BAplot=kron(B,A)*Hc';
disp('  ');

disp('  ');
ccc=input(' If you want to temporarily leave the program, specify "1": '); 
if isempty(ccc);ccc=0;end;
if ccc==1,
   disp('  ');
   disp(' You leave the interactive program and can manipulate/study current matrices ');
   disp(' To return to interactive program, type "return"');
   keyboard;
end;
disp('  ');

disp(' It is sometimes useful to RENORMALIZE components, e.g., such that')
disp(' the influence of the different A-mode components on the core is equalized');
disp(' This will usually lead to correlated components in that mode, and ');
disp(' to a more simple structure for the components in the OTHER modes, ');
disp(' ');