#
#   logReg package
#   Copyright (C) 2006  Jan de Leeuw <deleeuw@stat.ucla.edu>
#   UCLA Department of Statistics, Box 951554, Los Angeles, CA 90095-1554
#   
#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#   
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
###################################################################
#
# version 1.0.0    2006-03-27  first official release
#

logReg<-function(tab,x,b,vv="newt",third=FALSE,relax=FALSE,safe=FALSE,itmax=1000,eps=1e-6,ibnd=3,verbose=FALSE) {
nn<-rowSums(tab); n1<-tab[,1]; itel<-1; x<-as.matrix(x); m<-dim(x)[2]
vf<-crossprod(x,nn*x)/4; vs<-((eigen(vf)$values)[1])*diag(m)
if (third) kk<-bounds(nn,x)[ibnd]
repeat {
	xb<-as.vector(x%*%b); l<-1/(1+exp(-xb)); jj<-((2*l)-1)/(2*xb)
	u<-nn*l-n1; g<-colSums(u*x); mg<-max(abs(g))
	f<--(sum(n1*xb)+sum(nn*log(1-l)))
	cat("Iteration: ",formatC(itel,digits=6,width=6),
		"Function: ",formatC(f,digits=6,width=12,format="f"),
		"Maxgrad:  ",formatC(mg,digits=6,width=12,format="f"),
		"\n")
	if ((mg < eps) || (itel == itmax)) break
	vn<-crossprod(x,(nn*l*(1-l))*x); vj<-crossprod(x,(nn*jj)*x)
	if (vv=="newt") v<-vn
	if (vv=="jajo") v<-vj
	if (vv=="fixt") v<-vf
	if (vv=="scal") v<-vs
	if (!third) d<--solve(v,g)
		else d<-cubic(g,vn,kk,verbose=verbose)
	mfac<--sum(g*d)/sum(d*(vj%*%d)); step<-1
	if (safe)  step<-mfac
    if (relax) step<-2*step
    b<-b+step*d
	itel<-itel+1
	}
e<-max(1-step*eigen(solve(v,vn))$values)
return(list(b=b,g=g,v=v,l=l,e=e))
}

cubic<-function(b,a,k,off=1,eps=1e-6,itmax=100,verbose=FALSE) {
e<-eigen(a); lb<-e$values; ev<-e$vectors; bb<-colSums(b*ev); k2<-k/2; l2<-2/k
xold<-(l2*max(0,max(-lb)))+off; itel=1
repeat {
	ff<-sum((bb^2)/((lb+k2*xold)^2))-(xold^2)
	gg<--2*(k2*sum((bb^2)/((lb+k2*xold)^3))+xold)
	xnew<-xold-(ff/gg)
	if (verbose) cat("Secular ** ",
		"Function: ",formatC(ff,digits=6,width=12,format="f"),
		"Gradient: ",formatC(gg,digits=6,width=12,format="f"),
		"Oldval:   ",formatC(xold,digits=6,width=12,format="f"),
		"Newval:   ",formatC(xnew,digits=6,width=12,format="f"),
		"\n")
	if ((abs(ff) < eps) || (itel == itmax)) break
	xold<-xnew; itel<-itel+1
	}
if (xnew < (l2*max(0,max(-lb)))) stop("trouble")
y<--bb/(lb+k2*xnew)
return(as.vector(ev%*%y))
}


bounds<-function(n,z,x=rep(1,dim(z)[2]),eps=1e-6,itmax=100) {
b1<-sum(n*rowSums(z*z)^1.5)
b2<-eigen(crossprod(z,(n*sqrt(rowSums(z*z)))*z))$values[1]
itel<-1; x<-x/sqrt(sum(x*x)); fold<-sum(n*abs(z%*%x)^3)
repeat {
	h<-colSums(as.vector(n*((z%*%x)^2)*sign(z%*%x))*z); xup<-h/sqrt(sum(h^2)); fup<-sum(n*abs(z%*%xup)^3)
	if (((fup - fold) < eps) || (itel == itmax)) break
	itel<-itel+1; x<-xup; fold<-fup
	}
return(c(b1,b2,fup)*sqrt(3)/18)
}