n% splithalft3.m
% Procedure to compute splits, analyze them and compare them.
% Input:
%		X: original data matrix
%		n,m,p: sizes of data
%	   r1,r2,r3: sizes of core
%		centopt (see tucker3.m)
%		normopt (see tucker3.m)
%     renormmode (see tucker3.m)
%		wa_rel,wb_rel,wc_rel (see tucker3.m)
%
% uses phi.m c.m 


disp('  ');
disp(' This procedure performs a SPLIT-HALF analysis on X ');
disp('  ');
disp(' NOTE: ');
disp(' 1. In SPLIT-HALF analysis, the A-mode is taken as "replication mode"');
disp('    (which means that A-mode entitites are considered a random sample');
disp('    if this does not make sense, you should rearrange your data so that the A-mode is a replication mode)');
disp(' 2. The matrix X is modified by the program, so be careful in using X again after this analysis!');
disp('  ');

disp(' The splitting into two halves can be done randomly (default), or into odd vs. even sequence numbers');
cc=input(' If you prefer odd/even split, specify "1":  ');
if isempty(cc),cc=0;end;
if cc==0
   % create random splits
  w=rand(n,1);
  [ww,wi]=sort(w);
end;
if cc==1,wi=[1:2:n 2:2:n];end;
n1=ceil(n/2);
n2=n-n1;
X1=X(wi(1:n1),:);
X2=X(wi(n1+1:n),:);

