function [AS,BT,CU,K,S,T,U,f,f1,f2a,f2b,f2c,func]=varimcoco(A,B,C,H,wa_rel,wb_rel,wc_rel,rot1,rot2,rot3,nanal);
% [AS,BT,CU,K,S,T,U,f,f1,f2a,f2b,f2c,func]=varimcoco(A,B,C,H,wa_rel,wb_rel,wc_rel,rot1,rot2,rot3,nanal);
%
% program for varimax rotation of core and component matrix rotation to simple structure ORTHCOCO.M
% 
% see article Kiers, H.A.L. (1998), "Joint orthomax rotation of the core and component 
% matrices resulting from three-mode principal components analysis", Journal of Classification, 245-263. 
% 
%
% input: H = core array (r1 x r2r3, with second mode nested within third)
%        A, B, C = columnwise orthonormal component matrices
%			wa_rel,wb_rel,wc_rel: relative weights for component matrices
%			rot1,rot2,rot3: binary indicators whether (1) or not (0) to rotate this mode
%				(when omitted, all modes are rotated)
%			nanal=number of random starts (when omitted: 5)
%
% output: AS,BT,CU: rotated component matrices
% 			K: rotated core
% 			S, T, U: rotation matrices
%
% Notational differences with article: 
%		n,m,p 		are used for 	I,J,K in the paper
%		r1, r2, r3 	are used for 	P, Q, R in the paper
% 		H is used 	instead of 	G in the paper
% 
%
% uses orthmax2.m, permnew.m, ssq.m

conv=1e-6;					  % CONVERGENCE CRITERION
if nargin<10,nanal=5;end; % default 5 random starts
if nargin<8,rot1=1;rot2=1;rot3=1;end;
[n,r1]=size(A);
[m,r2]=size(B);
[p,r3]=size(C);
disp(' ');
fprintf('Results from %g random starts \n',nanal);
disp('  ');
% CHECK INPUT SPECIFICATION
sz=size(H);
if sz(1)~=r1 | sz(2)~=r2*r3,disp(' Sizes of A,B and C incompatible with H');beep;return;end;

gam1=1;gam2=1;gam3=1;
gama=1;gamb=1;gamc=1;

% Compute Natural Weights for core rotation: 
% ensure natural proportion for three core terms
w1=1/r2/r3;
w2=1/r1/r3;
w3=1/r1/r2;
ww=(w1+w2+w3);	 
w1=w1/ww;	% normalizes weights to sum of 1
w2=w2/ww;
w3=w3/ww;
ssh=ssq(H)^2/r1/r2/r3;	% mean sq of squared core elements 
			% makes contribution of qmax on uniform H equal to 1
w1=w1/ssh;	
w2=w2/ssh;
w3=w3/ssh;
v1=gam1*w1/(w1+w2+w3);
v2=gam2*w2/(w1+w2+w3);
v3=gam3*w3/(w1+w2+w3);

% Compute WEIGHTS for columnwise orthonormal component matrices;
% relative weights times natural weights (these are "natural" given columnwise orthonormality!)
wa=n/r1;
wb=m/r2;
wc=p/r3;
wa=wa_rel*wa;
wb=wb_rel*wb;
wc=wc_rel*wc;


% START LOOP FOR .. RUNS OF ALGORITHM
Ssup=[];Tsup=[];Usup=[];func=[];
for ana=1:nanal
  rand('seed',ana^2);
  if rot1==1,S=orth(rand(r1,r1)-.5);end;
  if rot2==1,T=orth(rand(r2,r2)-.5);end;
  if rot3==1,U=orth(rand(r3,r3)-.5);end;

% Initialize Y, AS, BT and CU
% multiply by weights^.25, see formula (7)

  Y=(w1+w2+w3)^(.25)*S*H*kron(U',T');
  AS=wa^(.25)*A*S';
  BT=wb^(.25)*B*T';
  CU=wc^(.25)*C*U';

% EVALUATE FUNCTION VALUE AT START

	% three-way orthomax part
  f1=ssq(Y'.*Y')-v1/(r2*r3)*sum((sum(Y'.*Y')).^2);  
  Y=permnew(Y,r1,r2,r3);	% Y of order r2 x r3r1
  f1=f1-v2/(r1*r3)*sum((sum(Y'.*Y')).^2);  
  Y=permnew(Y,r2,r3,r1); % Y of order r3 x r1r2
  f1=f1-v3/(r1*r2)*sum((sum(Y'.*Y')).^2);
  Y=permnew(Y,r3,r1,r2); % Y of order r1 x r2r3
	% orthomax for component matrices
  L=AS;L=L.*L;f2a=ssq(L)-gama/n*sum((sum(L)).^2);
  L=BT;L=L.*L;f2b=ssq(L)-gamb/m*sum((sum(L)).^2);
  L=CU;L=L.*L;f2c=ssq(L)-gamc/p*sum((sum(L)).^2);
  f2=f2a+f2b+f2c;
  f=f1+f2;
