metFold<-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]; itel<-1
delta_plus<-ifelse(delta>=0,delta,0); delta_min<-ifelse(delta<=0,-delta,0)
if (init=="svd") {
	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)
	}
if (init=="rnd") {
	x<-matrix(rnorm(n*p),n,p)
	y<-matrix(rnorm(m*p),m,p)
	}
d<-distRect(x,y); lold<-sum(w*(delta-d)^2)
repeat {
	dd<-ifelse(d>0,d,1)
	ww<-w*(1+(delta_min/dd)); wr<-rowSums(ww); wc=colSums(ww)
	v<-solve(diag(wc)+(1/m)-crossprod(ww,ww/wr))-(1/m)
	b<-w*delta_plus/dd; 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)
    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)) return(list(x=x,y=y,d=d,l=lnew))
	lold<-lnew; itel<-itel+1
	}
}


extFold<-function(
	delta,
	xbase,
	w=matrix(1,dim(delta)[1],dim(delta)[2]),
	p=2,
	verbose=FALSE,
	itmax=100,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
qbase<-qr(xbase); qx<-qr.Q(qbase); qr<-qr.R(qbase)
e<-delta^2; e<--0.5*(e-outer(rowSums(e)/m,colSums(e)/n,"+")+(sum(e)/(n*m)))
bx<-svd(crossprod(e,qx))$v[,1:p]; x<-qx%*%bx; y<-crossprod(e,x)
d<-distRect(x,y); lold<-sum(w*(delta-d)^2)
repeat {
	dd<-ifelse(d>0,d,1)
	ww<-w*(1+(delta_min/dd)); wr<-rowSums(ww); wc=colSums(ww)
	v<-solve(diag(wc)+(1/m)-crossprod(ww,ww/wr))-(1/m)
	b<-w*delta_plus/dd; 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)
    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)) return(list(x=x,y=y,d=d,l=lnew))
	lold<-lnew; itel<-itel+1
	}
}

initFold<-function(
	delta,
	w=matrix(1,dim(delta)[1],dim(delta)[2]),
	p=2,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
wr<-rowSums(w); wc=colSums(w)
v<-solve(diag(wc)+(1/m)-(t(w)%*%diag(1/wr)%*%w))-(1/m)
e<-delta^2; e<--0.5*(e-outer(rowSums(e)/m,colSums(e)/n,"+")+(sum(e)/(n*m)))
z<-svd(e,nu=p,nv=0); a<-z$u;  b<-cbind(crossprod(e,a),rep(1,m))
uu<-crossprod(w%*%b,a); vv<-crossprod(b,wc*b); ww<-crossprod(a,wr*a)
s<-diag(p); t<-rbind(s,rep(0,p)); x<-a%*%s; y<-b%*%t
d<-distRect(x,y); lold<-sum(w*(delta-d)^2)
repeat {
	bm<-as.matrix(w*delta/d); br<-rowSums(bm); bc<-colSums(bm)
	xraw<-(outer(br,rep(1,p))*x)-(bm%*%y); yraw<-(outer(bc,rep(1,p))*y)-(t(bm)%*%x)
	ysma<-v%*%(yraw+crossprod(w,xraw/wr)); xsma<-(xraw+(w%*%y))/wr
	dsma<-distRect(xsma,ysma); lsma<-sum(w*(delta-dsma)^2)
	tt<-crossprod(b,(wc*ysma)-crossprod(w,xsma)); ss<-crossprod(a,(wr*xsma)-w%*%ysma)
	rhs<-ss-crossprod(uu,solve(vv,tt)); lhs<-ww-crossprod(uu,solve(vv,uu))
	s<-solve(lhs,rhs); t<-solve(vv,tt-uu%*%s); x<-a%*%s; y<-b%*%t
	d<-distRect(x,y); lnew<-sum(w*(delta-d)^2)
	cat("Iteration: ",formatC(itel,digits=6,width=6),
	    "   Stress Old: ",formatC(lold,digits=6,width=12,format="f"),
	    "   Stress Free:",formatC(lsma,digits=6,width=12,format="f"),
	    "   Stress New:",formatC(lnew,digits=6,width=12,format="f"),
		"\n")
	if ((lold-lnew)<eps) return(list(x=x,y=y,d=d,l=lnew))
	lold<-lnew; itel<-itel+1
	}
}

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

krusFold<-function(ranks,p=2,eps1=1e-6,eps2=1e-6){
n<-dim(ranks)[1]; m<-dim(ranks)[2]; disp<-ranks; 
q<-m*(m+1)*(m-1)/12; s<-rep(0,n); t<-rep(q,n); lold<-Inf; x<-FALSE; y<-FALSE
repeat{
	uf<-metfold(disp,w=outer(1/t,rep(0,m),"+"),x=x,y=y,p=p,eps=eps1)
	d<-uf$d; x<-uf$x; y<-uf$y; luf<-uf$l
	for (i in 1:n) {
		ind<-as.vector(as.matrix(ranks[i,])); ord<-order(ind); dist<-d[i,]-mean(d[i,])
		dhat<-pava(dist[ord])[ind]; alpha=sum(dist^2)/sum(dist*dhat);dhat<-dhat*alpha
		s[i]<-sum((dist-dhat)^2); t[i]<-sum(dhat^2)
		disp[i,]<-dhat+mean(d[i,])
		}
	l<-sum(s/t)
	print(c(lold,luf,l))
	if ((lold-l)<eps2) return(list(x=x,y=y,d=d))
	lold<-l
	}
}


