require(gifi)

logitFold<-function(
	data,					# data (in data frame)
	ndim=2,                 # dimensionality (default 2)
	rank=ndim,           	# which quantification ranks (default all ndim)
	level="NO",          	# which quantification levels (default all nominal)            
	outeps=1e-4,
	ineps=1e-4,
	reg=1e-6,
	rule="dist",
	outmax=500,
	inmax=1,
	objplot=TRUE,
	voronoi=FALSE) {
nobj<-dim(data)[1]; nvar<-dim(data)[2]; pdev<-Inf
name<-deparse(substitute(data))
if (length(rank)==1) rank<-rep(rank,nvar)
if (length(level)==1) level<-rep(level,nvar)
xy<-homals(data,ndim=ndim,save=TRUE,outwrite=FALSE,silent=TRUE); x<-xy$x; yh<-xy$y
y<-matrix(0,0,ndim); wgt<-matrix(0,nobj,0); itel<-1
for (j in 1:nvar) {
	y<-rbind(y,yh[[j]])
	u<-as.factor(data[,j]); h<-levels(u); k<-length(h)
	wgt<-cbind(wgt,ifelse(is.na(outer(u,h,"==")),0,1))
	}
repeat {
	targ<-matrix(0,nobj,0); dist<-matrix(0,nobj,0); prob<-matrix(0,nobj,0)
	jind<-0; dev<-0; mdd<-0
	for (j in 1:nvar) {
		u<-as.factor(data[,j]); h<-levels(u); k<-length(h)
		gj<-ifelse(outer(u,h,"=="),1,0)
		gj<-ifelse(is.na(gj),0,gj)
		yj<-y[(jind+1):(jind+k),]
		if (rule=="dist")
			dj<-distRect(x,yj,reg)
		if (rule=="dist2")
			dj<-dist2Rect(x,yj)
		pj<-exp(-dj); pj<-pj/rowSums(pj)
		md<-gj-pj; mdd<-mdd+sum(md^2)
		targ<-cbind(targ,dj-(2*md))
		dist<-cbind(dist,dj); prob<-cbind(prob,pj)
		jind<-jind+k; dev<-dev-sum(gj*log(pj))
		}
	dev<-dev+mdd
	if (rule=="dist")
		uu<-smacofRect(targ,w=wgt,init=list(x,y),itmax=inmax,reg=reg,verbose=FALSE)
	if (rule=="dist2")
		uu<-alscalRect(targ,w=wgt,init=list(x,y),itmax=inmax,verbose=FALSE,ithist=FALSE)
	cat("Iteration: ",formatC(itel,digits=6,width=6),
	    "Deviance: ",formatC(dev+mdd,digits=6,width=12,format="f"),
	    "\n")
	if ((itel == outmax) || ((pdev-dev)<outeps)) break()
	itel<-itel+1; pdev<-dev
	x<-uu$x; y<-uu$y
	}
if (objplot) {
	pdf(file=paste("objplot_",rule,"_",name,".pdf",sep=""),encoding="MacRoman")
	plot(x,type="n",xlab="dimension 1",ylab="dimension 2")
	title(main=paste("Objects",name))
	text(x,as.character(labels(data)[[1]]))
	dev.off()
	}
if (voronoi) {
	pdf(file=paste("voronoi_",rule,"_",name,".pdf",sep=""),encoding="MacRoman"); jind<-0
	mm1<-max(rbind(x,y)); mm2<-min(rbind(x,y)); mm<-c(mm2,mm1)
	for (j in 1:nvar) {
		u<-as.factor(data[,j]); h<-levels(u); k<-length(h); yj<-y[(jind+1):(jind+k),]
		plot(deldir(yj[,1],yj[,2],rw=c(mm,mm)), wlines = "tess", wpoints = "none", 
			col = c(1, 2, 1), lty = c(1, 1),
			xlab="dimension 1",ylab="dimension 2")
		title(main=paste("Voronoi",name,as.character(labels(data)[[2]][j])))
		text(yj,as.character(h),col="red"); text(x,as.character(data[,j]))
		jind<-jind+k
		}
	dev.off()
	}
list(x=x,y=x,dev=dev)
}

