#
#   gifi package
#   Copyright (C) 2006  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 2.0.0, 	2006-05-04,		first release in this form
#          							with voronoi plots and
#									quicktime movies

homals<-function(dframe,   # data (in data-frame)
	sets=0,              	# list of vectors of indices
	ndim=2,              	# dimensionality (default 2)
	active=TRUE,            # which variables are active (single TRUE means all)
	rank=ndim,           	# which quantification ranks (default all ndim)
	level="NO",          	# which quantification levels (default all nominal)            
	starplots=FALSE,        # which starplots (default none)
	catplots=FALSE,         # which category plots (default none)
	trfplots=FALSE,         # which transformation plots (default none)
	lossplots=FALSE,        # which loss plots (default none)
	hullplots=FALSE,		# which hullplots (default none)
	spanplots=FALSE,		# which spanning tree plots (default none)
	labplots=FALSE,			# which labeled object score plots (default none)
	vecplots=FALSE,	        # which vector plots (default none)
	prjplots=FALSE,	        # which projection plots (default none)
	vorplots=FALSE,	        # which voronoi plots (default none)
	graphplot=FALSE,        # graphplot (default no)
	objplot=FALSE,          # object score plot (default no)
	objwrite=FALSE,        	# object scores written to file (default no)
	objlabels=FALSE,        # object score plot labeled (default no)
	silent=FALSE,           # iteration history to stdout
	outwrite=TRUE,          # iterations and quantifications to file
	offset=1.20,         	# offset for labeled plots (default 1.20)
	eps1=-Inf,           	# iteration precision eigenvalues (default 1e-6)
	eps2=1e-6,           	# iteration precision eigenvectors (default 1e-6)
	itermax=100,         	# maximum number of iterations (default 100)
	saveMe=TRUE,            # do we return the results
	movieplot=FALSE,        # movieplot (default no)
	movielab=FALSE,			# labels in movieplot (default no)
	moviequick=FALSE,       # quicktime movie (default no)
	timer=FALSE             # time the steps of program (default no)
)               
{            
if (timer) stime<-proc.time()
name<-deparse(substitute(dframe))
nobj<-dim(dframe)[1]; nvar<-dim(dframe)[2]; iter<-0; pops<-0; 
for (j in 1:nvar) dframe[,j]<-as.factor(dframe[,j]) 
if (length(sets) == 1) sets<-lapply(1:nvar,"c")
pca<-max(sapply(sets,length)) <= 1
nset<-length(sets)
if (pca) mis<-apply(dframe,1,function (x) sum(ifelse(is.na(x),0,1)))
	else mis<-apply(sapply(1:nset,
		function(s) apply(cbind(dframe[,sets[[s]]]),1,
		function (x) prod(ifelse(is.na(x),0,1)))),1,sum)
vname<-attr(dframe,"names") 
rname<-attr(dframe,"row.names")
if (ndim == 1) starplots<-catplots<-lossplots<-graphplot<-objplot<-F
if (length(active)==1) active<-rep(active,nvar)
if (length(starplots)==1) starplots<-rep(starplots,nvar)
if (length(catplots)==1) catplots<-rep(catplots,nvar)
if (length(trfplots)==1) trfplots<-rep(trfplots,nvar)
if (length(lossplots)==1) lossplots<-rep(lossplots,nvar)
if (length(hullplots)==1) hullplots<-rep(hullplots,nvar)
if (length(spanplots)==1) spanplots<-rep(spanplots,nvar)
if (length(labplots)==1) labplots<-rep(labplots,nvar)
if (length(vecplots)==1) vecplots<-rep(vecplots,nvar)
if (length(prjplots)==1) prjplots<-rep(prjplots,nvar)
if (length(vorplots)==1) vorplots<-rep(vorplots,nvar)
doPDF<-any(starplots,catplots,trfplots,lossplots,hullplots,spanplots,labplots,vecplots,prjplots,vorplots,graphplot,objplot)
if (length(rank)==1) rank<-rep(rank,nvar)
if (length(level)==1) level<-rep(level,nvar)
for (j in 1:nvar) {
k<-length(levels(dframe[,j]))
	if (rank[j] > min(ndim,k-1)) rank[j]<-min(ndim,k-1)
	}
if (outwrite)
	outfile<-file(paste(name,"out",sep="."),"w")
if (objwrite)
	objfile<-file(paste(name,"obj",sep="."),"w")
x<-cbind(orthogonalPolynomials(mis,1:nobj,ndim))
x<-normX(centerX(x,mis),mis)$q
y<-lapply(1:nvar, function(j) computeY(dframe[,j],x))
if (timer) itime<-proc.time()
if (outwrite)
	cat("Iterations:\n",file=outfile)
if (moviequick) system(paste("mkdir ",name,"_movie",sep=""))
repeat {
	iter<-iter+1
	totsum<-array(0.0,dim(x))
	if (pca) up.y<-pcaUpdateY(dframe,x,y,totsum,active,rank,level,nset)
		else up.y<-ccaUpdateY(dframe,x,y,totsum,active,rank,level,sets)
	y<-up.y$y; totsum<-up.y$totsum
	qv<-normX(centerX((1/mis)*totsum,mis),mis)
	z<-qv$q;r<-qv$r;ops=sum(r);aps<-sum(La.svd(crossprod(x,mis*z),0,0)$d)/ndim	
	if (outwrite)
		cat("Iteration: ",formatC(iter,digits=3,width=3)," Eigenvalues: ",
			formatC(r,digits=6,width=9,format="f"),
			" Gain: ",formatC(c(ops,aps),digits=6,width=9, format="f"),"\n",file=outfile)
	if (!silent)
		cat("Iteration: ",formatC(iter,digits=3,width=3)," Eigenvalues: ",
			formatC(r,digits=6,width=9,format="f")," Gain: ",
			formatC(c(ops,aps),digits=6,width=9, format="f"),"\n")
	if (movieplot) myMoviePlot(x,z,movielab)
	if (moviequick) {
	    setwd(paste(name,"_movie",sep=""))
		jpeg(paste("file",as.character(iter),".jpg",sep=""))
		myMoviePlot(x,z,movielab)
		dev.off()
		setwd("..")
		}
	if (((ops - pops) < eps1) || ((1.0 - aps) < eps2) || (iter == itermax)) break 
		else {x<-z; pops<-ops}
	}
if (timer) ctime<-proc.time()
if (moviequick) {
	setwd(paste(name,"_movie",sep=""))
	system(paste("ffmpeg -i file%d.jpg ",name,".mov",sep=""))
	system("mv *.mov ..")
	system("rm *.jpg")
	setwd("..")
	system(paste("rm -rf ",name,"_movie",sep=""))
	}
if (doPDF)
	pdf(file=paste(name,"pdf",sep="."),encoding="MacRoman")
xlim<-c(min(z[,1]),max(z[,1])); if (ndim > 1) ylim<-c(min(z[,2]),max(z[,2]))
if (graphplot) if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
	graphplot(name,dframe,z[,c(s,t)],s,t); 
	}
