#
#   voronoi 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-04-25,		first release
# version 1.0.1, 2006-04-26,		bug fix
# version 1.0.2, 2006-04-26, 		another bug fix
# version 1.1.0, 2006-04-26,        another parametrization
#

makeEdges<-function(a,b,fit="ip",verbose=FALSE) {
if (fit=="ip") c<--log(b) else c<-(-log(b)+rowSums(a^2))/2
n<-length(b); lns<-matrix(0,0,8)
for (i in 1:(n-1)) {
	for(j in (i+1):n) {
	    dd<-a[i,]-a[j,]; dc<-c[i]-c[j]; ss<-sum(dd^2)
	    if (is.nul(ss)) next() 
	    ee<-dc*dd/ss; ff<-c(-dd[2],dd[1])
		xlw<--Inf; xup<-Inf
		for (k in (1:n)[-c(i,j)]) {
		    dd<-a[i,]-a[k,]; dc<-c[i]-c[k]
		    mum<-sum(dd*ff); mom<-dc-sum(dd*ee)
		    if (is.nul(mum) & (mom > 0)) {
		    	xlw<-Inf; xup<--Inf
		    	}
		    if (mum>0) xlw<-max(xlw,mom/mum)
		    if (mum<0) xup<-min(xup,mom/mum)
		    if (verbose) {
				cat(formatC(i,digits=3,width=3),
					formatC(j,digits=3,width=3),
					formatC(k,digits=3,width=3),
					"mum ",formatC(mum,digits=4,width=8,format="f"),
					"mom ",formatC(mom,digits=4,width=8,format="f"),
					"xlw ",formatC(xlw,digits=4,width=8,format="f"),
					"xup ",formatC(xup,digits=4,width=8,format="f"),
					"\n")
					}
			}
		if (xlw<xup) lns<-rbind(lns,c(i,j,ee,ff,xlw,xup))
		}
	}
return(lns)
}

drawEdges<-function(a,b,fit="ip",far=1000,verbose=FALSE) {
lns<-makeEdges(a,b,fit,verbose=verbose)
p<-dim(lns)[1]; n<-dim(a)[1]
plot(a,type="n"); text(a,as.character(1:n))
for (i in 1:p) {
	ee<-lns[i,3:4]; ff<-lns[i,5:6]
	xlw<-lns[i,7]; if (xlw == -Inf) xlw<--far
	plw<-ee+xlw*ff
	xup<-lns[i,8]; if (xup == Inf) xup<-far
	pup<-ee+xup*ff
	lines(rbind(plw,pup),col="RED")
	}
}
	
is.nul<-function(x) {
return(abs(x)<1e-10)
}