source("diffEigen.R")
library(car)

myAnaCor<-function(data,ind=2:3,scal="no") {
nr<-nrow(data); nc<-ncol(data); nn<-nr*nc
mOff<-array(0,c(nr,nc,nn))
mRow<-array(0,c(nr,nr,nn))
mCol<-array(0,c(nc,nc,nn))
freq<-rep(0,nn)
k<-1
for (i in 1:nr) for (j in 1:nc) {
	mOff[i,j,k]<-1
	mRow[i,i,k]<-1
	mCol[j,j,k]<-1
	freq[k]<-data[i,j]
	k<-k+1
}
N<-sum(freq)
prop<-freq/N
rpar<-function(p) {
	res<-matrix(0,nr,nc)
	for (k in 1:nn) res<-res+p[k]*mOff[,,k]
	return(res)
	}
ppar<-function(p) {
	res<-matrix(0,nr,nr)
	for (k in 1:nn) res<-res+p[k]*mRow[,,k]
	return(res)
	}
qpar<-function(p) {
	res<-matrix(0,nc,nc)
	for (k in 1:nn) res<-res+p[k]*mCol[,,k]
	return(res)
	}
drpar<-function(p) return(mOff)
dppar<-function(p) return(mRow)
dqpar<-function(p) return(mCol)
gs<-gsvdDer(prop,ppar,qpar,rpar,dppar,dqpar,drpar,ind)
if (!(scal=="no")) gs<-gsvdScal(gs,scal=scal)
gdl<-gs$dl 
acovd<-gdl%*%(prop*t(gdl))/N
acovu<-array(0,c(2,2,nr))
gdx<-gs$dx
for (i in 1:nr) {
	gdxi<-gdx[i,,]
	gdxs<-drop(gdxi%*%prop)
	acovu[,,i]<-(gdxi%*%(prop*t(gdxi))-outer(gdxs,gdxs))/N
	}
acovv<-array(0,c(2,2,nc))
gdz<-gs$dz
for (i in 1:nc) {
	gdzi<-gdz[i,,]
	gdzs<-drop(gdzi%*%prop)
	acovv[,,i]<-(gdzi%*%(prop*t(gdzi))-outer(gdzs,gdzs))/N
	}
return(list(gs=gs,acovd=acovd,acovu=acovu,acovv=acovv,ind=ind))
}
	
myAnaPlot<-function(myAna,conf=.95){
rad<-sqrt(qchisq(conf,2))
gv<-myAna$gs$gv; gu<-myAna$gs$gu
nr<-nrow(gu); nc<-nrow(gv); ind<-myAna$ind
pdf("colPlot.pdf")
plot(gv[,ind],xlab="Dim 1",ylab="Dim 2",type="n")
for (i in 1:nc) ellipse(gv[i,ind],myAna$acovv[,,i],rad)
dev.off()
pdf("rowPlot.pdf")
plot(gu[,ind],xlab="Dim 1",ylab="Dim 2",type="n")
for (i in 1:nr) ellipse(gu[i,ind],myAna$acovu[,,i],rad)
dev.off()
}