function [Y,ev,A1,B1,C1,A2,B2,C2]=pcamean(X,n,m,p);
%[Y,ev,A1,B1,C1,A2,B2,C2]=pcamean(X,n,m,p);
% Program performs Pca of mean, aggregated over mode number indicated by aggregmode
%
% Input:
%	 X: data 
%	 n,m,p: array size
% Output:
%	ev = eigenvalues
%  for r components:
%    A1,B1,C1 solution based on varimax rotation of loadings 
%    A2,B2,C2 solution based on oblique Harris-Kaiser Independent Cluster rotation of loadings 
% 
% program uses normvari.m c.m nrm.m


% 	 aggregmode = mode over which means are computed
aggregmode=input(' Over which mode do you want to compute means? Enter 1,2 or 3 (for A,B, or C)?');
ev=[];
A1=[];B1=[];C1=[];
A2=[];B2=[];C2=[];
if aggregmode==1
   fprintf(' A %g by %g matrix for modes B x C will be computed \n',m,p);
   loadingsmode=input(' Which mode will serve as "variables mode"? Enter 2 (=B) or 3 (=C): ');
   % loadingsmode = mode which is to be considered as the one 
   %					for which you want to obtain simple structured loadings
   %	Note that the other mode then is the "scores" mode 
   %			with scores normalized to unit sums of squares
   Y=zeros(m,p);	% B-mode x C-mode
   disp('  ');
   if loadingsmode==2, load='P';disp(' B will contain loadings; C will contain scores');end;
   if loadingsmode~=2, load='Q';disp(' C will contain loadings; B will contain scores');end;
   disp('  ');
end;
if aggregmode==2
   fprintf(' A %g by %g matrix for modes C x A will be computed \n',p,n);
   loadingsmode=input(' Which mode will serve as "variables mode"? Enter 1 (=A) or 3 (=C):  ');
   X=permnew(X,n,m,p);	% m x p x n
   Y=zeros(p,n);		% C-mode x A-mode
   disp('  ');
   if loadingsmode~=1, load='P';disp(' C will contain loadings; A will contain scores');end;
   if loadingsmode==1, load='Q';disp(' A will contain loadings; C will contain scores');end;
   disp('  ');
end;
if aggregmode==3
   fprintf(' A %g by %g matrix for modes A x B will be computed \n',n,m);
   loadingsmode=input(' Which mode will serve as "variables mode"? Enter 1 (=A) or 2 (=B): ');
   X=permnew(X,n,m,p);	% m x p x n
   X=permnew(X,m,p,n);	% p x n x m
   Y=zeros(n,m);		% A-mode x B-mode
   disp('  ');
   if loadingsmode~=2, load='P';disp(' A will contain loadings; B will contain scores');end;
   if loadingsmode==2, load='Q';disp(' B will contain loadings; A will contain scores');end;
   disp('  ');
end;

mm=mean(X)';
Y(:)=mm;  % P-mode x Q-mode
[nr,nc]=size(Y);
if ssq(Y)<1e-30
   disp(' Means to analyze are exteremely small; probably data have been centered across the ');
   disp(' same mode as over which means are computed, which makes no sense.');
   disp(' Therefore, PCA of means is not carried out !');
   return;
end;
disp(' ');
cc=input(' If you want to standardize "variables", type "1":  ');
if isempty(cc),cc=0;end;
if cc==1
   if load=='Q',Y=nrm(c(Y));end;
   if load=='P',Y=nrm(c(Y'))';end;
   disp(' "Variables" have been standardized ');
end;
disp('  ');
[P,D,Q]=svd(Y,0); % P for rows  Q for columns
ev=diag(D).^2;
[rr,rrr]=size(ev);
disp(' Eigenvalues and cumulative percentages (Scree plot in Figure)');
writescr([ev cumsum(ev)/sum(ev)*100],'10.2');
plot([1:rr]',ev);
r=input(' How many components do you want to maintain for rotation? ');
disp(' ');
if load=='P'
  % varimax of P-mode loadings
  [P1,T]=normvari(P(:,1:r)*D(1:r,1:r));
  Q1=Q(:,1:r)*T;
  sg=diag(sign(sum(P1))); % reflect if necessary
  P1=P1*sg;
  Q1=Q1*sg;
  % Harris-Kaiser of P-mode loadings
  [P2,T]=normvari(P(:,1:r));
  Q2=Q(:,1:r)*D(1:r,1:r)*T;
  Ds=diag(sum(Q2.^2));
  P2=P2*Ds^.5;
  Q2=Q2/Ds^.5;
  sg=diag(sign(sum(P2))); % reflect if necessary
  P2=P2*sg;
  Q2=Q2*sg;
end;     
if load=='Q'
  % varimax of Q-mode loadings
  [Q1,T]=normvari(Q(:,1:r)*D(1:r,1:r));
  P1=P(:,1:r)*T;
  sg=diag(sign(sum(Q1))); % reflect if necessary
  P1=P1*sg;
  Q1=Q1*sg;
  % Harris-Kaiser of Q-mode loadings
  [Q2,T]=normvari(Q(:,1:r));
  P2=P(:,1:r)*D(1:r,1:r)*T;
  Ds=diag(sum(P2.^2));
  Q2=Q2*Ds^.5;
  P2=P2/Ds^.5;
  sg=diag(sign(sum(Q2))); % reflect if necessary
  P2=P2*sg;
  Q2=Q2*sg;
end;     

if aggregmode==1 
   B1=P1;B2=P2;C1=Q1;C2=Q2;
end;
if aggregmode==2
   C1=P1;C2=P2;A1=Q1;A2=Q2;
end;
if aggregmode==3
   A1=P1;A2=P2;B1=Q1;B2=Q2;
end;
disp(' A1, B1 and/or C1 contain solution based on varimax rotation of loadings');
disp(' A2, B2 and/or C2 contain solution based on oblique "HKIC" orthomax rotation of loadings');
disp(' ');