%  fprintf(' At start,  f = % 8.6f ',f);
%  fprintf('   (core: %8.6f; ',f1);
%  fprintf(' A %6.3f B %6.3f C %6.3f\n',f2a,f2b,f2c);
 
  % START OF ITERATIVE PART
  iter=0;
  fold=f-2*conv*abs(f);
  while abs(f-fold)>conv*abs(f)  
    iter=iter+1;
    fold=f;

    % UPDATING S
    if rot1==1
        % orthomax of     Y' = (w1+w2+w3)^.25 H_F ' with gamma = v1
	% + orthomax of   AS = wa^.25A 		    with gamma = gama
       
        [K,AS,S1]=orthmax2(Y',AS,v1,gama,conv);
        S=S1'*S;
    end;
    if rot1==0,K=Y';end;
    Y=permnew(K(1:r2*r3,:)',r1,r2,r3);		% Y of order r2 x r3r1

    % UPDATING T
    if rot2==1
        % orthom of   Y' = (w1+w2+w3)^.25 H_F ' with gamma=v2
	% + orthom of BT = wb^.25B              with gamma=gamb
       
        [K,BT,T1]=orthmax2(Y',BT,v2,gamb,conv);
        T=T1'*T;
    end;
    if rot2==0,K=Y';end;
    Y=permnew(K(1:r1*r3,:)',r2,r3,r1); 		% Y of order r3 x r1r2

    % UPDATING U
    if rot3==1
        % orthom of   Y' = (w1+w2+w3)^.25 H_F ' with gamma=v3
	% + orthom of CU = wc^.25C              with gamma=gamc
        [K,CU,U1]=orthmax2(Y',CU,v3,gamc,conv);
        U=U1'*U;
    end;
    if rot3==0,K=Y';end;
    Y=permnew(K(1:r1*r2,:)',r3,r1,r2); 		% Y of order r1 x r2r3

    % EVALUATE ORTHCOCO FUNCTION:
	% three-way orthomax part
    f1=ssq(Y'.*Y')-v1/(r2*r3)*sum((sum(Y'.*Y')).^2);  
    Y=permnew(Y,r1,r2,r3);	% Y of order r2 x r3r1
    f1=f1-v2/(r1*r3)*sum((sum(Y'.*Y')).^2);  
    Y=permnew(Y,r2,r3,r1); % Y of order r3 x r1r2
    f1=f1-v3/(r1*r2)*sum((sum(Y'.*Y')).^2);
    Y=permnew(Y,r3,r1,r2); % Y of order r1 x r2r3
	% orthomax for component matrices
    L=AS;L=L.*L;f2a=ssq(L)-gama/n*sum((sum(L)).^2);
    L=BT;L=L.*L;f2b=ssq(L)-gamb/m*sum((sum(L)).^2);
    L=CU;L=L.*L;f2c=ssq(L)-gamc/p*sum((sum(L)).^2);
    f2=f2a+f2b+f2c; 
    f=f1+f2;
%   fprintf(' At iter %2i,  f = % 8.6f ',iter,f);
%   fprintf('   (core: %8.6f; A,B,C: %6.3f)',f1,f2);
%   fprintf(' %6.3f %6.3f %6.3f\n',f2a,f2b,f2c);
  end;
  fprintf(' Run no. %g f= %6.3f ',ana,f);
  fprintf(' (core: %6.3f;',f1);
  fprintf(' A: %6.3f; B: %6.3f; C: %6.3f), %2i iters\n',f2a,f2b,f2c,iter);
  func(ana)=f;Ssup=[Ssup S(:)];Tsup=[Tsup T(:)];Usup=[Usup U(:)];
end; 	

% FINISH LOOP FOR ANALYSES

% REPORT BEST SOLUTION
[f,maxi]=max(func);
%disp(' Rotation matrices:');
S(:)=Ssup(:,maxi);T(:)=Tsup(:,maxi);U(:)=Usup(:,maxi);
disp('  ');
fprintf(' Three-way orthomax function value for best solution is %12.4f \n',f);
%disp(' Rotated core:');
AS=A*S';BT=B*T';CU=C*U';
% reflect to positive sums in A,B and C
sg=sign(sum(AS));sg(sg==0)=1;S=diag(sg)*S;AS=AS*diag(sg);
sg=sign(sum(BT));sg(sg==0)=1;T=diag(sg)*T;BT=BT*diag(sg);
sg=sign(sum(CU));sg(sg==0)=1;U=diag(sg)*U;CU=CU*diag(sg);
% compute core
K=S*H*kron(U',T');

% recompute final orthomax values
L=AS;L=L.*L;f2a=(ssq(L)-gama/n*sum((sum(L)).^2))*n/r1;
L=BT;L=L.*L;f2b=(ssq(L)-gamb/m*sum((sum(L)).^2))*m/r2;
L=CU;L=L.*L;f2c=(ssq(L)-gamc/p*sum((sum(L)).^2))*p/r3;
% three-way orthomax part
Y=((w1+w2+w3)^.25)*K;
f1=ssq(Y'.*Y')-v1/(r2*r3)*sum((sum(Y'.*Y')).^2);  
Y=permnew(Y,r1,r2,r3);	% Y of order r2 x r3r1
f1=f1-v2/(r1*r3)*sum((sum(Y'.*Y')).^2);  
Y=permnew(Y,r2,r3,r1); % Y of order r3 x r1r2
f1=f1-v3/(r1*r2)*sum((sum(Y'.*Y')).^2);
Y=permnew(Y,r3,r1,r2); % Y of order r1 x r2r3
disp(' ');
disp(' Varimax values of core and AS, BT and CU (all based on "natural weights")');
disp(' ')
disp('                Core             A               B               C ');
writescr([f1 f2a f2b f2c],'12.3');
disp(' ');
disp(' These simplicity values are based on "natural" weights and therefore comparable across matrices'); 
disp(' When multiplied by the relative weights, they give the contribution to the overall simplicity value'); 
disp(' (they are I^2/p, J^2/q or K^2/r, resp., times the sum of the variances of squared values)');