logitFold<-function(g,ndim=2,outeps=1e-6,ineps=1e-6,outmax=100,inmax=1) {
n<-dim(g)[1]; m<-dim(g)[2]; pdev<-Inf
xy<-homals(g,ndim=ndim,save=TRUE); x<-xy$x; yh<-xy$y
y<-matrix(0,0,ndim); for (j in 1:m) y<-rbind(y,yh[[j]]); itel<-1
repeat {
	targ<-matrix(0,n,0); dist<-matrix(0,n,0); prob<-matrix(0,n,0)
	jind<-0; dev<-0; mdd<-0
	for (j in 1:m) {
		u<-as.factor(g[,j]); h<-levels(u); k<-length(h)
		gj<-ifelse(outer(u,h,"=="),1,0)
		yj<-y[(jind+1):(jind+k),]; dj<-distRect(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))
		}
	sink("targ"); print(targ); sink()
	uu<-metFold(targ,itmax=inmax)
	cat("Iteration: ",formatC(itel,digits=6,width=6),
	    "Deviance: ",formatC(dev,digits=6,width=12,format="f"),
	    "\n")
	if ((itel == outmax) || ((pdev-dev)<outeps)) return(list(x=x,y=y))
	itel<-itel+1; pdev<-dev
	x<-uu$x; y<-uu$y
	}
}

binFold<-function(data,
	rownames=as.character(1:dim(data)[1]),
	rowlabs=as.character(1:dim(data)[1]),
	ndim=2,
	eps=1e-3,
	imax=10000,
	correct=FALSE,
	extreme=FALSE,
	save=FALSE,
	cutoff=1e-3,
	form="logit",
	offset=1.20) {
name<-deparse(substitute(data))
outfile<-file(paste(name,"out",sep="."),"w")
vlog<-function(a,b) ifelse(b>0,a*log(b),0)
ls<--2*length(which(!is.na(data)))*log(.5)
itel<-1
z<-ifelse(is.na(data),0,-4*(data-.5))
repeat{
	a<-apply(z,2,mean)
	z<-z-outer(rep(0,dim(z)[1]),a,"+")
	sv<-La.svd(z,nu=ndim,nv=ndim,method="dgesdd")
	x<-sv$u; y<-sv$d[1:ndim]*(sv$v)
	aa<-(x%*%y)+outer(rep(0,dim(z)[1]),a,"+")
	if (form=="logit") pr<-1/(1+exp(-aa)) else {
		pr<-pnorm(aa); pd<-dnorm(aa)
		}
	if (correct)
		tb<-table(as.vector(data),as.vector(ifelse(pr>.5,1,0)))
	if (extreme)
		{
		pe<-pr[which(!is.na(data))]
		pp<-length(which((pe>(1-cutoff))|(pe<cutoff)))/length(pe)
		}
	lk<-vlog(data,pr)+vlog(1-data,1-pr)
	lt<--2*sum(lk[which(!is.na(data))])
	cat("Iteration: ",formatC(itel,digits=6,width=6),
	    " Deviance: ",formatC(ls,digits=6,width=12,format="f"),
	    " ==>",formatC(lt,digits=6,width=12,format="f"))
	if (correct)
	    cat(" Correct: ",formatC(sum(diag(tb))/sum(tb),digits=6,width=8,format="f"))
	if (extreme)
		cat(" Extreme: ",formatC(pp,digits=6,width=8 ,format="f"))
	cat("\n")
	if (((ls-lt)<eps) || (itel==imax)) break()
	itel<-itel+1
	ls<-lt
	z<-ifelse(is.na(data),aa,aa-4*(pr-data))
	}
radius<-mean(rowSums(x^2))
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]))
pdf(file=paste(paste(name,"_row",sep=""),"pdf",sep="."),encoding="MacRoman")
plot(x,type="n",main=paste("Ideal point plot for",name), 
	xlab=paste("dimension",1),ylab=paste("dimension",2),xlim=offset*xx,ylim=offset*yy)
text(x,rownames,cex=.5,col="red")
dev.off()
pdf(file=paste(paste(name,"_col",sep=""),"pdf",sep="."),encoding="MacRoman")
plot(x,type="n",main=paste("Roll call plot for",name), 
	xlab=paste("dimension",1),ylab=paste("dimension",2),xlim=offset*xx,ylim=offset*yy)
text(x,rowlabs,cex=.5,col="green")
for (i in 1:dim(y)[2]) {
    intercept<--a[i]/y[2,i]; slope<--y[1,i]/y[2,i]
	abline(intercept,slope)
	u<-slope*intercept; v<-1+(slope^2); w<-(radius*v)-(intercept^2)
	if (w>0) { 
	    x1<--(u+sqrt(w))/v; x2<--(u-sqrt(w))/v;
		y1<-intercept+(slope*x1); y2<-intercept+(slope*x2)
		text(x1,y1,as.character(i),col="red")
		text(x2,y2,as.character(i),col="red")
		}
	}
dev.off()
close(outfile)
if (save) list(intercepts=a,rowpoints=x,columnpoints=y,probabilities=pr)
}