#
# elegant functions to minimize sstress
#
# version 0.0.1 2006-02-06 (first attempt, symmetric case)
# version 0.0.2 2006-02-07 (added rectangular case)
# version 0.1.0 2006-02-09 (fixed rectangular case)
# version 0.1.1 2006-02-09 (got rid of large matrix multiplication)
# version 0.1.2 2006-02-09 (improved stepsize)
# version 0.2.0 2006-02-09 (introduced linear restrictions in symmetric case)
# version 0.2.1 2006-02-20 (put in alternative cheaper bounds)
# version 0.2.2 2006-02-22 (no weights, no matrix storage)
#

require("mdsutils")

elegantSym<-function(
	delta,
	p=2,
	w=1,
	bnd="norm",
	init="none",
	verbose=TRUE,
	itmax=500,
	eps=1e-6) {
n<-dim(delta)[1]; itel<-1
if (!is.matrix(init)) x<-torgerson(sqrt(delta),p=p) else x<-init
if (!is.matrix(w)) sbd<-2*n
	else if (bnd == "norm") sbd<-4*max(rowSums(w))
		else sbd<-symBound(w)
d<-dist2Sym(x); r<-delta-d; sstress=sum(w*r*r)/2; psstress<-sstress
repeat {
    v<--w*r; diag(v)<--rowSums(v); x<-x+t(qr.solve(x,v))/sbd
	d<-dist2Sym(x); r<-delta-d; sstress=sum(w*r*r)/2
	if (verbose) {
		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))
}

elegantRect<-function(
	delta,
	p=2,
	w=1,
	bnd="norm",
	init="torg",
	verbose=TRUE,
	itmax=500,
	eps=1e-6) {
n<-dim(delta)[1]; m<-dim(delta)[2]; itel<-1
if (!is.matrix(w)) sbd<-n+m+2 
	else if (bnd == "norm") sbd<-4*max(c(rowSums(w),colSums(w)))
		else sbd<-rectBound(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<-delta-d
sstress=sum(w*r*r)/2; psstress<-sstress
repeat {
    v<--w*r; vr<--rowSums(v); vc<--colSums(v); c<-crossprod(x)+crossprod(y)
    u<-solve(c,t(x*vr)+crossprod(y,t(v))); v<-solve(c,t(y*vc)+crossprod(x,v))
    x<-x+t(u)/sbd; y<-y+t(v)/sbd
	d<-dist2Rect(x,y); r<-delta-d; rr<-rowSums(w*r); sstress=sum(w*r*r)/2
	if (verbose) {
		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))
}

elegantGen<-function() {}

elegantThreeSym<-function() {}

elegantLinSym<-function(
	delta,
	w=1,
	z="ident",
	p=2,
	init="torg",
	verbose=TRUE,
	itmax=500,
	eps=1e-6) {
n<-dim(delta)[1]; itel<-1; 
if (is.character(z)) z<-diag(n); zz<-qr(z); z<-qr.Q(zz)
if (!is.matrix(w)) sbd<-2*n
	else if (bnd == "norm") sbd<-4*max(rowSums(w))
		else sbd<-symBound(w)
if (!is.matrix(init)) x<-torgerson(sqrt(delta),p=p) else x<-init
a<-crossprod(z,x); x<-z%*%a; d<-dist2Sym(x); r<-delta-d
sstress=sum(w*r*r)/2; psstress<-sstress
repeat {
    v<--w*r; diag(v)<--rowSums(v); v<-crossprod(z,v)%*%z; a<-a+t(qr.solve(a,v))/sbd
	x<-z%*%a; d<-dist2Sym(x); r<-delta-d; sstress=sum(w*r*r)/2
	if (verbose) {
		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))
}

symBound<-function(w,eps=1e-6) {
n<-dim(w)[1]; x<-diag(n)/sqrt(n); z<-matrix(0,n,n); plbd<--Inf
repeat {
    d<-diag(x); c<-outer(d,d,"+")-2*x; z<--w*c
	diag(z)<--rowSums(z)
	lbd<-sum(z*x)
	x<-z/sqrt(sum(z^2))	
	if ((lbd-plbd)<eps) break()
	plbd<-lbd
	}
return(lbd)
}

rectBound<-function(w,eps=1e-6) {
n<-dim(w)[1]; m<-dim(w)[2]; nm<-n+m; 
nn<-1:n; no<-(n+1):nm; x<-diag(nm)/sqrt(nm); plbd<--Inf
repeat {
    z<-matrix(0,nm,nm)
    d<-diag(x); c<-outer(d,d,"+")-2*x
    z[nn,no]<--w*c[nn,no]; z[no,nn]<--t(w)*c[no,nn]
	diag(z)<--rowSums(z)
	lbd<-sum(z*x); print(lbd)
	x<-z/sqrt(sum(z^2))	
	if ((lbd-plbd)<eps) break()
	plbd<-lbd
	}
return(lbd)
}