if (objplot) if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
	objplot(name,rname,objlabels,offset,z[,c(s,t)],s,t); 
	}
ylist<-alist<-clist<-ulist<-NULL
for (j in 1:nvar) {
gg<-dframe[,j]; c<-computeY(gg,z); d<-as.vector(table(gg))
lst<-restrictY(d,c,rank[j],level[j])
y<-lst$y; a<-lst$a; u<-lst$z
ylist<-c(ylist,list(y)); alist<-c(alist,list(a)); clist<-c(clist,list(c)); ulist<-c(ulist,list(u))
if (outwrite) {
	writeHead(name,vname[j],ndim,active[j],rank[j],level[j],outfile)
	writeY(dframe[,j],y,"Y",outfile)
	writeY(dframe[,j],c,"C",outfile)
	writeY(dframe[,j],u,"Z",outfile)
	writeA(t(a),outfile)
	}
if (starplots[j]) {
	if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		starplot(name,vname[j],dframe[,j],y[,c(s,t)],z[,c(s,t)],s,t)}}
if (catplots[j]) {
	if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		catplot(name,vname[j],dframe[,j],y[,c(s,t)],xlim,ylim,s,t)}}
if (trfplots[j]) trfplot(name,vname[j],dframe[,j],y)
if (lossplots[j]) {
	if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		lossplot(name,vname[j],dframe[,j],y[,c(s,t)],z[,c(s,t)],s,t)}}
