#
# alscal package
# functions to minimize sstress
#
# version 0.0.1  only metric, one matrix
#
# version 0.1.0 	02-06-06 (added elegant algorithm)
# version 0.1.1 	02-06-06 (elegant algorithm with power iterations)
# version 0.1.2     02-07-06 (removed elegant to separate package)


require("mdsutils")
require("lowpoly")

alscalRect<-function(
	delta,
	w="none",
	p=2,
	init="svd",
	verbose=FALSE,
	ithist=TRUE,
	itmax=100,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
if (!is.matrix(w)) w<-matrix(1,n,m)
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))
}

alscalSym<-function(
	delta,
	w=matrix(1,dim(delta)[1],dim(delta)[1]),
	p=2,
	init="torg",
	verbose=FALSE,
	ithist=TRUE,
	itmax=100,
	eps=1e-6) {
n<-dim(delta)[1]; itel<-1
nn<-1:n; pp<-1:p; wr<-rowSums(w)
if (is.matrix(init)) x<-init
else {
	e<--0.5*(delta-outer(rowSums(delta)/n,colSums(delta)/n,"+")+(sum(delta)/(n*n)))
	z<-svd(e,nu=p,nv=0); x<-t(sqrt(z$d[1:p])*t(z$u))
	}
d<-dist2Sym(x); r<-delta-d; rr<-rowSums(w*r)
sstress=sum(w*r*r)/2; psstress<-sstress
repeat {
    for (i in nn) for (s in pp) {
	    u<-x[i,s]-x[,s]; ww<-w[i,]
		a4<-wr[i]; a3<-2*sum(ww*u)
		a2<-(4*sum(ww*u*u))-(2*rr[i])
		a1<--4*sum(ww*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)
		d[,i]<-d[,i]+(2*th*u)+(th*th); d[i,i]<-d[i,i]-(2*th*th)
		r<-delta-d; rr<-rowSums(w*r); sstress=sum(w*r*r)/2
		if (verbose) print(sstress)
		}
	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,d=d,sstress=sstress))
}