#
#   aspect package
#   Copyright (C) 2005  Jan de Leeuw <deleeuw@stat.ucla.edu>
#   UCLA Department of Statistics, Box 951554, Los Angeles, CA 90095-1554
#   
#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#   
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
###################################################################
#
# version 0.1.00, 2006-04-15,		first official release
#
# todo
# -- monotone regression
# -- nonmonotone criteria
# -- other aspects
# -- use splines, getting rid of expandFrame
#

maxAspect<-function(data,aspect=aspectSum,itmax=100,eps=1e-6,extra=1) {
m<-dim(data)[2]; n<-dim(data)[1]; r<-diag(m); fold<--Inf; itel<-1
ncat<-sapply(1:m,function(j) length(table(data[,j])))
ccat<-c(0,cumsum(ncat)); y<-list(); for (j in 1:m) y<-c(y,list(1:ncat[j]))
burt<-burtTable(data); d<-diag(burt)
for (j in 1:m) {
    indj<-(ccat[j]+1):ccat[j+1]
	dj<-d[indj]
	y[[j]]<-y[[j]]-sum(dj*y[[j]])/n
	y[[j]]<-y[[j]]/sqrt(sum(dj*y[[j]]*y[[j]])/n)
	}
repeat {
	for (j in 1:m) {
		indj<-(ccat[j]+1):ccat[j+1]
		for (l in 1:m) {
			indl<-(ccat[l]+1):ccat[l+1]
			r[j,l]<-sum(y[[j]]*(burt[indj,indl]%*%y[[l]]))/n
			}
		}
	a<-aspect(r,extra); f<-a$f; g<-a$g
	print(f)
	for (j in 1:m) {
		indj<-(ccat[j]+1):ccat[j+1]; y[[j]]<-rep(0,ncat[j]); dj<-d[indj]
		for (l in 1:m) {
			indl<-(ccat[l]+1):ccat[l+1]
			if (j != l) y[[j]]<-y[[j]]+g[j,l]*burt[indj,indl]%*%y[[l]]
			}
		y[[j]]<-y[[j]]/dj
		y[[j]]<-y[[j]]-sum(dj*y[[j]])/n
		y[[j]]<-y[[j]]/sqrt(sum(dj*y[[j]]*y[[j]])/n)
		}
	if (((f-fold) < eps) || (itel == itmax)) break
	itel<-itel+1; fold<-f
	}
return(list(f=f,y=y,r=r,e=eigen(r,only.values=TRUE)$values))
}


aspectSum<-function(r,pow) {
m<-dim(r)[1]
list(f=sum(r^pow),g=pow*r^(pow-1))
}

aspectAbs<-function(r,pow) {
list(f=sum(abs(r)^pow)/2,g=pow*sign(r)*(abs(r)^(pow-1)))
}

aspectSMC<-function(r,nvar) {
m<-dim(r)[1]; ind<-which(1:m!=nvar)
rn<-as.vector(r[nvar,ind]); rr<-r[ind,ind]
b<-solve(rr,rn)
g<-diag(m); g[nvar,ind]<--b; g[ind,nvar]<--b
g[ind,ind]<-outer(b,b)
list(f=sum(b*rn),g=-g)
}

aspectSumSMC<-function(r,extra) {
m<-dim(r)[1]; f<-0; g<-matrix(0,m,m)
for (j in 1:m) {
	a<-aspectSMC(r,j)
	f<-f+(a$f); g<-g+(a$g)
	}
list(f=f,g=g)
}

aspectEigen<-function(r,p) {
s<-eigen(r); v<-s$vectors[,1:p]
list(f=sum(s$values[1:p]),g=v%*%t(v))
}

aspectDeterminant<-function(r,extra) {
list(f=-log(det(r)),g=-solve(r))
}

lineals<-function(data,itmax=100,eps=1e-6) {
m<-dim(data)[2]; n<-dim(data)[1]; r<-diag(m); t<-diag(m); fold<-Inf; itel<-1
ncat<-sapply(1:m,function(j) length(table(data[,j])))
ccat<-c(0,cumsum(ncat)); y<-list(); for (j in 1:m) y<-c(y,list(1:ncat[j]))
burt<-burtTable(data); d<-diag(burt)
for (j in 1:m) {
	indj<-(ccat[j]+1):ccat[j+1]; dj<-d[indj]
	y[[j]]<-y[[j]]-sum(dj*y[[j]])/n
	y[[j]]<-y[[j]]/sqrt(sum(dj*y[[j]]*y[[j]]))
	}
repeat {
	f<-0
	for (j in 1:m) {
		indj<-(ccat[j]+1):ccat[j+1]; yj<-y[[j]]
		for (l in 1:m) {
			indl<-(ccat[l]+1):ccat[l+1]; dl<-d[indl]; yl<-y[[l]]
			r[j,l]<-sum(burt[indj,indl]*outer(yj,yl))
			c<-burt[indj,indl]%*%diag(1/pmax(1,dl))%*%burt[indl,indj]
			t[j,l]<-sum(c*outer(yj,yj))			
			f<-f+(t[j,l]-r[j,l]^2)
			}
		}
	print(f)
	for (j in 1:m) {
		indj<-(ccat[j]+1):ccat[j+1]; nc<-ncat[j]; yj<-y[[j]]
		dj<-d[indj]; c<-matrix(0,nc,nc)
		for (l in 1:m) {
			if (j != l) {
				indl<-(ccat[l]+1):ccat[l+1]; dl<-d[indl]; yl<-y[[l]]
				u<-burt[indj,indl]%*%(diag(1/pmax(1,dl))-2*outer(yl,yl))%*%burt[indl,indj] 
				c<-c+u 
				}
			}
		e<-eigen(c/sqrt(outer(dj,dj)))
		y[[j]]<-e$vectors[,nc]/sqrt(dj)
		}
	if (((fold-f) < eps) || (itel == itmax)) break
	itel<-itel+1; fold<-f
	}
return(list(f=f,y=y,rr=r^2,t=t))
}


burtTable<-function(data) {
nvar<-dim(data)[2]; ncat<-sapply(1:nvar,function(j) length(table(data[,j])))
tcat<-sum(ncat); first<-cumsum(c(1,ncat))[1:nvar]
burt<-matrix(0,tcat,tcat)
for (i in 1:nvar) {
	ii<-first[i]:(first[i]+ncat[i]-1)
	if (i < nvar) {
		for (j in (i+1):nvar) {
			jj<-first[j]:(first[j]+ncat[j]-1)
			burt[jj,ii]<-t(burt[ii,jj]<-as.matrix(table(data[,i],data[,j])))
			}
		}
	burt[ii,ii]<-as.matrix(table(data[,i],data[,i]))
	}
return(burt)
}

listRep<-function(x,m) {
y<-list()
for (j in 1:m) y<-c(y,list(x))
return(y)
}