if (hullplots[j]) {
	if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		hullplot(name,vname[j],dframe[,j],y[,c(s,t)],z[,c(s,t)],s,t)}}
if (spanplots[j]) {
	require("ape",quietly=T)
	if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		spanplot(name,vname[j],dframe[,j],y[,c(s,t)],z[,c(s,t)],s,t)}}
if (labplots[j]) {
		if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		labplot(name,vname[j],dframe[,j],offset,z[,c(s,t)],s,t)}
		}
if (vorplots[j]) {
		if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		vorplot(name,vname[j],dframe[,j],offset,y[,c(s,t)],z[,c(s,t)],s,t)}
		}
if (vecplots[j] && (rank[j]==1)) {
		if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
		vecplot(name,vname[j],dframe[,j],z[,c(s,t)],y[,c(s,t)],a[c(s,t)],s,t,offset,objlabels)}}
if (prjplots[j] && (rank[j]==1)) {
		if (ndim > 1) for (s in 1:(ndim-1)) for (t in (s+1):ndim) {
			prjplot(name,vname[j],dframe[,j],z[,c(s,t)],y[,c(s,t)],a[c(s,t)],s,t,offset,objlabels)}}
			}
if (doPDF) dev.off()
if (objwrite){
	write(z,file=objfile,ncolumns=ndim)
	close(objfile)}
if (timer && outwrite) {otime<-proc.time()
	cat("Input time:     ",formatC(itime-stime,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Compute time:   ",formatC(ctime-itime,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Iteration time: ",formatC((ctime-itime)/iter,digits=6,width=10,format="f"),"\n",file=outfile)
	cat("Output time:    ",formatC(otime-ctime,digits=6,width=10,format="f"),"\n",file=outfile)}
if (outwrite) close(outfile)
if (saveMe) list(x=z,y=ylist,c=clist,a=alist,z=ulist)
}

sumSet<-function(g,n,p,y,set,active) {
z<-array(0.0,c(n,p))
for (j in set)
	if (active[j]){	    
	    gg<-g[,j]; ii<-which(!is.na(gg)); 
		z[ii,]<-z[ii,]+y[[j]][gg[ii],]}
return(z)
}

pcaUpdateY<-function(dframe,x,y,totsum,active,rank,level,nset){
for (j in 1:nset){
	if (active[j]){	    
		gg<-dframe[,j]; ycen<-computeY(gg,x); d<-as.vector(table(gg)); ii<-which(!is.na(gg))
		yy<-restrictY(d,ycen,rank[j],level[j])$y
		y[[j]]<-yy
		totsum[ii,]<-totsum[ii,]+yy[gg[ii],]
		}
	}
return(list(y=y,totsum=totsum))
}

ccaUpdateY<-function(dframe,x,y,totsum,active,rank,level,sets){
nobj<-dim(x)[1]; ndim<-dim(x)[2]; nset<-length(sets)
for (l in 1:nset) {
	indi<-sets[[l]]
	ss<-sumSet(dframe,nobj,ndim,y,indi,active)
	for (j in indi) {
		if (active[j]){	
			gg<-dframe[,j]; yy<-y[[j]]; ii<-which(!is.na(gg)); d<-as.vector(table(gg))
			ss<-ss-yy[gg[ii],]
			ycen<-computeY(dframe[,j],x-ss)
			yy<-restrictY(d,ycen,rank[j],level[j])$y
			ss<-ss+yy[gg[ii],]
			y[[j]]<-yy
			}
		}
	totsum<-totsum+ss
	}
return(list(y=y,totsum=totsum))
}

computeY<-function(g,x) apply(x, 2, function(z) tapply(z,g,mean))

restrictY<-function(d,y,r,level) {
if (sum(y^2) == 0) return(y)
switch(level,
"NO"=return(nominalY(d,y,r)),
"OR"=return(ordinalY(d,y,r)),
"NU"=return(numericalY(d,y,r)),
"PO"=return(polynomialY(d,y,r)))
}

nominalY<-function(d,y,r) {
qq<-La.svd(sqrt(d)*y,r,r)
zz<-(1/sqrt(d))*qq$u
aa<-qq$d[1:r]*qq$vt
list(y=zz%*%aa,z=zz,a=aa)
}

