#
# these are functions for multidimensional scaling. they handle weights, 
# metric and non-metric, constraints on the configaration,
# and various types of individual differences models. 
#
# version 0.5.01 2005-10-08
# version 0.6.00 2006-02-22 (unified two versions)
# version 0.7.00 2007-05-11 (added spherical constraints)
#

require("pava")
require("mdsutils")


smacofCons<-function(
	diss,
	constraints,
	p,
	wgths=initWeights(diss),
	metric=TRUE,
	ties="primary",
	extra=FALSE,	
	verbose=FALSE,
	modulus=1,
	itmax=100,
	eps=1e-6) {
n<-attr(diss,"Size"); m<-length(diss); dhat<-normDiss(diss,wgths)
w<-vmat(wgths); v<-myGenInv(w); itel<-1;
x<-constraints(matrix(rnorm(n*p),n,p),w,extra)
d<-dist(x); lb<-sum(wgths*d*dhat)/sum(wgths*d^2); x<-lb*x; d<-lb*d
sold<-sum(wgths*(dhat-d)^2)
repeat {
	b<-bmat(dhat,wgths,d); y<-v%*%b%*%x
	y<-constraints(y,w,extra); e<-dist(y)
	ssma<-sum(wgths*(dhat-e)^2)
	if (!metric) {
	    if ((itel%%modulus) == 0) {
			if (ties=="primary") daux<-monregP(diss,e,wgths)
			if (ties=="secondary") daux<-monregS(diss,e,wgths)
			if (ties=="tertiary") daux<-monregT(diss,e,wgths)
			dhat<-normDiss(daux,wgths)
			}
		}
	snon<-sum(wgths*(dhat-e)^2)
	if (verbose) cat("Iteration: ",formatC(itel,width=3, format="d")," Stress: ",
		formatC(c(sold,ssma,snon),digits=8,width=12,format="f"),"\n")
	if (((sold-snon)<eps) || (itel == itmax)) break()
	x<-y; d<-e; sold<-snon; itel<-itel+1	
	}
list(dhat=dhat,d=e,x=y)
}

uniqueC<-function(x,w,extra) {
n<-dim(x)[1]; p<-dim(x)[2]-n
return(cbind(x[,1:p],diag(diag(w%*%x[,p+(1:n)])/diag(w))))
}

linearC<-function(x,w,extra) {
return(extra%*%solve(crossprod(extra,w%*%extra),crossprod(extra,w%*%x)))
}

diagC<-function(x,w,extra) {
return(extra%*%diag(colSums(extra*(w%*%x))/colSums(extra*(w%*%extra))))
}