smacofRect<-function(
	delta,
	w=matrix(1,dim(delta)[1],dim(delta)[2]),
	p=2,
	init="svd",
	verbose=FALSE,
	itmax=100,
	reg=1e-6,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
delta<-ifelse(is.na(delta),0,delta)
delta_plus<-ifelse(delta>=0,delta,0); delta_min<-ifelse(delta<=0,-delta,0)
if (is.list(init))
	{
	x<-init[[1]]; y<-init[[2]]
	}
else {
	e<-delta_plus^2; e<--0.5*(e-outer(rowSums(e)/m,colSums(e)/n,"+")+(sum(e)/(n*m)))
	z<-svd(e,nu=p,nv=0); x<-z$u;  y<-crossprod(e,x)
	}
d<-distRect(x,y,reg); lold<-sum(w*(delta-d)^2)
repeat {
	ww<-w*(1+(delta_min/d)); wr<-rowSums(ww); wc<-colSums(ww)
	v<-solve(diag(wc)+(1/m)-crossprod(ww,ww/wr))-(1/m)
	b<-w*delta_plus/d; br<-rowSums(b); bc<-colSums(b)
	xraw<-(br*x)-(b%*%y); yraw<-(bc*y)-crossprod(b,x)
	y<-v%*%(yraw+crossprod(ww,xraw/wr)); x<-(xraw+(ww%*%y))/wr; d<-distRect(x,y,reg)
    lnew<-sum(w*(delta-d)^2)
	if (verbose) {
		cat("Iteration: ",formatC(itel,digits=6,width=6),
	    	"   Stress:   ",formatC(lold,digits=6,width=12,format="f"),
	    	" ==>",formatC(lnew,digits=6,width=12,format="f"),"\n")
	    }
	if (((lold-lnew)<eps) || (itel==itmax)) break()
	lold<-lnew; itel<-itel+1
	}
return(list(x=x,y=y,d=d,l=lnew))
}

met2Fold<-function(
	delta,
	w=matrix(1,dim(delta)[1],dim(delta)[2]),
	p=2,
	init="svd",
	verbose=FALSE,
	itmax=100,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; nn<-1:n; mm<-n+(1:m); itel<-1
delta<-ifelse(is.na(delta),0,delta); c<-matrix(0,n+m,n+m)
if (is.list(init))
	{
	x<-init[[1]]; y<-init[[2]]
	}
else {
	e<--0.5*(delta-outer(rowSums(delta)/m,colSums(delta)/n,"+")+(sum(delta)/(n*m)))
	z<-svd(e,nu=p,nv=0); x<-z$u;  y<-crossprod(e,x)
	}
d<-dist2Rect(x,y); lold<-sum(w*(delta-d)^2); bnd<-n+m+2
b<-w*delta; br<-rowSums(b); bc<-colSums(b)
repeat {
    bd<-w*d; bdr<-rowSums(bd); bdc<-colSums(bd)
    c[nn,mm]<-(x%*%t(y))-(b-bd)/bnd
	c[mm,nn]<-(y%*%t(x))-t(b-bd)/bnd
	c[nn,nn]<-(x%*%t(x))+diag(br-bdr)/bnd
	c[mm,mm]<-(y%*%t(y))+diag(bc-bdc)/bnd
    ei<-eigen(c); z<-ei$vectors[,1:p]
	z<-z*outer(rep(1,n+m),sqrt(ei$values[1:p]))
    x<-z[nn,]; y<-z[mm,]; d<-dist2Rect(x,y); lnew<-sum(w*(delta-d)^2)
	if (verbose) {
		cat("Iteration: ",formatC(itel,digits=6,width=6),
	    	"   Stress:   ",formatC(lold,digits=6,width=12,format="f"),
	    	" ==>",formatC(lnew,digits=6,width=12,format="f"),"\n")
	    }
	if (((lold-lnew)<eps) || (itel==itmax)) break()
	lold<-lnew; itel<-itel+1
	}
return(list(x=x,y=y,d=d,l=lnew))
}