ccc=0;
disp('  ');
disp(' You will now enter the split half procedure with the options specified so far');
disp(' However, you may want to modify certain choices here (rather than rerunning the full tucker3.m)');
ccc=input(' Do you want to modify certain choices? If so, Enter "1": ');
disp('  ');
if ccc==1
  disp(' Centering:'); 
  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 ');
  fprintf(' Your centering choice is %g \n',centopt);
  cc=input(' If you want to change it, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1
     centopt=input(' Specify Centering option: ');
  end;
  disp('  ');
  disp(' Normalizing:');
  disp(' 1= within A-mode ');
  disp(' 2= within B-mode ');
  disp(' 3= within C-mode ');
  fprintf(' Your normalizing choice is %g \n',normopt);
  cc=input(' If you want to change it, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1,norm=input(' Specify Normalizing option ');end;
  disp('  ');
  if renormmode==1,disp(' Your choice was to renormalize A-mode ');end;
  if renormmode==2,disp(' Your choice was to renormalize B-mode ');end;
  if renormmode==3,disp(' Your choice was to renormalize C-mode ');end;
  if renormmode~=1 & renormmode~=2 & renormmode~=3,disp(' Your choice was not to renormalize solution');end;
  cc=input(' If you want to change it, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1
     renormmode=input(' Which mode do you want to renormalize? Enter 1 (=A), 2 (=B) or 3 (=C): ');
  end;
  disp('  ');
  fprintf(' Numbers of components for A, B and C are: %g, %g, %g \n',r1,r2,r3);
  cc=input(' If you want to change these, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1,
    disp('  ');
    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?  ');
  end;
  disp('  ');
  disp(' Additional random starts?');
  fprintf(' Your choice was to use %g additional random starts in the analysis\n',addanal);
  cc=input(' If you want to change it, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1,addanal=input(' Specify number of additional random starts ');end;
  disp('  ');
  fprintf(' Relative simplicity rotation weights for A, B and C are: %g, %g, %g \n',wa_rel,wb_rel,wc_rel);
  cc=input(' If you want to change these, specify "1" ');
  if isempty(cc),cc=0;end;
  if cc==1,
    wa_rel=input(' Specify relative weight for simplicity of A-mode  ');
    wb_rel=input(' Specify relative weight for simplicity of B-mode  ');
    wc_rel=input(' Specify relative weight for simplicity of C-mode  ');
  end;
  disp('  ');
end;

disp(' Analysis of FULL DATA');
disp('  ');
disp('  ');
% Full data analysis (reference)
% Full data Preprocessing 
cc=centopt;
if cc==1 | cc==12 | cc==13
   X=cent3(X,n,m,p,1);
   disp(' X has been centered across A-mode');
end;
if cc==2 | cc==12 | cc==23
   X=cent3(X,n,m,p,2);
   disp(' X has been centered across B-mode');
end;
if cc==3 | cc==13 | cc==23
   X=cent3(X,n,m,p,3);
   disp(' X has been centered across C-mode');
end;
cc=normopt;
disp('  ');
if cc==1 
   X=norm3(X,n,m,p,1);
   disp(' X has been normalized within A-mode');
end;
if cc==2 
   X=norm3(X,n,m,p,2);
   disp(' X has been normalized within B-mode');
end;
if cc==3
   X=norm3(X,n,m,p,3);
   disp(' X has been normalized within C-mode');
end;
% Full data analysis 
[A,B,C,H,f,iter,fp,La,Lb,Lc]=tuck3abk(X,n,m,p,r1,r2,r3,0,1e-6);
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(X,n,m,p,r1,r2,r3,1,1e-6);
   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;
[A,B,C,H]=renormsol(A,B,C,H,renormmode); 
% Full data Rotation
[AS,BT,CU,K,S,T,U,f,f1,f2a,f2b,f2c,func]=varimcoco(A,B,C,H,wa_rel,wb_rel,wc_rel);
disp(' Backnormalize A,B,C, and H');
if renormmode==1
   Ds=diag(sum(AS.^2).^.5);
   AS=AS/Ds;
   K=Ds*K;
end;
if renormmode==2
   Ds=diag(sum(BT.^2).^.5);
   BT=BT/Ds;
   K=K*kron(eye(r3),Ds);
end;
if renormmode==3
   Ds=diag(sum(CU.^2).^.5);
   CU=CU/Ds;
   K=K*kron(Ds,eye(r2));
end;
Afull=AS;
Bfull=BT;
Cfull=CU;
Kfull=K

% Split 1 data analysis (reference)
% Split 1 Preprocessing 
disp(' Analysis of SPLIT 1');
disp('  ');
disp('  ');
n_orig=n;
n=n1;
X_orig=X;
X=X1;
cc=centopt;
if cc==1 | cc==12 | cc==13
   X=cent3(X,n,m,p,1);
   disp(' X has been centered across A-mode');
end;
if cc==2 | cc==12 | cc==23
   X=cent3(X,n,m,p,2);
   disp(' X has been centered across B-mode');
end;
if cc==3 | cc==13 | cc==23
   X=cent3(X,n,m,p,3);
   disp(' X has been centered across C-mode');
end;
cc=normopt;
disp('  ');
if cc==1 
   X=norm3(X,n,m,p,1);
   disp(' X has been normalized within A-mode');
end;
if cc==2 
   X=norm3(X,n,m,p,2);
   disp(' X has been normalized within B-mode');
end;
if cc==3
   X=norm3(X,n,m,p,3);
   disp(' X has been normalized within C-mode');
end;
Xs1=X;
% Split 1 analysis 
[A,B,C,H,f,iter,fp,La,Lb,Lc]=tuck3abk(X,n,m,p,r1,r2,r3,0,1e-6);
func1=ones(1+addanal,1);
func1=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(X,n,m,p,r1,r2,r3,1,1e-6);
   func1(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(func1','6.2');disp('  ');
end;
disp(' No Postprocessing used in splits; taken care off by reference rotation!');
As1=A;
Bs1=B;
Cs1=C;
Ks1=H;

% Split 2 data analysis (reference)
% Split 2 Preprocessing 
disp(' Analysis of SPLIT 2');
n=n2;
X=X2;
cc=centopt;
if cc==1 | cc==12 | cc==13
   X=cent3(X,n,m,p,1);
   disp(' X has been centered across A-mode');
end;
if cc==2 | cc==12 | cc==23
   X=cent3(X,n,m,p,2);
   disp(' X has been centered across B-mode');
end;
if cc==3 | cc==13 | cc==23
   X=cent3(X,n,m,p,3);
   disp(' X has been centered across C-mode');
end;
cc=normopt;
disp('  ');
if cc==1 
   X=norm3(X,n,m,p,1);
   disp(' X has been normalized within A-mode');
end;
if cc==2 
   X=norm3(X,n,m,p,2);
   disp(' X has been normalized within B-mode');
end;
if cc==3
   X=norm3(X,n,m,p,3);
   disp(' X has been normalized within C-mode');
end;
Xs2=X;
% Split 2 analysis 
[A,B,C,H,f,iter,fp,La,Lb,Lc]=tuck3abk(X,n,m,p,r1,r2,r3,0,1e-6);
func2=ones(1+addanal,1);
func2=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(X,n,m,p,r1,r2,r3,1,1e-6);
   func2(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(func2','6.2');disp('  ');
end;
disp(' No Postprocessing used in splits; taken care off by reference rotation!');
As2=A;
Bs2=B;
Cs2=C;
Ks2=H;

n=n_orig;
X=X_orig;

% Compute stability indices:
disp(' Rotate towards overall solution');
disp(' ');
Ss1=As1\Afull(wi(1:n1),:);
Ss2=As2\Afull(wi(n1+1:n),:);
Ts1=Bs1\Bfull;
Ts2=Bs2\Bfull;
Us1=Cs1\Cfull;
Us2=Cs2\Cfull;

As1=As1*Ss1;
As2=As2*Ss2;
Bs1=Bs1*Ts1;
Bs2=Bs2*Ts2;
Cs1=Cs1*Us1;
Cs2=Cs2*Us2;
Ks1=Ss1\Ks1/kron(Us1',Ts1');
Ks2=Ss2\Ks2/kron(Us2',Ts2');

disp('  ');
disp(' RESULTS stability analysis');
disp(' Both splits are analyzed and rotated toards solution for full data');
disp(' One can compare complete outcomes for different splits');
disp(' or inspect correlation/congruence coefficients between corresponding columns of component matrices');
disp('  ');
disp(' Split-half solutions for B');
disp(' B''s next to each other, separated by column of 0''s ');
writescr([Bs1 zeros(m,1) Bs2],'6.2');
disp('  ');
disp(' Split-half solutions for C');
disp(' C''s next to each other, separated by column of 0''s ');
writescr([Cs1 zeros(p,1) Cs2],'6.2');
disp('  ');
disp(' Congruences for A in splits and in appropriate part of Afull');
disp('          comp.no.         split1          split2');
writescr([[1:r1]' diag(phi(Afull(wi(1:n1),:),As1)) diag(phi(Afull(wi(n1+1:n),:),As2))],'7.2');

disp(' Correlations for A in splits and in appropriate part of Afull');
disp('          comp.no.         split1          split2');
clear c;
writescr([[1:r1]' diag(phi(c(Afull(wi(1:n1),:)),c(As1))) diag(phi(c(Afull(wi(n1+1:n),:)),c(As2)))],'7.2');
disp('  ');
disp(' Congruence values for B-mode component matrix')
writescr(diag(phi(Bs1,Bs2)),'6.2')
disp(' Congruence values for C-mode component matrix');
writescr(diag(phi(Cs1,Cs2)),'6.2')

disp(' Relatively strong stability check of Core:');
disp(' (simply comparing core matrices for two splits)');
disp('  ');
disp(' Core for split1');
writescr(Ks1,'8.2');
disp(' Core for split2');
writescr(Ks2,'8.2');

disp('  ');
disp(' Weaker but Sufficient stability check of Core:');
disp(' (computing cores in two splits, using full data solutions for A,B and C) ');
disp('  ');
A1=Afull(wi(1:n1),:);
A2=Afull(wi(n1+1:n),:);
Hss1=A1\Xs1;
Hss1=Bfull\permnew(Hss1,r1,m,p);
Hss1=Cfull\permnew(Hss1,r2,p,r1);
Hss1=permnew(Hss1,r3,r1,r2);
Hss2=A2\Xs2;
Hss2=Bfull\permnew(Hss2,r1,m,p);
Hss2=Cfull\permnew(Hss2,r2,p,r1);
Hss2=permnew(Hss2,r3,r1,r2);
disp(' Core for split1');
writescr(Hss1,'8.2');
disp(' Core for split2');
writescr(Hss2,'8.2');
disp(' Core for full data');
writescr(Kfull,'8.2');

disp(' NOTE:');
disp(' For advanced analysis one can further study the results of the splits and ');
disp(' of the full data set');
disp(' Results for A, B, C, and Core of full data are given in Afull, Bfull, Cfull and Kfull');
disp(' Results for splits are given in As1,Bs1,Cs1,Ks1 and As2,Bs2,Cs2,Ks2');
disp(' Results for core in weak stability check are given in Hss1 and Hss2');
disp(' ');