ordinalY<-function(d,y,r,itermax=100,eps=1e-6) {
qq<-La.svd(sqrt(d)*y,r,r)
a1=cbind(qq$vt[1,])
z1=cbind(orthogonalPolynomials(d,1:length(d),1))
if (r > 1) {
	a2<-cbind(qq$vt[2:r,])
	z2<-cbind(qq$u[,2:r])}
iter<-1; sold<-Inf
repeat{
	if (r == 1) ytilde<-y else ytilde<-y-(z2%*%t(a2))
	z1<-ytilde%*%a1/sum(a1*2)
	zo<-pava(z1,d)
	qo<-sum(d*zo^2)
	if (qo > 0) ao<-crossprod(ytilde,d*zo)/qo else ao<-a1
	so<-sum(d*(ytilde-zo%*%t(ao))^2)
	zp<-pava(-z1,d)
	qp<-sum(d*zp^2)
    if (qp > 0) ap<-crossprod(ytilde,d*zp)/qp else ap<-a1
	sp<-sum(d*(ytilde-zp%*%t(ap))^2)
	if (so < sp) {a1<-ao; z1<-zo; snew<-so}
		else {a1<-ap; z1<-zp; snew<-sp}
	if (r > 1) {
    	ytilde=y-(z1%*%t(a1))
    	qq<-La.svd(sqrt(d)*ytilde,r-1,r-1)
    	z2<-(1/sqrt(d))*qq$u
    	a2<-t(qq$d[1:(r-1)]*qq$vt)
    	snew<-sum(d*(ytilde-z2%*%t(a2))^2)}
	if ((iter == itermax) || ((sold - snew) < eps)) break()
		else {iter<-iter+1; sold<-snew}
}
if (r > 1) z<-cbind(z1,z2) else z<-cbind(z1)
z<-weightedGramSchmidt(z,d)$pol
a<-crossprod(y,d*z)
list(yhat=z%*%t(a),z=z,a=a)
}

numericalY<-function(d,y,r) {
z0<-orthogonalPolynomials(d,1:length(d),1)
a0<-as.vector(crossprod(z0*d,y))
if (r == 1)
	return(list(y=z0%o%a0,z=cbind(z0),a=rbind(a0)))
else {
	yy<-y-z0%o%a0
	qq<-La.svd(sqrt(d)*yy,r,r)
	zz<-cbind(z0,array((1/sqrt(d))*qq$u[,1:(r-1)],c(dim(y)[1],r-1)))
	aa<-rbind(a0,array((qq$d[1:r]*qq$vt)[1:(r-1),],c(r-1,dim(y)[2])))
	return(list(y=zz%*%aa,z=zz,a=aa))}
}

polynomialY<-function(d,y,r) {
k<-length(d)
zz<-orthogonalPolynomials(d,1:k,min(r,k-1))
aa<-crossprod(zz,d*y)
list(y=zz%*%aa,z=zz,a=aa)
}

