# X is a list of matrices, Z is a list of matrices, y is a list of vectors

twoStepLS<-function(x,z,y,w="ident") {
br<-firstStep(x,y)
gm<-secondStep(z,br$b,w)
}

firstStep<-function(x,y) {
b<-list(); r<-list(); m<-length(x)
for (k in 1:m) {
	lsf<-lsfit(x[[k]],y[[k]])
	b<-appendList(lsf$coef); r<-appendList(lsf$residuals)
	}
return(list(b=b,r=r))
}

secondStep<-function(z,b,w) {
c<-matrix(0,p,p); d<-rep(0,p); m<-length(z); p<-dim(z[[1]])[2]
for (k in 1:m) {
    if (length(w)>1) ww<-w[[k]]
    	else if (w=="ident") ww<-diag(dim(z[[k]])[1])
    		else ww<-crossprod(x[[k]])
	c<-c+crossprod(z[[k]],ww%*%z[[k]])
	d<-d+crossprod(z[[k]],ww%*%b[[k]])
	}
return(solve(c,d))
}

suffStatML<-function(x,y) {
c<-list(); d<-list()
}

balancedML<-function(x,y,z,oneRes=TRUE) {
m<-length(z); p<-dim(z[[1]])[2]
b<-firstStep(x,y)$b
gm<-secondStep(z,b,w="ident"); w<-matrix(0,p,p)
for (k in 1:m) {
	r<-b[[k]]-z[[k]]%*%gm
	w<-w+outer(r,r)
	}
w<-w/m
gm<-secondStep(z,b,w=w)
}

appendList<-function(x,a) {
return(c(x,list(a)))
}

mkInteract<-function(x,z) {
}

mkMixed<-function(x,z) {
}

randomIntercept<-function(x,z) {
}

useEMforML<-function() {}

useALSforML<-function() {}

matUBTwo<-function(a,b) {
d<-chol(a+b); c<-solve(d); z<-crossprod(c,a%*%c); e<-eigen(z)
ep<-pmax(e$values,1-e$values); k<-e$vectors
return(crossprod(sqrt(ep)*crossprod(k,d)))
}

matUB<-function(x) {
s<-x[[1]]
for (i in 2:length(x)) s<-matUBTwo(s,x[[i]])
return(s)
}

matFromUP<-function(a,b,s=a+b) {
sa<-s-a; sb<-s-b; 
d<-chol(sa+sb); c<-solve(d); z<-crossprod(c,sa%*%c); e<-eigen(z)
ep<-pmin(e$values,1-e$values); k<-e$vectors
return(s-crossprod(sqrt(ep)*crossprod(k,d)))
}