alscalRect<-function(
	delta,
	w=matrix(1,dim(delta)[1],dim(delta)[2]),
	p=2,
	init="svd",
	verbose=FALSE,
	ithist=TRUE,
	itmax=100,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
nn<-1:n; mm<-1:m; pp<-1:p; wr<-rowSums(w); wc<-colSums(w)
if (is.list(init))
	{
	x<-init[[1]]; y<-init[[2]]
	}
else {
	e<--0.5*(delta-outer(rowSums(delta)/m,colSums(delta)/n,"+")+(sum(delta)/(n*m)))
	z<-svd(e,nu=p,nv=0); x<-z$u;  y<-crossprod(e,x)
	}
d<-dist2Rect(x,y); r<-w*(delta-d); t<-r*(delta-d)
rr<-rowSums(r); rc<-colSums(r)
tr<-rowSums(t); tc<-colSums(t); sstress<-sum(tc); psstress=sstress
repeat {
    for (i in nn) for (s in pp) {
	    u<-x[i,s]-y[,s]; ww<-w[i,]
		a4<-wr[i]; a3<-4*sum(ww*u)
		a2<-(4*sum(ww*u*u))-(2*rr[i])
		a1<--4*sum(r[i,]*u); a0<-sstress
		mn<-minimizeQuartic(c(a0,a1,a2,a3,a4)); th<-mn[1]
		x[i,s]<-x[i,s]+th; d[i,]<-d[i,]+(2*th*u)+(th*th)
		r[i,]<-ww*(delta[i,]-d[i,]); t[i,]<-r[i,]*(delta[i,]-d[i,])
		rr[i]<-sum(r[i,]); tr[i]<-sum(t[i,]); sstress=sum(tr)
		if (verbose) print(sstress)
		}
	rc<-colSums(r); tc<-colSums(t)
    for (j in mm) for (s in pp) {
	    u<-x[,s]-y[j,s]; ww<-w[,j]
		a4<-wc[j]; a3<--4*sum(ww*u)
		a2<-(4*sum(ww*u*u))-(2*rc[j])
		a1<-4*sum(r[,j]*u); a0<-sstress
		mn<-minimizeQuartic(c(a0,a1,a2,a3,a4)); th<-mn[1]
		y[j,s]<-y[j,s]+th; d[,j]<-d[,j]-(2*th*u)+(th*th)
		r[,j]<-ww*(delta[,j]-d[,j]); t[,j]<-r[,j]*(delta[,j]-d[,j])
		rc[j]<-sum(r[,j]); tc[j]<-sum(t[,j]); sstress=sum(tc)
		if (verbose) print(sstress)
		}
	rr<-rowSums(r); tr<-rowSums(t)
	if (ithist) 	{
		cat("Iteration: ",formatC(itel,digits=6,width=6),
			"   SStress:   ",formatC(psstress,digits=6,width=12,format="f"),
			" ==>",formatC(sstress,digits=6,width=12,format="f"),"\n")
		}
	if (((psstress-sstress)<eps) || (itel==itmax)) break()
	psstress<-sstress; itel<-itel+1
}
return(list(x=x,y=y,d=d,sstress=sstress))
}


minimizeQuartic<-function(a) {
#
# D.J. Jeffrey: Formulae, algorithms, and quartic extrema. 
# Mathematics Magazine 70, pp349 - 356 (1997).
# http://www.apmaths.uwo.ca/~djeffrey/Offprints/matmag97.ps 
#
a0<-a[1]; a1<-a[2]; a2<-a[3]; a3<-a[4]; a4<-a[5]
if (!is.pos(a4)) stop("no minimum")
b2<-(a2/(3*a4))-((a3*a3)/(8*a4*a4))
b1<-(a1/(2*a4))+((a3*a3*a3)/(16*a4*a4*a4))-((a2*a3)/(4*a4*a4))
if (is.nul(b1)) {
    bm<-min(0,b2)
	return(c(sqrt(-3*bm/2),-9*bm*bm/4))
	}
b0<-((a2*a3*a3)/(16*a4*a4))-((a1*a3)/(4*a4))-((3*a3*a3*a3*a3)/(256*a4*a4*a4))
ss<-(b1*b1)+(b2*b2*b2)+sqrt(as.complex((b1*b1*b1*b1)+(2*b1*b1*b2*b2*b2)))
s3<-ss^(1/3); kf<-s3+(b2*b2/s3)+b2; kf<-Re(kf)
gg<--(3/4)*(kf-b2)*(kf-(3*b2)); ff<-(a4*gg)+a0+b0; th<--b1/kf-(a3/(4*a4))
return(c(th,ff))
}

is.nul<-function(x) {
return(abs(x) < 1e-10)
}

is.pos<-function(x) {
return(x > 1e-10)
}

is.neg<-function(x) {
return(x < -1e-10)
}

distRect<-function(x,y,reg) {
return(sqrt(reg+outer(rowSums(x^2),rowSums(y^2),"+")-2*x%*%t(y)))
}

dist2Rect<-function(x,y) {
return(outer(rowSums(x^2),rowSums(y^2),"+")-2*x%*%t(y))
}