/* This program performs a permutation test on the Procrstes distance between mean
configurations of two groups. Specimens must initially be aligned by GPA. The dataset
should consist of coordinate variables and a single character variable defining group. */

data three; set two; sps=substr(cat,1,2); drop cat; run;
proc sort data=three; by sps; run;						

proc means data=three noprint; by sps; var x1-x60; output out=means; run;

data four; set means; if _STAT_='MEAN'; drop _TYPE_ _FREQ_ _STAT_; run;

proc iml;

k=3;										* k=2,3 for 2D or 3D data;
nreps=1000;									* sets the number of permutations;

use three; read all into X; 
read all var _CHAR_ into SPS;
use four; read all into MNS;

N=NROW(X);
L=NCOL(X);
ind=design(SPS); ind=ind[+,];
p=(L/k);

MN1=MNS[1,]; MN1=shape(MN1,p,k);
MN2=MNS[2,]; MN2=shape(MN2,p,k);
I=I(p);	CaP=j(p,p,1/p);
s1=sqrt(trace((I-CaP)*MN1*MN1`*(I-CaP)));
s2=sqrt(trace((I-CaP)*MN2*MN2`*(I-CaP)));
MN1prime=(I-CaP)*MN1/s1;
MN2prime=(I-CaP)*MN2/s2;
cov2=MN1prime`*MN2prime;
CALL SVD (U,sig,V,cov2);
sig2=abs(sig); S=sig/sig2;
S=S`; S=diag(S);
H=V*S*U`;
MN2star=MN2prime*H;	

y=MN1prime-MN2star;
d=sqrt(trace(y`*y));
rho=2*arsin(d/2);

A=SPS[1,]; B=SPS[ind[1,1]+1,]; 
print 'Procrustes distance between ' A 'and   ' B ':  ' rho;

count=0;
small=min(ind);

do loop=1 to nreps;

rand=j(N,1,0); rand=normal(rand);
ord=rank(rand);
T=j(N,L); do i=1 to N; T[i,]=X[ord[i,],]; end;

grp1=T[1:ind[1,1],];
grp2=T[ind[1,1]+1:N,];

M1=grp1[+,]/NROW(grp1);
M2=grp2[+,]/NROW(grp2); 

M1=shape(M1,p,k);
M2=shape(M2,p,k);
I=I(p);	CaP=j(p,p,1/p);
s1=sqrt(trace((I-CaP)*M1*M1`*(I-CaP)));
s2=sqrt(trace((I-CaP)*M2*M2`*(I-CaP)));
M1prime=(I-CaP)*M1/s1;
M2prime=(I-CaP)*M2/s2;
cov2=M1prime`*M2prime;
CALL SVD (U,sig,V,cov2);
sig2=abs(sig); S=sig/sig2;
S=S`; S=diag(S);
H=V*S*U`;
M2star=M2prime*H;	

y=M1prime-M2star;				
d=sqrt(trace(y`*y));			
rhotemp=2*arsin(d/2);			

if rhotemp >= rho then do; count=count+1; end;

/*
if loop=1 then do; create temp from rhotemp; append from rhotemp; end;
if loop>1 then do; edit temp; append from rhotemp; end;
*/

end;

alpha=count/nreps;

print 'Number of permutations:' nreps;
print 'Probability that group means are the same:' alpha;

abort;