writeHead<-function(name,vname,p,a,r,c,ofile) {
s<-NULL; sl<-35+10*p; for (i in 1:sl) s<-paste(s,"*",sep="")
cat("\n",formatC(s,format="s"),"\n",file=ofile,sep="")
cat(formatC(vname,format="s"),
if (a) "(Active)," else "(Passive),",
"Rank =", formatC(r,width=1),"and",
"Level is", 
switch(c,"NO"="nominal","OR"="ordinal","NU"="numerical","PO"="polynomial"),
"\n",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
}

writeY<-function(g,y,t,ofile) {
d<-as.vector(table(g));l<-levels(g)
s<-NULL; sl<-35+10*dim(y)[2]; for (i in 1:sl) s<-paste(s,"*",sep="")
u<-NULL; sl<-35+10*dim(y)[2]; for (i in 1:sl) u<-paste(u,"-",sep="")
switch(t,"Z"=cat("Lower Rank Quantifications\n",file=ofile),"C"=cat("Category Centroids\n",file=ofile),
"Y"=cat("Rank-restricted Category Quantifications\n",file=ofile))
cat(formatC(u,format="s"),"\n",file=ofile)
for (k in 1:length(d))
cat(formatC(l[k],width=10,format="s"),
formatC(d[k],width=5)," *** ",
formatC(y[k,],digits=6,width=10,format="f")," *** ",
formatC(sum(d[k]*y[k,]^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
cat(formatC(" ",format="s",width=10),
formatC(sum(d),width=5)," *** ",
formatC(d%*%(y^2),digits=6,width=10,format="f")," *** ",
formatC(sum(d%*%(y^2)),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
}

writeA<-function(a,ofile) {
s<-NULL; sl<-35+10*dim(a)[2]; for (i in 1:sl) s<-paste(s,"*",sep="")
u<-NULL; sl<-35+10*dim(a)[2]; for (i in 1:sl) u<-paste(u,"-",sep="")
cat("Category Loadings\n",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
for (k in 1:dim(a)[1])
cat(formatC(" ",format="s",width=10),formatC(k,width=5)," *** ",
formatC(a[k,],digits=6,width=10,format="f")," *** ",
formatC(sum(a[k,]^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(u,format="s"),"\n",file=ofile)
cat(formatC(" ",format="s",width=15)," *** ",
formatC(apply(a^2,2,sum),digits=6,width=10,format="f")," *** ",
formatC(sum(a^2),digits=6,width=10,format="f"),"\n",sep="",file=ofile)
cat(formatC(s,format="s"),"\n",file=ofile)
} 
catplot<-function(name,vname,g,y,xx,yy,s,t) {
plot(y,type="l",xlim=xx,ylim=yy,main=paste("Category plot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t))
text(y,levels(g))
abline(h=0)
abline(v=0)
}

trfplot<-function(name,vname,g,y) {
yy<-c(min(y),max(y))
p<-dim(y)[2]
r<-rainbow(p-1)
plot(y[,1],ylim=yy,type="l",
main=paste("Transformation plot for",name,":",vname),
xlab="original",ylab="transformed")
text(y[,1],levels(g))
if (p > 1)
	for (s in 2:p) {lines(y[,s],col=r[s-1]); text(y[,s],levels(g))}
abline(h=0)
}

starplot<-function(name,vname,g,y,x,s,t) {
plot(x,col="GREEN",pch=8,main=paste("Starplot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t))
z<-computeY(g, x)
points(z,type="n")
text(z,levels(g),col="RED")
for (i in 1:length(g)) lines(rbind(x[i,],z[g[i],]),col="BLUE")
}

spanplot<-function(name,vname,g,y,x,s,t) {
plot(x,col="GREEN",pch=8,main=paste("Spanplot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t))
lev<-levels(g)
rb<-rainbow(length(lev))
for (k in lev) {
	ind<-which(k==g)
	n<-length(ind)
	mm<-mst(dist(x[ind,]))
	for (i in 1:n) {
		jnd<-which(1==as.vector(mm[i,]))
		sapply(jnd,function(r) lines(rbind(x[ind[i],],x[ind[r],]),col=rb[which(lev==k)]))
		}
	}
}

hullplot<-function(name,vname,g,y,x,s,t) {
plot(x,col="GREEN",pch=8,main=paste("Hullplot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t))
for (k in levels(g)) {
	ind<-which(k==g)
	lst<-ind[chull(x[ind,])]
	lines(x[c(lst,lst[1]),])
	text(x[lst,],k) 
	}
}

lossplot<-function(name,vname,g,y,x,s,t) {
z<-computeY(g, x); k<-dim(z)[1]
xx<-c(min(c(z[,1],y[,1])),max(c(z[,1],y[,1])))
yy<-c(min(c(z[,2],y[,2])),max(c(z[,2],y[,2])))
plot(y,type="n",main=paste("Lossplot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=xx,ylim=yy)
text(y,levels(g),col="RED");lines(y,col="RED")
points(z,type="n")
text(z,levels(g),col="GREEN");lines(z,col="GREEN")
for (i in 1:k) lines(rbind(y[i,],z[i,]),col="BLUE")
abline(h=0)
abline(v=0)
}

labplot<-function(name,vname,g,offset,x,s,t) {
xx<-c(min(x[,1]),max(x[,1]))
yy<-c(min(x[,2]),max(x[,2]))
plot(x,type="n",main=paste("Labeled plot for",name,":",vname),
xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=offset*xx,ylim=offset*yy)
text(x,as.vector(g))
}

vorplot<-function(name,vname,g,offset,y,x,s,t) {
z<-rbind(x,y)
xx<-c(min(z[,1]),max(z[,1]))
yy<-c(min(z[,2]),max(z[,2]))
plot(x,type="n",main=paste("Voronoi plot for",name,":",vname),
	xlab=paste("dimension",s),ylab=paste("dimension",t),
	xlim=offset*xx,ylim=offset*yy)
drawEdges(y,far=1000)
text(x,as.vector(g))
}

vecplot<-function(name,vname,g,x,y,a,s,t,offset,objlabel) {
xx<-c(min(x[,1]),max(x[,1]))
yy<-c(min(x[,2]),max(x[,2]))
zz<-c(min(xx[1],yy[1]),max(xx[2],yy[2]))
if (objlabel) {
	plot(x,type="n",main=paste("Vectorplot for",name,":",vname),
	xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=offset*zz,ylim=offset*zz)
	text(x,as.vector(g))}
else {
	plot(x,col="GREEN",pch=8,main=paste("Vectorplot for",name,":",vname),
	xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=zz,ylim=zz)
	}
text(y,levels(g),col="RED")
slope=a[2]/a[1]
abline(coef=c(0,slope))
for (i in 1:length(g)) 
	{
	xs<-xe<-x[i,]; xe[1]<-(xs[1]+(xs[2]*slope))/(1+(slope^2)); xe[2]<-slope*xe[1]
	lines(rbind(xs,xe),col="BLUE")
	}
abline(h=0)
abline(v=0)
}

prjplot<-function(name,vname,g,x,y,a,s,t,offset,objlabel) {
xx<-c(min(x[,1]),max(x[,1]))
yy<-c(min(x[,2]),max(x[,2]))
zz<-c(min(xx[1],yy[1]),max(xx[2],yy[2]))
if (objlabel) {
	plot(x,type="n",main=paste("Projection plot for",name,":",vname),
	xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=offset*zz,ylim=offset*zz)
	text(x,as.vector(g))}
else {
	plot(x,col="GREEN",pch=8,main=paste("Projection Plot for",name,":",vname),
	xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=zz,ylim=zz)
	}
text(y,levels(g),col="RED")
slope=a[2]/a[1]
abline(coef=c(0,slope))
slope=-a[1]/a[2]
for (j in 1:dim(y)[1])
	abline(coef=c(y[j,2]-slope*y[j,1],slope),col="RED")
for (i in 1:dim(x)[1])
	{
	j<-g[i]; icpt<-y[j,2]-slope*y[j,1]
	u<-(x[i,1]+slope*(x[i,2]-icpt))/(1+slope^2)
	lines(rbind(x[i,],c(u,icpt+slope*u)),col="BLUE")
	}
abline(h=0)
abline(v=0)
}

graphplot<-function(name,g,x,s,t) {
plot(x,col="GREEN",pch=8,main=paste("Graphplot for",name),
xlab=paste("dimension",s),ylab=paste("dimension",t))
for (j in 1:dim(g)[2]){
	y<-computeY(g[,j], x)
	points(y,col="RED",pch=16)
	for (i in 1:dim(g)[1]) lines(rbind(x[i,],y[g[i,j],]))}
}

objplot<-function(name,rname,objlabel,offset,x,s,t) {
xx<-c(min(x[,1]),max(x[,1]))
yy<-c(min(x[,2]),max(x[,2]))
if (objlabel) {
plot(x,type="n",main=paste("Object score plot for",name),
xlab=paste("dimension",s),ylab=paste("dimension",t),xlim=offset*xx,ylim=offset*yy)
text(x,rname)}
else {
plot(x,pch=8,main=paste("Object score plot for",name),
xlab=paste("dimension",s),ylab=paste("dimension",t)) }
}

myMoviePlot<-function(x,y,movielab=FALSE) {
plot(x,xlab="dimension 1",ylab="dimension 2",type="n")
if (movielab) {
	text(x,as.character(1:dim(x)[1]),pch="+",col="GREEN")
	text(y,as.character(1:dim(x)[1]),pch="+",col="RED")
	}
else {
	points(x,pch="+",col="GREEN")
	points(y,pch="+",col="RED")
	}
sapply(1:dim(x)[1],function(i) lines(matrix(c(x[i,1:2],y[i,1:2]),2,2,byrow=TRUE)))
}

weightedGramSchmidt<-function(x,w) {
ss<-NULL; 
for (j in 1:dim(x)[2]) {
if (j > 1) {xx<-x[,1:(j-1)]; x[,j]<-x[,j]-xx%*%(crossprod(xx,(w*x[,j])))}
s<-sqrt(sum(w*x[,j]^2)); ss<-c(ss,s); x[,j]<-x[,j]/s;}
list(pol=x,fac=ss)
}

orthogonalPolynomials<-function(w,x,p) {
z<-weightedGramSchmidt(outer(x,0:p,"^"),w)$pol[,2:(p+1)]
}

centerX<-function(x,w) apply(x,2,function(z) z-weighted.mean(z,w))

normX<-function(x,w) {
qq<-La.svd(sqrt(w)*x); list(q=(1/sqrt(w))*(qq$u),r=qq$d)}

makeEdges<-function(a) {
c<-(rowSums(a^2))/2
n<-dim(a)[1]; lns<-matrix(0,0,8)
for (i in 1:(n-1)) {
	for(j in (i+1):n) {
	    dd<-a[i,]-a[j,]; dc<-c[i]-c[j]; ss<-sum(dd^2)
	    if (is.nul(ss)) next() 
	    ee<-dc*dd/ss; ff<-c(-dd[2],dd[1])
		xlw<--Inf; xup<-Inf
		for (k in (1:n)[-c(i,j)]) {
		    dd<-a[i,]-a[k,]; dc<-c[i]-c[k]
		    mum<-sum(dd*ff); mom<-dc-sum(dd*ee)
		    if (is.nul(mum) & (mom > 0)) {
		    	xlw<-Inf; xup<--Inf
		    	}
		    if (mum>0) xlw<-max(xlw,mom/mum)
		    if (mum<0) xup<-min(xup,mom/mum)
			}
		if (xlw<xup) lns<-rbind(lns,c(i,j,ee,ff,xlw,xup))
		}
	}
return(lns)
}

drawEdges<-function(a,far=1000) {
lns<-makeEdges(a)
p<-dim(lns)[1]; n<-dim(a)[1]
for (i in 1:p) {
	ee<-lns[i,3:4]; ff<-lns[i,5:6]
	xlw<-lns[i,7]; if (xlw == -Inf) xlw<--far
	plw<-ee+xlw*ff
	xup<-lns[i,8]; if (xup == Inf) xup<-far
	pup<-ee+xup*ff
	lines(rbind(plw,pup),col="RED")
	}
}
	
is.nul<-function(x) {
return(abs(x)<1e-12)
}

pava<-function(x,w=rep(1,length(x)),block=weighted.mean){
nblock<-n<-length(x); blocklist<-array(1:n,c(n,2)); blockvalues<-x; active<-1
repeat{
	if (!isUpSatisfied(blockvalues,active)) {
		blockmerge<-mergeBlockup(blocklist,blockvalues,x,w,active,block)
		blockvalues<-blockmerge$v; blocklist<-blockmerge$l
		nblock<-nblock-1
		while (!isDownSatisfied(blockvalues,active)) {
			blockmerge<-mergeBlockup(blocklist,blockvalues,x,w,active-1,block)
			blockvalues<-blockmerge$v; blocklist<-blockmerge$l; 
			nblock<-nblock-1; active<-active-1;
			}
		}
	else if (active == nblock) break() else active<-active+1
	}	
putBack(n,blocklist,blockvalues)
}

mergeBlockup<-function(blocklist,blockvalues,x,w,i,block){
n<-length(blockvalues); nn<-1:n; ii<-which(i+1!=nn)
blocklist[i,]<-c(blocklist[i,1],blocklist[i+1,2])
indi<-blocklist[i,1]:blocklist[i+1,2]
blockvalues[i]<-block(x[indi],w[indi])
blocklist<-blocklist[ii,]
if (length(ii) == 1) dim(blocklist)<-c(1,2)
blockvalues<-blockvalues[ii]
list(v=blockvalues,l=blocklist)
}

putBack<-function(n,blocklist,blockvalues){
x<-rep(0,n);nb<-length(blockvalues)
for (i in 1:nb) {
		x[blocklist[i,1]:blocklist[i,2]]<-blockvalues[i]}
return(x)
}

isUpSatisfied<-function(x,i) (i == length(x))||(x[i]<=x[i+1])

isDownSatisfied<-function(x,i) (i == 1)||(x[i-1]<=x[i])