#
#	Fest.S
#
#	S function empty.space()
#	Computes estimates of the empty space function
#
#	$Revision: 4.5 $	$Date: 2001/11/06 06:59:55 $
#
"Fest" <- 	
"empty.space" <-
function(X, eps = NULL, r=NULL, breaks=NULL) {
#
#	pp:		point pattern (an object of class 'ppp')
#	eps:		raster grid mesh size for distance transform
#				(unless specified by pp$window)
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
# First discretise
	dwin <- as.mask(X$window, eps)
        dX <- ppp(X$x, X$y, window=dwin)
#        
# histogram breakpoints 
#
        breaks <- handle.r.b.args(r, breaks, dwin, eps)
#
#  compute distances and censoring distances
	if(X$window$type == "rectangle") {
                # original data were in a rectangle
                # output of exactdt() is sufficient
		e <- exactdt(dX)
		dist <- e$d
		bdry <- e$b
	} else {
                # window is irregular..
          
                # Distance transform & boundary distance for all pixels
		e <- exactdt(dX)
		b <- bdist.pixels(dX$window, coords=F)
                # select only those pixels inside mask
		dist <- e$d[dwin$m]
		bdry <- b[dwin$m]
	} 
# censoring indicators
	d <- (dist <= bdry)
#  observed distances
	o <- pmin(dist, bdry)
#        
#
# calculate Kaplan-Meier and border corrected estimates
	result <- km.rs(o, bdry, d, breaks)

# also calculate UNCORRECTED e.d.f. !!!! use with care
        hh <- hist(dist,breaks=breaks$val,plot=F)$counts
        edf <- cumsum(hh)/sum(hh)
        result$raw <- edf

# append theoretical value for Poisson
        lambda <- X$n/area.owin(X$window)
        result$theo <- 1 - exp( - lambda * pi * result$r^2)

# neaten up and return        
        result$breaks <- NULL
        
	return(data.frame(result))
}

	
.First.lib <- function(lib, pkg) {
    library.dynam("spatstat", pkg, lib)
    cat("spatstat\nversion: 1.0-1\n")
}

#
#	Gest.S
#
#	Compute estimates of nearest neighbour distance distribution function G
#
#	$Revision: 4.3 $	$Date: 2001/11/06 07:00:19 $
#
################################################################################
#
"Gest" <-
"nearest.neighbour" <-
function(X, r=NULL, breaks=NULL, ...) {
#	X		point pattern (of class ppp)
#				(unless specified by X$window)
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
	verifyclass(X, "ppp")
        
#  compute nearest neighbour distances
	nnd <- nndist(X$x, X$y)
		
#  determine breakpoints for r values
        breaks <- handle.r.b.args(r, breaks, X$window)

#  UNCORRECTED e.d.f. of nearest neighbour distances: use with care
        hh <- hist(nnd,breaks=breaks$val,plot=F)$counts
        edf <- cumsum(hh)/sum(hh)
	
#  distance to boundary
        bdry <- bdist.points(X)

#  observations
	o <- pmin(nnd,bdry)
#  censoring indicators
	d <- (nnd <= bdry)
#
# calculate Kaplan-Meier and border correction (Reduced Sample) estimators
	result <- km.rs(o, bdry, d, breaks)
#        
# append uncorrected e.d.f.        
        result$raw <- edf
# append theoretical value for Poisson
        lambda <- X$n/area.owin(X$window)
        result$theo <- 1 - exp( - lambda * pi * result$r^2)

# neaten up and return        
        result$breaks <- NULL

	return(data.frame(result))
}	

nndist <- function(x, y) {
	#  computes the vector of nearest-neighbour distances 
	#  for the pattern of points (x[i],y[i])
	#	
#  matrix of squared distances between all pairs of points
	sq <- function(a, b) { (a-b)^2 }
	squd <-  outer(x, x, sq) + outer(y, y, sq)
#  reset diagonal to a large value so it is excluded from minimum
	diag(squd) <- Inf
#  nearest neighbour distances
	nnd <- sqrt(apply(squd,1,min))
	nnd
}
#	Gmulti.S
#
#	Compute estimates of nearest neighbour distance distribution functions
#	for multitype point patterns
#
#	S functions:	
#		Gcross                G_{ij}
#		Gdot		      G_{i\bullet}
#		Gmulti	              (generic)
#
#	$Revision: 4.4 $	$Date: 2001/11/06 10:55:56 $
#
################################################################################

"Gcross" <-		
function(X, i=1, j=2, r=NULL, breaks=NULL, ...) {
#	computes G_{ij} estimates
#
#	X		marked point pattern (of class 'ppp')
#	i,j		the two mark values to be compared
#  
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
	X <- as.ppp(X)
	if(!is.marked(X))
		stop("point pattern has no \'marks\'")
	window <- X$window
# 
	I <- (X$marks == i)
	if(sum(I) == 0) stop("No points are of type i")
	J <- (X$marks == j)
	if(sum(J) == 0) stop("No points are of type j")
#
	Gmulti(X, I, J, r, breaks)
}	

"Gdot" <- 	
function(X, i=1, r=NULL, breaks=NULL, ...) {
#  Computes estimate of 
#      G_{i\bullet}(t) = 
#  P( a further point of pattern in B(0,t)| a type i point at 0 )
#	
#	X		marked point pattern (of class ppp)
#  
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
	X <- as.ppp(X)

	if(!is.marked(X))
		stop("point pattern has no \'marks\'")
#
	I <- (X$marks == i)
	if(sum(I) == 0) stop("No points are of type i")
	J <- rep(T, X$n)	# i.e. all points
# 
	Gmulti(X, I, J, r, breaks)
}	

	
##########

"Gmulti" <- 	
function(X, I, J, r=NULL, breaks=NULL, ...) {
#
#  engine for computing the estimate of G_{ij} or G_{i\bullet}
#  depending on selection of I, J
#  
#	X		marked point pattern (of class ppp)
#	
#	I,J		logical vectors of length equal to the number of points
#			and identifying the two subsets of points to be
#			compared.
#  
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
	verifyclass(X, "ppp")
	if(!is.marked(X))
		stop("point pattern has no \'marks\'")
	window <- X$window
        brks <- handle.r.b.args(r, breaks, window)$val        
# Extract points of group I and J
	if(!is.logical(I) || !is.logical(J))
		stop("I and J must be logical vectors")
	if(length(I) != X$n || length(J) != X$n)
		stop("length of I or J does not equal the number of points")
	if(sum(I) == 0) stop("No points satisfy condition I")
	xI <- X$x[I]
	yI <- X$y[I]
	if(sum(J) == 0) stop("No points satisfy condition J")
	xJ <- X$x[J]
	yJ <- X$y[J]
#  compute squared distance from each type i point to
#  the nearest other point of any type
	xIm <- matrix(xI, nrow=length(xI), ncol=length(xJ))
	yIm <- matrix(yI, nrow=length(yI), ncol=length(yJ))
	xJm <- matrix(xJ, nrow=length(xJ), ncol=length(xI))
	yJm <- matrix(yJ, nrow=length(yJ), ncol=length(yI))
	squd <-  (xIm - t(xJm))^2 + (yIm - t(yJm))^2
# identify which pairs in 'squd' correspond to identical points
	dag <- (diag(1, X$n) != 0)
	dag <- dag[I,J]
#  reset this 'diagonal' to a large value
	oo <- (diff(window$xrange) + diff(window$yrange) )^2
	squd[dag] <- oo
#  "type I to type J" nearest neighbour distances
	nnd <- sqrt(apply(squd,1,min))
#  distance to boundary from each type i point
        bdry <- bdist.points(X[I, ])
#  observations
	o <- pmin(nnd,bdry)
#  censoring indicators
	d <- (nnd <= bdry)
#
# calculate
	result <- km.rs(o, bdry, d, brks)
        result$breaks <- NULL

#  UNCORRECTED e.d.f. of I-to-J nearest neighbour distances: use with care
        hh <- hist(nnd,breaks=brks,plot=F)$counts
        result$raw <- cumsum(hh)/sum(hh)

# theoretical value for marked Poisson processes
        lamJ <- sum(J)/area.owin(window)
        result$theo <- 1 - exp( - lamJ * pi * result$r^2)
        
	return(data.frame(result))
}	


#	Jest.S
#
#	Usual invocation to compute J function
#	if F and G are not required 
#
#	$Revision: 4.4 $	$Date: 2001/11/13 01:25:49 $
#
#
#
"Jest" <-
function(X, eps=NULL, r=NULL, breaks=NULL) {
        X <- as.ppp(X)
        brks <-  handle.r.b.args(r, breaks, X$window)$val
	FF <- Fest(X, eps, breaks=brks)
	G <- Gest(X, breaks=brks)
        ratio <- function(a, b, c) {
          result <- a/b
          result[ b == 0 ] <- c
          result
        }
	Jrs <- ratio(1-G$rs, 1-FF$rs, NA)
	Jkm <- ratio(1-G$km, 1-FF$km, NA)
	Jun <- ratio(1-G$raw, 1-FF$raw, NA)
        theo <- rep(1, length(FF$r))

        rslt <- data.frame(J=Jkm,r=FF$r,km=Jkm,rs=Jrs,un=Jun,theo=theo)
        attr(rslt, "F") <- FF
        attr(rslt, "G") <- G
        rslt
}

#	Jmulti.S
#
#	Usual invocations to compute multitype J function(s)
#	if F and G are not required 
#
#	$Revision: 4.3 $	$Date: 2001/11/12 10:40:04 $
#
#
#
"Jcross" <-
function(X, i=1, j=2, eps=NULL, r=NULL, breaks=NULL) {
#
#       multitype J function J_{ij}(r)
#  
#	X:		point pattern (an object of class 'ppp')
#       i, j:           types for which J_{i,j}(r) is calculated  
#	eps:		raster grid mesh size for distance transform
#				(unless specified by X$window)
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
        X <- as.ppp(X)
        if(!is.marked(X))
          stop("point pattern has no \'marks\'")
        I <- (X$marks == i)
        J <- (X$marks == j)
        result <- Jmulti(X, I, J, eps, r, breaks)
	return(result)
}

"Jdot" <-
function(X, i=1, eps=NULL, r=NULL, breaks=NULL) {
#  
#    multitype J function J_{i\dot}(r)
#  
#	X:		point pattern (an object of class 'ppp')
#       i:              mark i for which we calculate J_{i\cdot}(r)  
#	eps:		raster grid mesh size for distance transform
#				(unless specified by X$window)
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#
        X <- as.ppp(X)
        if(!is.marked(X))
          stop("point pattern has no \'marks\'")
        I <- (X$marks == i)
        J <- rep(T, X$n)
        result <- Jmulti(X, I, J, eps, r, breaks)
	return(result)
}

"Jmulti" <- 	
function(X, I, J, eps=NULL, r=NULL, breaks=NULL) {
#  
#    multitype J function (generic engine)
#  
#	X		marked point pattern (of class ppp)
#	
#	I,J		logical vectors of length equal to the number of points
#			and identifying the two subsets of points to be
#			compared.
#  
#	eps:		raster grid mesh size for distance transform
#				(unless specified by X$window)
#  
#       r:              (optional) values of argument r  
#	breaks:		(optional) breakpoints for argument r
#  
#
        X <- as.ppp(X)
        brks <- handle.r.b.args(r, breaks, X$window, eps)$val
	FJ <- Fest(X[J], eps, breaks=brks)
	GIJ <- Gmulti(X, I, J, breaks=brks)
        ratio <- function(a, b, c) {
          result <- a/b
          result[ b == 0 ] <- c
          result
        }
	Jkm <- ratio(1-GIJ$km, 1-FJ$km, NA)
	Jrs <- ratio(1-GIJ$rs, 1-FJ$rs, NA)
	Jun <- ratio(1-GIJ$raw, 1-FJ$raw, NA)
        theo <- rep(1, length(FJ$r))
        
	result <- data.frame(J=Jkm,r=FJ$r,km=Jkm,rs=Jrs,un=Jun,theo=theo)
        attr(result, "G") <- GIJ
        attr(result, "F") <- FJ
        result
}
#
#	Kest.S		Estimation of K function
#
#	$Revision: 4.5 $	$Date: 2001/11/06 10:27:52 $
#
#	Border correction only.
#	Arbitrary window.
#       Fastest implementation, using 'hist' and 'reduced.sample'
#
#
# -------- functions ----------------------------------------
#	Kest()		compute estimate of K by border correction
#
#	pairdist()	compute matrix of distances between 	
#			  each pair of data points
#
#       Kount()         internal routine
#
# -------- standard arguments ------------------------------	
#	pp		point pattern (of class 'ppp')
#
#	r		distance values at which to compute K	
#
# -------- standard outputs ------------------------------	
#	border:		K function estimated by border method
#			using standard formula (denominator = count of points)
#
#       bord.modif:	K function estimated by border method
#			using modified formula 
#			(denominator = area of eroded window
#
#	r:		same as input
#
# -------- example -----------------------------------------
#
#	x <- runif(50)
#	y <- runif(50)
#	X <- ppp(x,y,c(0,1),c(0,1))
# either:
#	r <- seq(0,0.25,length=50)
#	k <- Kest(X, r)$border
#	plot(r, k, type="l")
#	lines(r, pi * r^2, lty=2)
# or:
#       Kdat <- Kest(X)
#       plot(border ~ r, type="l", data=Kdat)
# ------------------------------------------------------------------------

"pairdist"<-
function(x, y)
{
	if(length(x) != length(y))
		stop("lengths of x and y do not match")
	n <- length(x)
	xx <- matrix(rep(x, n), nrow = n)
	yy <- matrix(rep(y, n), nrow = n)
	d <- sqrt((xx - t(xx))^2 + (yy - t(yy))^2)
	return(d)
}
	
"Kest"<-
function(X, r=NULL, breaks=NULL, slow=F, both=F, ...)
{
	verifyclass(X, "ppp")

	npoints <- X$n
	area <- area.owin(X$window)
	lambda <- npoints/area
	lambda2 <- (npoints * (npoints - 1))/(area^2)

        breaks <- handle.r.b.args(r, breaks, X$window)
        r <- breaks$r
        
        # pairwise distance
	d <- pairdist(X$x, X$y)
        diag(d) <- Inf

        # distances to boundary
	b <- bdist.points(X)

        # count 
        RS <- Kount(d, b, breaks, slow)
        
        K.b1 <- RS$numerator/(lambda * RS$denom.count)
        
        if(!both)
          return(data.frame(border=K.b1, r=r, theo= pi * r^2))
        
        denom.area <- eroded.areas(X$window, r)
        K.b2 <- RS$numerator/(lambda2 * denom.area)

        return(data.frame(border=K.b1, bord.modif=K.b2, r=r,
                          theo=pi * r^2))
}
	
Kount <- function(d, b, breaks, slow=F) {
  #
  # "internal" routine to compute border-correction estimate of K or Kij
  #
  # d : matrix of pairwise distances
  #                  (to exclude diagonal entries, set diag(d) = Inf)
  # b : column vector of distances to window boundary
  # breaks : breakpts object
  #

  if(slow) { ########## slow ##############
          
       r <- breaks$r
       
       nr <- length(r)
       numerator <- numeric(nr)
       denom.count <- numeric(nr)

       for(i in 1:nr) {
         close <- (d <= r[i])
         nclose <- apply(close, 1, sum) # assumes diag(d) set to Inf
         bok <- (b > r[i])
         numerator[i] <- sum(nclose[bok])
         denom.count[i] <- sum(bok)
       }
	
  } else { ############# fast ####################

        # determine which distances d_{ij} were observed without censoring
        bb <- matrix(b, nrow=nrow(d), ncol=ncol(d))
        uncen <- (d <= bb)
        #
        # histogram of noncensored distances
        nco <- hist(d[uncen], breaks=breaks$val,plot=F,probability=F)$counts
        # histogram of censoring times for noncensored distances
        ncc <- hist(bb[uncen], breaks=breaks$val,plot=F,probability=F)$counts
        # histogram of censoring times (yes, this is a different total size)
        cen <- hist(b, breaks=breaks$val,plot=F,probability=F)$counts
        # go
        RS <- reduced.sample(nco, cen, ncc, show=T)
        # extract results
        numerator <- RS$numerator
        denom.count <- RS$denominator
        # check
        if(length(numerator) != breaks$ncells)
          stop("internal error: length(numerator) != breaks$ncells")
        if(length(denom.count) != breaks$ncells)
          stop("internal error: length(denom.count) != breaks$ncells")
             
  }
  
  return(list(numerator=numerator, denom.count=denom.count))
}
#
#	Kinhom.S	Estimation of K function for inhomogeneous patterns
#
#	$Revision: 1.1 $	$Date: 2001/11/21 09:07:22 $
#
#	Kinhom()	compute estimate of K_inhom
#
#       Currently uses border method and slow code...
#                       
#       Reference:
#            Non- and semiparametric estimation of interaction
#	     in inhomogeneous point patterns
#            A.Baddeley, J.Moller, R.Waagepetersen
#            Statistica Neerlandica 54 (2000) 329--350.
#
"Kinhom"<-
  function (X, lambda, r = NULL, breaks = NULL) {
    verifyclass(X, "ppp")
    breaks <- handle.r.b.args(r, breaks, X$window)
    r <- breaks$r
    if (!is.vector(lambda))
        stop("The argument \'lambda\' should be a vector")
    if (length(lambda) != X$n) 
        stop("The length of the vector \'lambda\' should equal the number of data points")
    d <- pairdist(X$x, X$y)
    diag(d) <- Inf
    b <- bdist.points(X)
    w <- 1/lambda
    RS <- Kwtsum(d, b, w, w, breaks)
    K <- RS$numerator/RS$denominator
    return(data.frame(K = K, r = r, theo = pi*r**2, border=K))
}


	
Kwtsum <- function(d, b, wi, wj, breaks) {
  #
  # "internal" routine to compute border-correction estimate of Kinhom
  #
  # d : matrix of pairwise distances
  #                  (to exclude diagonal entries, set diag(d) = Inf)
  # b : column vector of distances to window boundary
  # wi : vector of weights on rows of d (estimate of 1/lambda_i)
  # wj : vector of weights on columns of d (estimate of 1/lambda_j)
  # 
  # breaks : breakpts object
  #

  if(length(wi) != nrow(d))
    stop("length of \'wi\' doesn\'t match number of rows of \'d\'")
  
  if(length(wj) != ncol(d))
    stop("length of \'wj\' doesn\'t match number of columns of \'d\'")
  
  r <- breaks$r
  nr <- length(r)
  numerator <- numeric(nr)
  denominator <- numeric(nr)

  wiwj <- outer(wi, wj, "*")
       
  for(k in 1:nr) {
    close <- (d <= r[k]) # assuming d[i,j] = Inf when X[i] == X[j]
    numi <- matrowsum(close * wiwj)
    bok <- (b > r[k]) 
    numerator[k] <- sum(numi[bok])
    denominator[k] <- sum(wi[bok])
  }
	
  return(list(RS=numerator/denominator,
              numerator=numerator,
              denominator=denominator))
}
#
#	Kmulti.S		
#
#	Compute estimates of cross-type K functions
#	for multitype point patterns
#
#	$Revision: 4.4 $	$Date: 2001/11/06 11:05:09 $
#
#	Border correction only.
#	Arbitrary window.
#
# -------- functions ----------------------------------------
#	Kcross()	cross-type K function K_{ij}
#                       between types i and j
#
#	Kdot()          K_{i\bullet}
#                       between type i and all points regardless of type
#
#       Kmulti()        (generic)
#
#	crossdist()	compute matrix of distances between 	
#			  each pair of data points 
#			  in two separate lists of points
#
# -------- standard arguments ------------------------------	
#	X		point pattern (of class 'ppp')
#				including 'marks' vector
#	r		distance values at which to compute K	
#
# -------- standard outputs ------------------------------	
#	border:		cross K function estimated by border method
#			using standard formula (denominator = count of points)
#
#       bord.modif:	cross K function estimated by border method
#			using modified formula 
#			(denominator = area of eroded window
#
#	r:		same as input
#
#
# -------- example -----------------------------------------
#
#	x <- runif(50)
#	y <- runif(50)
#	m <- c(rep(1,25),rep(2,25))
#	X <- ppp(x,y,c(0,1),c(0,1),marks=m)
# either:
#	r <- seq(0,0.5,length=50)
#	k <- Kcross(X, 1, 2, r)$border
#	plot(r, k, type="l")
#	lines(r, pi * r^2, lty=2)
# or:
#       Kdat <- Kcross(X, 1, 2)
#       plot(border ~ r, type="l", data=Kdat)
#
# ------------------------------------------------------------------------

"crossdist"<-
function(x1, y1, x2, y2)
{
        # returns matrix[i,j] = distance from (x1[i],y1[i]) to (x2[j],y2[j])
	if(length(x1) != length(y1))
		stop("lengths of x1 and y1 do not match")
	if(length(x2) != length(y2))
		stop("lengths of x2 and y2 do not match")
	n1 <- length(x1)
	n2 <- length(x2)
	X1 <- matrix(rep(x1, n2), ncol = n2)
	Y1 <- matrix(rep(y1, n2), ncol = n2)
	X2 <- matrix(rep(x2, n1), ncol = n1)
	Y2 <- matrix(rep(y2, n1), ncol = n1)
	d <- sqrt((X1 - t(X2))^2 + (Y1 - t(Y2))^2)
	return(d)
}

"Kcross" <- 
function(X, i=1, j=2, r=NULL, breaks=NULL, ...) {
	verifyclass(X, "ppp")
	if(!is.marked(X))
		stop("point pattern has no marks (no component 'marks')")
	
	I <- (X$marks == i)
	J <- (X$marks == j)
	
	if(!any(I)) stop(paste("No points have mark i =", i))
	if(!any(J)) stop(paste("No points have mark j =", j))
	
	Kmulti(X, I, J, r, breaks)
}

"Kdot" <- 
function(X, i=1, r=NULL, breaks=NULL, ...) {
	verifyclass(X, "ppp")
	if(!is.marked(X))
		stop("point pattern has no marks (no component 'marks')")
	
	I <- (X$marks == i)
	J <- rep(T, X$n)  # i.e. all points
	
	if(!any(I)) stop(paste("No points have mark i =", i))
	
	Kmulti(X, I, J, r, breaks)
}


"Kmulti"<-
function(X, I, J, r=NULL, breaks=NULL, ...)
{
	verifyclass(X, "ppp")

	npoints <- X$n
	x <- X$x
	y <- X$y
	area <- area.owin(X$window)

        breaks <- handle.r.b.args(r, breaks, X$window)
        r <- breaks$r
        
	if(!is.logical(I) || !is.logical(J))
		stop("I and J must be logical vectors")
	if(length(I) != npoints || length(J) != npoints)
	     stop("The length of I and J must equal \
 the number of points in the pattern")
	
	if(!any(I)) stop("no points satisfy I")
	if(!any(J)) stop("no points satisfy J")
		
	nI <- sum(I)
	nJ <- sum(J)
	lambdaI <- nI/area
	lambdaJ <- nJ/area

# interpoint distances		
	d <- crossdist(x[I], y[I], x[J], y[J])
# distances to boundary	
	b <- (bdist.points(X))[I]
        
# interpoint distances that refer to identical points
# are excluded from consideration by setting them to infinity        
        common <- I & J
        if(any(common)) {
          Irow <- cumsum(I)
          Jcol <- cumsum(J)
          icommon <- (1:npoints)[common]
          for(i in icommon)
            d[Irow[i], Jcol[i]] <- Inf
        }
          
# border correction estimator          
	RS <- Kount(d, b, breaks, slow=F)  
        
        K.b1 <- RS$numerator/(lambdaJ * RS$denom.count)

# variant 
        denom.area <- eroded.areas(X$window, r)
        K.b2 <- RS$numerator/(denom.area * nI * nJ)

# theoretical value for marked Poisson
        theo <- pi * r^2
        
        return(data.frame(border=K.b1, bord.modif = K.b2, r=r, theo=theo))
          
}
#
#	affine.S
#
#	$Revision: 1.1 $	$Date: 2001/11/21 09:07:22 $
#

affinexy <- function(X, mat=diag(c(1,1)), vec=c(0,0)) {
  # Y = M X + V
  ans <- mat %*% rbind(X$x, X$y) + matrix(vec, nrow=2, ncol=length(X$x))
  list(x = ans[1,],
       y = ans[2,])
}

"affine.owin" <- function(W,  mat=diag(c(1,1)), vec=c(0,0)) {
  verifyclass(W, "owin")
  # Inspect the determinant
  detmat <- det(mat)
  if(abs(detmat) < .Machine$double.eps)
    stop("Matrix of linear transformation is singular")
  #
  switch(W$type,
         rectangle={
           # convert rectangle to polygon
           P <- owin(W$xrange, W$yrange, poly=
                     list(x=W$xrange[c(1,2,2,1)],
                          y=W$yrange[c(1,1,2,2)]))
           # call polygonal case
           return(affine.owin(P, mat, vec))
         },
         polygonal={
           # First transform the polygonal boundaries
           bdry <- lapply(W$bdry, affinexy, mat=mat, vec=vec)
           # If determinant < 0, traverse polygons in reverse direction
           if(detmat < 0)
             bdry <- lapply(bdry, reverse.xypolygon)
           # Compute bounding box of new polygons
           xr <- range(unlist(lapply(bdry, function(a) a$x)))
           yr <- range(unlist(lapply(bdry, function(a) a$y)))
           # wrap up
           return(owin(xr, yr, poly=bdry))
         },
         mask={
           stop("Sorry, \'affine.owin\' is not yet implemented for masks")
         },
         stop("Unrecognised window type")
         )
}

"affine.ppp" <- function(X, mat=diag(c(1,1)), vec=c(0,0)) {
  verifyclass(X, "ppp")
  r <- affinexy(X, mat, vec)
  w <- affine.owin(X$window, mat, vec)
  return(ppp(r$x, r$y, window=w, marks=X$marks))
}


"affine" <- function(X, ...) {
  UseMethod("affine")
}

### ---------------------- shift ----------------------------------

"shift" <- function(X, ...) {
  UseMethod("shift")
}

shiftxy <- function(X, vec=c(0,0)) {
  list(x = X$x + vec[1],
       y = X$y + vec[2])
}

"shift.owin" <- function(W,  vec=c(0,0)) {
  verifyclass(W, "owin")
  # Shift the bounding box
  W$xrange <- W$xrange + vec[1]
  W$yrange <- W$yrange + vec[2]
  switch(W$type,
         rectangle={
         },
         polygonal={
           # Shift the polygonal boundaries
           W$bdry <- lapply(W$bdry, shiftxy, vec=vec)
         },
         mask={
           # Shift the pixel coordinates
           W$xcol <- W$xcol + vec[1]
           W$yrow <- W$yrow + vec[2]
           # That's all --- the mask entries are unchanged
         },
         stop("Unrecognised window type")
         )
  return(W)
}

"shift.ppp" <- function(X, vec=c(0,0)) {
  verifyclass(X, "ppp")
  r <- shiftxy(X, vec)
  w <- shift.owin(X$window, vec)
  return(ppp(r$x, r$y, window=w, marks=X$marks))
}


#
#
#   allstats.R
#
#   $Revision: 1.5 $   $Date: 2001/11/06 07:09:17 $
#
#
allstats <- function(pp,dataname=NULL,verb=F) {
#
# Function allstats --- to calculate the F, G, K, and J functions
# for an unmarked point pattern.
#
  verifyclass(pp,"ppp")
  if(is.marked(pp))
    stop("This function is applicable only to unmarked patterns.\n")

# get sensible r values
  brks <- handle.r.b.args(r=NULL, breaks=NULL, pp$window, eps=NULL)
  r    <- brks$r

# initialise  
  fns <- list()
  titles <- list()

# estimate F, G and J 
  if(verb) cat("Calculating F, G, J ...")
  Jout <- Jest(pp,eps=NULL,breaks=brks)
  if(verb) cat("ok.\n")

# extract empty space function F
  Fout <- attr(Jout, "F")
  fns[[1]] <- Fout[c("r","km","rs","raw","hazard","theo")]
  titles[[1]] <- "F function"
  if(verb) cat("F done.\n")

# extract Nearest neighbour distance distribution function G
  Gout <- attr(Jout, "G")
  fns[[2]] <- Gout[c("r","km","rs","hazard","theo")]
  titles[[2]] <- "G function"
  if(verb) cat("G done.\n")

# extract J function
  fns[[3]] <- Jout[c("r","km","rs","un","theo")]
  titles[[3]] <- "J function"
  if(verb) cat("J done.\n")

# compute second moment function K
  fns[[4]] <- Kest(pp,eps=NULL,breaks=brks)[c("r","border","theo")]
  titles[[4]] <- "K function"
  if(verb) cat("K done.\n")

# wrap into 'fasp' object
  
  deform <- list(cbind(km,theo)~r,cbind(km,theo)~r,
               cbind(km,theo)~r,cbind(border,theo)~r)
  witch <- matrix(1:4,2,2,byrow=T)

  if(is.null(dataname))
    dataname <- deparse(substitute(pp))
  title <- paste("Four summary functions for ",
              	dataname,".",sep="")

  rslt <- list(fns=fns,titles=titles,default.formula=deform,which=witch,
             dataname=dataname,title=title)
  class(rslt) <- "fasp"
  rslt
}
#
#      alltypes.R
#
#   $Revision: 1.4 $   $Date: 2001/11/25 00:32:04 $
#
#
alltypes <- function(pp, fun="K", dataname=NULL,verb=F) {
#
# Function 'alltypes' --- calculates a summary function for
# each type, or each pair of types, in a multitype point pattern
#
  verifyclass(pp,"ppp")

# validate 'fun'
  switch(fun, F={}, G={}, J={}, K={},
         stop("Unrecognized function name: ",fun,".\n"))
  wrong <- function(...) {stop("Internal error!")}

# list all possible types  
  if(!is.marked(pp)) {
    um <- 1
    nm <- 1
  } else {
    if(!is.factor(pp$marks))
      stop("the marks must be a factor")
    um <- levels(pp$marks)
    nm <- length(um)
  }
  
# get sensible 'r' values
  brks <- handle.r.b.args(r=NULL, breaks=NULL, pp$window, eps=NULL)
  r    <- brks$r

  # select appropriate statistics
  F1 <- switch(fun,F=Fest,G=Gest,J=Jest,K=Kest, wrong)
  F2 <- switch(fun,F={},G=Gcross,J=Jcross,K=Kcross, wrong)
  RF  <- switch(fun,F=c("r","km","rs","raw","hazard","theo"),
                   G=c("r","km","rs","hazard", "theo"),
                   J=c("r","km","rs","un","theo"),
                   K=c("r","border","theo"),
                   NA)
  deform <- switch(fun,F=cbind(km,theo)~r,
                      G=cbind(km,theo)~r,
                      J=cbind(km,theo)~r,
                      K=cbind(border,theo)~r,
                      NA)

# initialise 'fasp' object
  fns  <- list()
  
  if(fun=="F") {
    witch <- matrix(1:nm,ncol=1,nrow=nm)
    titles <- if(nm > 1) as.list(paste("mark =", um)) else list("")
  } else {
    witch <- matrix(1:(nm^2),ncol=nm,nrow=nm,byrow=T)
    titles <- if(nm > 1)
      as.list(paste("(", um[t(row(witch))], ",", um[t(col(witch))], ")", sep=""))
    else
      list("")
  }

  # compute function array
  k   <- 0

  for(i in 1:nrow(witch)) {
	for(j in 1:ncol(witch)) {
          if(verb) cat("i =",i,"j =",j,"\n")
          k <- k+1
          fns[[k]] <-
            if(nm == 1) # univariate pattern
              F1(pp,eps=NULL,breaks=brks)[RF]
            else if(fun=="F" | i==j) # F_i or G_ii, J_ii, K_ii
              F1(pp[pp$marks==um[i]], eps=NULL,breaks=brks)[RF]
            else 
              F2(pp,um[i],um[j], eps=NULL, breaks=brks)[RF]
        }
      }

  # wrap up into 'fasp' object
  if(is.null(dataname)) dataname <- deparse(substitute(pp))

  if(nm > 1)
	title <- paste("Array of ",fun," functions for ",
              	dataname,".",sep="")
  else
	title <- paste(fun," function for ",dataname,".",sep="")

  rslt <- list(fns=fns,titles=titles,default.formula=deform,which=witch,
             dataname=dataname,title=title)
  class(rslt) <- "fasp"
  rslt
}
#
#	breakpts.S
#
#	A simple class definition for the specification
#       of histogram breakpoints in the special form we need them.
#
#	even.breaks()
#
#	$Revision: 1.5 $	$Date: 2001/11/16 09:03:00 $
#
#
#       Other functions in this directory use the standard Splus function
#	hist() to compute histograms of distance values.
#       One argument of hist() is the vector 'breaks'
#	of breakpoints for the histogram cells. 
#
#       The breakpoints must
#            (a) span the range of the data
#            (b) be given in increasing order
#            (c) satisfy breaks[2] = 0,
#
#	The function make.even.breaks() will create suitable breakpoints.
#
#       Condition (c) means that the first histogram cell has
#       *right* endpoint equal to 0.
#
#       Since all our distance values are nonnegative, the effect of (c) is
#       that the first histogram cell counts the distance values which are
#       exactly equal to 0. Hence F(0), the probability P{X = 0},
#       is estimated without a discretisation bias.
#
#	We assume the histograms have followed the default counting rule
#	in hist(), which is such that the k-th entry of the histogram
#	counts the number of data values in 
#		I_k = ( breaks[k],breaks[k+1] ]	for k > 1
#		I_1 = [ breaks[1],breaks[2]   ]
#
#	The implementations of estimators of c.d.f's in this directory
#       produce vectors of length = length(breaks)-1
#       with value[k] = estimate of F(breaks[k+1]),
#       i.e. value[k] is an estimate of the c.d.f. at the RIGHT endpoint
#       of the kth histogram cell.
#
#       An object of class 'breakpts' contains:
#
#              $val     the actual breakpoints
#              $max     the maximum value (= last breakpoint)
#              $ncells  total number of histogram cells
#              $r       right endpoints, r = val[-1]
#              $even    logical = TRUE if cells known to be evenly spaced
#              $npos    number of histogram cells on the positive halfline
#                        = length(val) - 2,
#                       or NULL if cells not evenly spaced
#              $step    histogram cell width
#                       or NULL if cells not evenly spaced
#       
# --------------------------------------------------------------------
breakpts <- function(val, maxi, even=F, npos=NULL, step=NULL) {
  out <- list(val=val, max=maxi, ncells=length(val)-1, r = val[-1],
              even=even, npos=npos, step=step)
  class(out) <- "breakpts"
  out
}

"make.even.breaks" <- 
function(bmax, npos, bstep) {
        if(missing(bstep) && missing(npos))
          stop("Must specify either \'bstep\' or \'npos\'")
        if(!missing(npos)) {
          bstep <- bmax/npos
          val <- seq(0, bmax, length=npos+1)
          val <- c(-bstep,val)
        } else {
          val <- seq(0, bmax, by=bstep)
          val <- c(-bstep,val)
          npos <- length(val) - 2
        }
        breakpts(val, bmax, T, npos, bstep)
}

"as.breakpts" <- function(...) {

  XL <- list(...)

  if(length(XL) == 1) {
    # single argument
    X <- XL[[1]]

    if(!is.null(class(X)) && class(X) == "breakpts")
    # X already in correct form
      return(X)
  
    if(is.vector(X) && length(X) > 2) {
    # it's a vector
      if(X[2] != 0)
        stop("breakpoints do not satisfy breaks[2] = 0")
      steps <- diff(X)
      if(all(steps == mean(steps)))
        # equally spaced
        return(breakpts(X, max(X), T, length(X)-2, steps[1]))
      else
        # unknown spacing
        return(breakpts(X, max(X), F))
    }
  } else {

    # There are multiple arguments.
  
    # exactly two arguments - interpret as even.breaks()
    if(length(XL) == 2)
      return(make.even.breaks(XL[[1]], XL[[2]]))

    # two arguments 'max' and 'npos'
  
    if(!is.null(XL$max) && !is.null(XL$npos))
      return(make.even.breaks(XL$max, XL$npos))

    # otherwise
    stop("Don't know how to convert these data to breakpoints")
  }
  # never reached
}


check.hist.lengths <- function(hist, breaks) {
  verifyclass(breaks, "breakpts")
  nh <- length(hist)
  nb <- breaks$ncells
  if(nh != nb)
    stop(paste("Length of histogram =", nh,
               "not equal to number of histogram cells =", nb))
}

breakpts.from.r <- function(r) {
        if(r[1] != 0)
          stop("First r value must be 0")
        dr <- r[2] - r[1]
        b <- c(-dr, r)
        return(as.breakpts(b))
}

handle.r.b.args <- function(r=NULL, breaks=NULL, window, eps=NULL) {

        if(!is.null(r) && !is.null(breaks))
          stop("Do not specify both \'r\' and \'breaks\'")
  
        if(!is.null(breaks)) {
          breaks <- as.breakpts(breaks)
        } else if(!is.null(r)) {
          breaks <- breakpts.from.r(r)
	} else {
          # both 'r' and 'breaks' are missing
          if(is.null(eps)) {
            if(!is.null(window$xstep))
              eps <- window$xstep
            else 
              eps <- diff(window$xrange)/100
          }
          # warning(paste("step size for argument \'r\' defaults to", eps/4))
          breaks <- make.even.breaks( diameter(window), bstep=eps/4)
        }

        return(breaks)
}
#
#
#	classes.S
#
#	$Revision: 1.3 $	$Date: 2000/07/11 10:53:19 $
#
#	Generic utilities for classes
#
#	These are supposed to work under BOTH Splus3.4 and Splus5.1
#
#
#--------------------------------------------------------------------------

verifyclass <- function(X, C, N=deparse(substitute(X)), fatal=T) {
  if(is.null(class(X)) || class(X) != C) {
    if(fatal) {
        gripe <- paste("argument \'", N,
                       "\' is not of class \'", C, "\'", sep="")
	stop(gripe)
    } else 
	return(F)
  }
  return(T)
}

#--------------------------------------------------------------------------

checkfields <- function(X, L) {
	  # X is a list, L is a vector of strings
	  # Checks for presence of field named L[i] for all i
	return(all(!is.na(match(L,names(X)))))
}

getfields <- function(X, L, fatal=T) {
	  # X is a list, L is a vector of strings
	  # Extracts all fields with names L[i] from list X
	  # Checks for presence of all desired fields
	  # Returns the sublist of X with fields named L[i]
	absent <- is.na(match(L, names(X)))
	if(any(absent)) {
		gripe <- paste("Needed the following components:",
				paste(L, collapse=", "),
				"\nThese ones were missing: ",
				paste(L[absent], collapse=", "))
		if(fatal)
			stop(gripe)
		else 
			warning(gripe)
	} 
	return(X[L[!absent]])
}



#
#       conspire.S
#
#  $Revision: 1.1 $    $Date: 2001/11/05 05:50:40 $
#
#

conspire <- function(indata, formula, subset=NULL, lty=NULL, col=NULL) {

  lhs <- formula[[2]]
  rhs <- formula[[3]]

  # evaluate expression a in data frame b
  evaluate <- function(a,b) {
    if(exists("is.R") && is.R())
      eval(a, envir=b)
    else
      eval(a, local=b)
  }
  
  lhsdata <- evaluate(lhs, indata)
  rhsdata <- evaluate(rhs, indata)
  
  if(is.vector(lhsdata))
    lhsdata <- matrix(lhsdata, ncol=1)

  if(!is.vector(rhsdata))
    stop("rhs of formula seems not to be a vector")

  if(!is.null(subset)) {
    keep <- if(is.character(subset))
		evaluate(parse(text=subset), indata)
            else
                evaluate(subset, indata)
    lhsdata <- lhsdata[keep, , drop=F]
    rhsdata <- rhsdata[keep]
  }

  ylim <- range(lhsdata,na.rm=T)
  xlim <- range(rhsdata,na.rm=T)

  # work out how to label the plot
  yname <- paste(lhs)
  if(length(yname) > 1 && yname[[1]] == "cbind")
    ylab <- paste(yname[-1], collapse=" , ")
  else
    ylab <- as.character(formula)[2]
  
  xlab <- as.character(formula)[3]

  plot(xlim, ylim, type="n", xlab=xlab, ylab=ylab)

  nplots <- ncol(lhsdata)

  if(is.null(lty))
    lty <- 1:nplots
  if(is.null(col))
    col <- rep(1,nplots)
  
  for(i in 1:nplots)
    lines(rhsdata, lhsdata[,i], lty=lty[i], col=col[i])
}
#
#	distbdry.S		Distance to boundary
#
#	$Revision: 4.2 $	$Date: 2001/08/07 09:22:52 $
#
# -------- functions ----------------------------------------
#
#	bdist.points()
#                       compute vector of distances 
#			from each point of point pattern
#                       to boundary of window
#
#       bdist.pixels()
#                       compute matrix of distances from each pixel
#                       to boundary of window
#
#       erode.mask()    erode the window mask by a distance r
#                       [yields a new window]
#
#       erode.owin()    erode the window by a distance r
#                       [yields a new window]
#
#
# 
"bdist.points"<-
function(X)
{
	verifyclass(X, "ppp") 
	x <- X$x
	y <- X$y
	window <- X$window
        switch(window$type,
               rectangle = {
		xmin <- min(window$xrange)
		xmax <- max(window$xrange)
		ymin <- min(window$yrange)
		ymax <- max(window$yrange)
		result <- pmin(x - xmin, xmax - x, y - ymin, ymax - y)
               },
               polygonal = {
                 xy <- cbind(x,y)
                 result <- rep(Inf, X$n)
                 bdry <- window$bdry
                 for(i in seq(bdry)) {
                   polly <- bdry[[i]]
                   nsegs <- length(polly$x)
                   for(j in 1:nsegs) {
                     j1 <- if(j < nsegs) j + 1 else 1
                     seg <- c(polly$x[j],  polly$y[j],
                              polly$x[j1], polly$y[j1])
                     result <- pmin(result, distppl(xy, seg))
                   }
                 }
               },
               mask = {
                 b <- bdist.pixels(window, coords=F)
                 loc <- nearest.raster.point(x,y,window)
                 result <- numeric(X$n)
                 for(i in 1:(X$n)) 
			result[i] <- b[loc$row[i],loc$col[i]]
               },
               stop("Unrecognised window type", window$type)
               )
        return(result)
}

"bdist.pixels" <- function(w, ..., coords=T) {
	verifyclass(w, "owin")

        masque <- as.mask(w, ...)
        
        switch(w$type,
               mask = {
                 neg <- complement.owin(masque)
                 m <- exactPdt(neg)
                 b <- pmin(m$d,m$b)
               },
               rectangle = {
                 x <- raster.x(masque)
                 y <- raster.y(masque)
                 xmin <- w$xrange[1]
                 xmax <- w$xrange[2]
                 ymin <- w$yrange[1]
                 ymax <- w$yrange[2]
                 b <- pmin(x - xmin, xmax - x, y - ymin, ymax - y)
               },
               polygonal = {
                 # set up pixel raster
                 x <- as.vector(raster.x(masque))
                 y <- as.vector(raster.y(masque))
                 b <- rep(0, length(x))
                 # test each pixel in/out, analytically
                 inside <- inside.owin(x, y, w)
                 # compute distances for these pixels
                 xy <- cbind(x[inside], y[inside])
                 dxy <- rep(Inf, sum(inside))
                 bdry <- w$bdry
                 for(i in seq(bdry)) {
                   polly <- bdry[[i]]
                   nsegs <- length(polly$x)
                   for(j in 1:nsegs) {
                     j1 <- if(j < nsegs) j + 1 else 1
                     seg <- c(polly$x[j],  polly$y[j],
                              polly$x[j1], polly$y[j1])
                     dxy <- pmin(dxy, distppl(xy, seg))
                   }
                 }
                 b[inside] <- dxy
               },
               stop("unrecognised window type", w$type)
               )

        # reshape it
        b <- matrix(b, nrow=masque$dim[1], ncol=masque$dim[2])
        
        if(coords)
          # return in a format which can be plotted by image(), persp() etc
          return(list(x=masque$xcol, y=masque$yrow, z=t(b)))
        else
          # return matrix (for internal use by package)
          return(b)
} 

erode.mask <- function(w, r) {
  # erode a binary image mask without changing any other entries
  
	verifyclass(w, "owin")
        if(w$type != "mask")
          stop("window w is not of type \'mask\'")
        
	bb <- bdist.pixels(w, coords=F)

        if(r > max(bb))
          warning("eroded mask is empty")

        w$m <- (bb >= r)
        return(w)
}

        
erode.owin <- function(w, r, shrink.frame=T, ...) {
	verifyclass(w, "owin")

        if(2 * r >= max(diff(w$xrange), diff(w$yrange)))
          stop("erosion distance r too large for frame of window")

        # compute the dimensions of the eroded frame
        exr <- w$xrange + c(r, -r)
        eyr <- w$yrange + c(r, -r)
        ebox <- list(x=exr[c(1,2,2,1)], y=eyr[c(1,1,2,2)])

        if(w$type == "rectangle") {
          # result is a smaller rectangle
          if(shrink.frame)
            return(owin(exr, eyr))  # type 'rectangle' 
          else
            return(owin(w$xrange, w$yrange, poly=ebox)) # type 'polygonal'
        }

        # otherwise erode the window in pixel image form
        if(w$type != "mask") 
          w <- as.mask(w, ...)
        
        wnew <- erode.mask(w, r)

        if(shrink.frame) {
          # trim off some rows & columns of pixel raster
          keepcol <- (wnew$xcol >= exr[1] & wnew$xcol <= exr[2])
          keeprow <- (wnew$yrow >= eyr[1] & wnew$yrow <= eyr[2])
          wnew$xcol <- w$xcol[keepcol]
          wnew$yrow <- w$yrow[keepcol]
          wnew$dim <- c(sum(keeprow), sum(keepcol))
          wnew$m <- wnew$m[keeprow, keepcol]
        }

        return(wnew)
}	
#
#	dummy.S
#
#	Utilities for generating patterns of dummy points
#
#       $Revision: 4.3 $     $Date: 2001/12/11 15:29:17 $
#
#	corners()	corners of window
#	gridcenters()	points of a rectangular grid
#	stratrand()	random points in each tile of a rectangular grid
#	spokes()	Rolf's 'spokes' arrangement
#	
#	concatxy()	concatenate any lists of x, y coordinates
#
#	default.dummy()	Default action to create a dummy pattern
#		
	
corners <- function(window) {
	window <- as.owin(window)
	x <- window$xrange[c(1,2,1,2)]
	y <- window$yrange[c(1,1,2,2)]
	return(list(x=x, y=y))
}

gridcenters <-	
gridcentres <- function(window, nx, ny) {
	window <- as.owin(window)
	xr <- window$xrange
	yr <- window$yrange
	x <- seq(xr[1], xr[2], length = 2 * nx + 1)[2 * (1:nx)]
	y <- seq(yr[1], yr[2], length = 2 * ny + 1)[2 * (1:ny)]
	x <- rep(x, ny)
	y <- rep(y, rep(nx, ny))
	return(list(x=x, y=y))
}

stratrand <- function(window,nx,ny, k=1) {
	
	# divide window into an nx * ny grid of tiles
	# and place k points at random in each tile
	
	window <- as.owin(window)

	wide  <- diff(window$xrange)/(nx - 1)
	high  <- diff(window$yrange)/(ny - 1)
        cent <- gridcentres(window, nx, ny)
	cx <- rep(cent$x, k)
	cy <- rep(cent$y, k)
	n <- nx * ny * k
	x <- cx + runif(n, min = -wide/2, max = wide/2)
	y <- cy + runif(n, min = -high/2, max = high/2)
	return(list(x=x,y=y))
}

tilecentroids <- function (W, nx, ny)
{
  W <- as.owin(W)
  if(W$type == "rectangle")
    return(gridcentres(W, nx, ny))
  else {
    # approximate
    W   <- as.mask(W)
    xx  <- as.vector(raster.x(W)[W$m])
    yy  <- as.vector(raster.y(W)[W$m])
    pid <- gridindex(xx,yy,W$xrange,W$yrange,nx,nx)$index
    x   <- tapply(xx,pid,mean)
    y   <- tapply(yy,pid,mean)
    return(list(x=x,y=y))
  }
}

tilemiddles <- function (W, nx, ny)
{
    if(W$type == "rectangle")
      return(gridcentres(W, nx, ny))
    
    # pixel approximation to window
    # This matches the pixel approximation used to compute tile areas
    # and ensures that dummy points are generated only inside those tiles
    # that have nonzero digital area
    M   <- as.mask(W)
    xx  <- as.vector(raster.x(M))[M$m]
    yy  <- as.vector(raster.y(M))[M$m]
    # label the pixels by tile index
    pid <- gridindex(xx,yy,W$xrange,W$yrange,nx,ny)$index
    pidvalues <- levels(as.factor(pid))
    # Latter is for comparison with the output of tapply( ..., pid, ...)
    
    ######## 1st try: centroids of tiles
    # (usually inside tile, e.g. if it's convex, but not always inside)
    cx <- tapply(xx, pid, mean)
    cy <- tapply(yy, pid, mean)
    # centroids which are inside window
    cok <- inside.owin(cx, cy, W)
    x <- cx[cok]
    y <- cy[cok]
    # tiles for which this strategy did not work
    nbg <- pidvalues[!cok]
    
    ######### 2nd try: middle point in list of pixels in each tile
    # (always inside tile, by construction)
    todo <- pid %in% nbg
    xx <- xx[todo]
    yy <- yy[todo]
    pid <- pid[todo]
    middle <- function(v) { n <- length(v); mid <- ceiling(n/2); v[mid]}
    mx   <- tapply(xx,pid,middle)
    my   <- tapply(yy,pid,middle)
    
    x <- c( x, mx)
    y <- c( y, my)
    return(list(x=x,y=y))
}

spokes <- function(x, y, nrad = 3, nper = 3, fctr = 1.5) {
	#
	# Rolf Turner's "spokes" arrangement
	#
	# Places dummy points on radii of circles 
	# emanating from each data point x[i], y[i]
	#
	#       nrad:    number of radii from each data point
	#       nper:	 number of dummy points per radius
	#       fctr:	 length of largest radius = fctr * M
	#                where M is mean nearest neighbour distance in data
	#
	M <- mean(nndist(x,y))
	lrad  <- fctr * M / nper
	theta <- 2 * pi * (1:nrad)/nrad
	cs    <- cos(theta)
	sn    <- sin(theta)
	xt    <- lrad * as.vector((1:nper) %o% cs)
	yt    <- lrad * as.vector((1:nper) %o% sn)
	xd    <- as.vector(outer(x, xt, "+"))
	yd    <- as.vector(outer(y, yt, "+"))
	
	return(list(x=xd,y=yd))
}
	
# concatenate any number of list(x,y) into a list(x,y)
		
concatxy <- function(...) {
	x <- unlist(lapply(list(...), function(w) {w$x}))
	y <- unlist(lapply(list(...), function(w) {w$y}))
	if(length(x) != length(y))
		stop("Internal error: lengths of x and y unequal")
	return(list(x=x,y=y))
}

#------------------------------------------------------------

default.dummy <- function(X, nx=NULL, ny=NULL, random=F) {
	# default action to create dummy points.
	# regular grid of (order) nx * ny tiles
	# plus corner points of window frame.
	# 
	X <- as.ppp(X)
	win <- X$window
        # default dimensions
        if(is.null(nx)) nx <- default.ngrid(X)
        if(is.null(ny)) ny <- nx
        # make dummy points
        dummy <- if(random) 
                  stratrand(win, nx, ny, 1)
                else
                  tilemiddles(win, nx, ny)
        # add corner points
	dummy <- concatxy(dummy, corners(win))
	# restrict to window (if it's a mask)
	ok <- inside.owin(dummy$x, dummy$y, win)
	x <- dummy$x[ok]
	y <- dummy$y[ok]
	ppp(x, y, window=win)
}


default.ngrid <- function(X) {
	# default dimensions of rectangular grid
        # (for dummy points, tile weights etc)
  X <- as.ppp(X)
  min(30, 10 * ceiling(X$n/10))
}
#
#	exactPdt.S
#	S function exactPdt() for exact distance transform of pixel image
#
#	$Revision: 4.1 $	$Date: 2001/08/07 09:48:04 $
#

"exactPdt"<-
function(im)
{
#
        verifyclass(im, "owin")
        if(im$type != "mask")
          stop("Input must be a window of type \'mask\'")
#	
	nr <- im$dim[1]
	nc <- im$dim[2]
# pad out the input image with a margin of width 1 on all sides
	x <- im$m
	x <- cbind(F, x, F)
	x <- rbind(F, x, F)
#	
	res <- .C("ps_exact_dt_S",
		as.double(im$xrange[1]),
		as.double(im$yrange[1]),
		as.double(im$xrange[2]),
		as.double(im$yrange[2]),
		nr = as.integer(nr),
		nc = as.integer(nc),
		as.logical(t(x)),
		dd = as.double (matrix(0, ncol = nc + 2, nrow = nr + 2)),
		rr = as.integer(matrix(0, ncol = nc + 2, nrow = nr + 2)),
		cc = as.integer(matrix(0, ncol = nc + 2, nrow = nr + 2)),
		bb = as.double (matrix(0, ncol = nc + 2, nrow = nr + 2))
		)
	dist <- matrix(res$dd, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]
	rows <- matrix(res$rr, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]
	cols <- matrix(res$cc, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]
	bdist<- matrix(res$bb, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]

        # convert from C to S
        rows <- rows + 1
        cols <- cols + 1
        
	out <- im
	out$m <- NULL
	out <- append(out, list(d=dist,row=rows,col=cols,b=bdist))
	invisible(out)
}
#
#	exactdt.S
#	S function exactdt() for exact distance transform
#
#	$Revision: 4.1 $	$Date: 2001/08/07 09:22:52 $
#

"exactdt"<-
function(X, ...)
{
	verifyclass(X, "ppp")

        w <- as.mask(X$window, ...)

#
	nr <- w$dim[1]
	nc <- w$dim[2]
#	
	res <- .C("exact_dt_S",
		as.double(X$x),
		as.double(X$y),
		as.integer(X$n),
		as.double(w$xrange[1]),
		as.double(w$yrange[1]),
		as.double(w$xrange[2]),
		as.double(w$yrange[2]),
		nr = as.integer(nr),
		nc = as.integer(nc),
		d = as.double(matrix(0, ncol = nc + 2, nrow = nr + 2)),
		i = as.integer(matrix(0, ncol = nc + 2, nrow = nr + 2)),
		b = as.double(matrix(0, ncol = nc + 2, nrow = nr + 2)))
	dist <- matrix(res$d, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]
	inde <- matrix(res$i, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]
	bdry <- matrix(res$b, ncol = nc + 2, byrow = T)[2:(nr + 1), 2:(nc +1)]

        inde <- inde + 1    # convert from C to S indexing
        
	result <- list(d = dist, i = inde, b = bdry)
	append(result, w)
}
#
#	fasp.R
#
#	$Revision: 1.3 $	$Date: 2002/01/18 06:52:33 $
#
#
#-----------------------------------------------------------------------------
#

"[.fasp" <-
"subset.fasp" <-
  function(x, I, J, drop, ...) {

        verifyclass(x, "fasp")
        
        m <- nrow(x$which)
        n <- ncol(x$which)
        
        if(missing(I)) I <- 1:m
        if(missing(J)) J <- 1:n

        # determine index subset for lists 'fns', 'titles' etc
        included <- rep(F, length(x$fns))
        w <- as.vector(x$which[I,J])
        included[w] <- T

        # determine positions in shortened lists
        newk <- cumsum(included)

        # assign result
        Y <- list()
        Y$fns <- x$fns[included]
        Y$titles <- x$titles[included]
        xdf <- x$default.formula
        Y$default.formula <- if(any(class(xdf)=="formula")) xdf else xdf[included]
        oldwhich <- x$which[I,J,drop=F]
        newwhich <- newk[oldwhich]
        Y$which <- matrix(newwhich, ncol=ncol(oldwhich), nrow=nrow(oldwhich))
        Y$dataname <- x$dataname
        Y$title <- x$title
        class(Y) <- "fasp"
        
        return(Y)
}
#
#
#   formulae.S
#
#   Functions for manipulating model formulae
#
#	$Revision: 1.6 $	$Date: 2002/01/18 06:57:18 $
#
#   identical.formulae()
#          Check whether two formulae are identical
#
#   termsinformula()
#          Extract the terms from a formula
#
#   sympoly()
#          Create a symbolic polynomial formula
#
#   polynom()
#          Analogue of poly() but without dynamic orthonormalisation
#
# -------------------------------------------------------------------
#	

identical.formulae <- function(a, b) {
   if(is.null(a))
      return(is.null(b))
   if(is.null(b))
      return(is.null(a))
   if(class(a) != "formula" || class(b) != "formula")
      stop("one of the arguments is not a formula")
   ac <- as.character(a)
   bc <- as.character(b)
   return( (length(ac) == length(bc)) && all(ac == bc))
}

termsinformula <- function(x) {
  if(is.null(x)) return(character(0))
  if(class(x) != "formula")
    stop("argument is not a formula")
  attr(terms(x), "term.labels")
}

sympoly <- function(x,y,n) {

   if(nargs()<2) stop("Degree must be supplied.")
   if(nargs()==2) n <- y
   eps <- abs(n%%1)
   if(eps > 0.000001 | n <= 0) stop("Degree must be a positive integer")
   
   x <- deparse(substitute(x))
   temp <- NULL
   left <- "I("
   rght <- ")"
   if(nargs()==2) {
	for(i in 1:n) {
		xhat <- if(i==1) "" else paste("^",i,sep="")
		temp <- c(temp,paste(left,x,xhat,rght,sep=""))
	}
   }
   else {
	y <- deparse(substitute(y))
	for(i in 1:n) {
		for(j in 0:i) {
			k <- i-j
			xhat <- if(k<=1) "" else paste("^",k,sep="")
			yhat <- if(j<=1) "" else paste("^",j,sep="")
			xbit <- if(k>0) x else ""
			ybit <- if(j>0) y else ""
			star <- if(j*k>0) "*" else ""
			term <- paste(left,xbit,xhat,star,ybit,yhat,rght,sep="")
			temp <- c(temp,term)
		}
	}
      }
   as.formula(paste("~",paste(temp,collapse="+")))
 }


polynom <- function(x, ...) {
  rest <- list(...)
  # degree not given
  if(length(rest) == 0)
    stop("degree of polynomial must be given")
  #call with single variable and degree
  if(length(rest) == 1) {
    degree <- ..1
    if((degree %% 1) != 0 || length(degree) != 1 || degree < 1)
      stop("degree of polynomial should be a positive integer")

    # compute values
    result <- outer(x, 1:degree, "^")

    # compute column names - the hard part !
    namex <- deparse(substitute(x))
    # check whether it needs to be parenthesised
    if(!is.name(substitute(x))) 
      namex <- paste("(", namex, ")", sep="")
    # column names
    namepowers <- if(degree == 1) namex else 
                       c(namex, paste(namex, "^", 2:degree, sep=""))
    namepowers <- paste("[", namepowers, "]", sep="")
    # stick them on
    dimnames(result) <- list(NULL, namepowers)
    return(result)
  }
  # call with two variables and degree
  if(length(rest) == 2) {

    y <- ..1
    degree <- ..2

    # list of exponents of x and y, in nice order
    xexp <- yexp <- numeric()
    for(i in 1:degree) {
      xexp <- c(xexp, i:0)
      yexp <- c(yexp, 0:i)
    }
    nterms <- length(xexp)
    
    # compute 

    result <- matrix(, nrow=length(x), ncol=nterms)
    for(i in 1:nterms) 
      result[, i] <- x^xexp[i] * y^yexp[i]

    #  names of these terms
    
    namex <- deparse(substitute(x))
    # namey <- deparse(substitute(..1)) ### seems not to work in R
    zzz <- as.list(match.call())
    namey <- deparse(zzz[[3]])

    # check whether they need to be parenthesised
    # if so, add parentheses
    if(!is.name(substitute(x))) 
      namex <- paste("(", namex, ")", sep="")
    if(!is.name(zzz[[3]])) 
      namey <- paste("(", namey, ")", sep="")

    nameXexp <- c("", namex, paste(namex, "^", 2:degree, sep=""))
    nameYexp <- c("", namey, paste(namey, "^", 2:degree, sep=""))

    # make the term names
       
    termnames <- paste(nameXexp[xexp + 1],
                       ifelse(xexp > 0 & yexp > 0, ".", ""),
                       nameYexp[yexp + 1],
                       sep="")
    termnames <- paste("[", termnames, "]", sep="")

    dimnames(result) <- list(NULL, termnames)
    # 
    return(result)
  }
  stop("Can't deal with more than 2 variables yet")
}
#
#
#    geyer.S
#
#    $Revision: 1.2 $	$Date: 2001/07/27 03:56:05 $
#
#    Geyer's saturation process
#
#    Geyer()    create an instance of Geyer's saturation process
#                 [an object of class 'interact']
#
# ------------------------------------------------------------------
#    Note: if you want to imitate this, remember that 'pairsat.family'
#    expects the saturation parameter 'sat' to be called $par$saturate
#    in this 'interact' object.
# -------------------------------------------------------------------
#	

Geyer <- function(r, sat) {
  out <- 
  list(
         name     = "Geyer saturation process",
         family    = pairsat.family,
         pot      = function(d, par) {
                         ifelse(d <= par$r, 1, 0)  # same as for Strauss
                    },
         par      = list(r = r, saturate=sat),
         parnames = c("interaction distance","saturation parameter"),
         init     = function(self) {
                      r <- self$par$r
                      sat <- self$par$saturate
                      if(!is.numeric(r) || length(r) != 1 || r <= 0)
                       stop("interaction distance r must be a positive number")
                      if(!is.numeric(sat) || length(sat) != 1 || sat < 1)
                       stop("saturation parameter sat must be a number >= 1")
                    },
         update = NULL,  # default OK
         print = NULL    # default OK
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#	interact.S
#
#
#	$Revision: 1.6 $	$Date: 2002/01/18 06:46:36 $
#
#	Class 'interact' representing the interpoint interaction
#               of a point process model
#              (e.g. Strauss process with a given threshold r)
#
#       Class 'isf' representing a generic interaction structure
#              (e.g. pairwise interactions)
#
#	These do NOT specify the "trend" part of the model,
#	only the "interaction" component.
#
#               The analogy is:
#
#                       glm()             mpl()
#
#                       model formula     trend formula
#
#                       family            interaction
#
#               That is, the 'systematic' trend part of a point process
#               model is specified by a 'trend' formula argument to mpl(),
#               and the interpoint interaction is specified as an 'interact'
#               object.
#
#       You only need to know about these classes if you want to
#       implement a new point process model.
#
#       THE DISTINCTION:
#       An object of class 'isf' describes an interaction structure
#       e.g. pairwise interaction, triple interaction,
#       pairwise-with-saturation, Dirichlet interaction.
#       Think of it as determining the "order" of interaction
#       but not the specific interaction potential function.
#
#       An object of class 'interact' completely defines the interpoint
#       interactions in a specific point process model, except for the
#       regular parameters of the interaction, which are to be estimated
#       by mpl() or otherwise. An 'interact' object specifies the values
#       of all the 'nuisance' or 'irregular' parameters. An example
#       is the Strauss process with a given, fixed threshold r
#       but with the parameters beta and gamma undetermined.
#
#       DETAILS:
#
#       An object of class 'isf' contains the following:
#
#	     $name               Name of the interaction structure         
#                                        e.g. "pairwise"
#
#	     $print		 How to 'print()' this object
#				 [A function; invoked by the 'print' method
#                                 'print.isf()']
#
#            $eval               A function which evaluates the canonical
#                                sufficient statistic for an interaction
#                                of this general class (e.g. any pairwise
#                                interaction.)
#
#       If lambda(u,X) denotes the conditional intensity at a point u
#       for the point pattern X, then we assume
#                  log lambda(u, X) = theta . S(u,X)
#       where theta is the vector of regular parameters,
#       and we call S(u,X) the sufficient statistic.
#
#       A typical calling sequence for the $eval function is
#
#            (f$eval)(X, U, E, potentials, potargs, correction)
#
#       where X is the data point pattern, U is the list of points u
#       at which the sufficient statistic S(u,X) is to be evaluated,
#       E is a logical matrix equivalent to (X[i] == U[j]),
#       $potentials defines the specific potential function(s) and
#       $potargs contains any nuisance/irregular parameters of these
#       potentials [the $potargs are passed to the $potentials without
#       needing to be understood by $eval.]
#       $correction is the name of the edge correction method.
#
#
#       An object of class 'interact' contains the following:
#
#
#            $name               Name of the specific potential
#                                        e.g. "Strauss"
#
#            $family              Object of class "isf" describing
#                                the interaction structure
#
#            $pot	         The interaction potential function(s)
#                                -- usually a function or list of functions.
#                                (passed as an argument to $family$eval)
#
#            $par                list of any nuisance/irregular parameters
#                                (passed as an argument to $family$eval)
#
#            $parnames           vector of long names/descriptions
#                                of the parameters in 'par'
#
#            $init()             initialisation action
#                                or NULL indicating none required
#
#            $update()           A function to modify $par
#                                [Invoked by 'update.interact()']
#                                or NULL indicating a default action
#
#	     $print		 How to 'print()' this object
#				 [Invoked by 'print' method 'print.interact()']
#                                or NULL indicating a default action
#
# --------------------------------------------------------------------------

print.isf <- function(x, ...) {
  verifyclass(x, "isf")
  if(!is.null(x$print))
    (x$print)(x)
  invisible(NULL)
}

print.interact <- function(x, ...) {
  verifyclass(x, "interact")
  if(!is.null(x$print))
     (x$print)(x)
  else {
    # default
    print.isf(x$family)
    cat(paste("Interaction:", x$name, "\n"))
    # just print the parameter names and their values
    cat(paste(x$parnames, ":\t", x$par, "\n", sep=""))
  }
  invisible(NULL)
}

update.interact <- function(object, ...) {
  verifyclass(object, "interact")
  if(!is.null(object$update))
    (object$update)(object, ...)
  else {
    # Default
    # just match the arguments in "..."
    # with those in object$par and update them
    want <- list(...)
    m <- match(names(want),names(object$par))
    nbg <- is.na(m)
    if(any(nbg)) {
        which <- paste((names(want))[nbg])
        warning(paste("Arguments not matched: ", which))
    }
    m <- m[!nbg]
    object$par[m] <- u
    # call object's own initialisation routine
    if(!is.null(object$init))
      (object$init)(object)
    object
  }    
}

#
#	kmrs.S
#
#	S code for Kaplan-Meier and reduced sample
#	estimates of a distribution function
#	from _histograms_ of censored data.
#
#	kaplan.meier()
#	reduced.sample()
#       km.rs()
#
#	$Revision: 3.3 $	$Date: 2000/06/22 13:19:43 $
#
#	The functions in this file produce vectors `km' and `rs'
#	where km[k] and rs[k] are estimates of F(breaks[k+1]),
#	i.e. an estimate of the c.d.f. at the RIGHT endpoint of the interval.
#

"kaplan.meier" <-
function(obs, nco, breaks) {
#	obs: histogram of all observations : min(T_i,C_i)
#	nco: histogram of noncensored observations : T_i such that T_i <= C_i
# 	breaks: breakpoints (vector or 'breakpts' object, see breaks.S)
#
        breaks <- as.breakpts(breaks)

	n <- length(obs)
	if(n != length(nco)) 
		stop("lengths of histograms do not match")
	check.hist.lengths(nco, breaks)
#
#	
#   reverse cumulative histogram of observations
	d <- cumsum(obs[n:1])[n:1]
#
#  product integrand
	s <- ifelse(d > 0, 1 - nco/d, 1)
#
	km <- 1 - cumprod(s)
#  km has length n;  km[i] is an estimate of F(r) for r=breaks[i+1]
#	
	widths <- diff(breaks$val)
	lambda <-  - log(ifelse(s > 0, s, 1))/widths 
#  lambda has length n; lambda[i] is an estimate of
#  the average of \lambda(r) over the interval (breaks[i],breaks[i+1]).
#	
	return(list(km=km, lambda=lambda))
}

"reduced.sample" <-
function(nco, cen, ncc, show=F)
#	nco: histogram of noncensored observations: T_i such that T_i <= C_i
#	cen: histogram of all censoring times: C_i
#	ncc: histogram of censoring times for noncensored obs:
#		C_i such that T_i <= C_i
#
#	Then nco[k] = #{i: T_i <= C_i, T_i \in I_k}
#	     cen[k] = #{i: C_i \in I_k}
#	     ncc[k] = #{i: T_i <= C_i, C_i \in I_k}.
#
{
	n <- length(nco)
	if(n != length(cen) || n != length(ncc))
		stop("histogram lengths do not match")
#
#	denominator: reverse cumulative histogram of censoring times
#		denom(r) = #{i : C_i >= r}
#	We compute 
#		cc[k] = #{i: C_i > breaks[k]}	
#	except that > becomes >= for k=0.
#
	cc <- cumsum(cen[n:1])[n:1]
#
#
#	numerator
#	#{i: T_i <= r <= C_i }
#	= #{i: T_i <= r, T_i <= C_i} - #{i: C_i < r, T_i <= C_i}
#	We compute
#		u[k] = #{i: T_i <= C_i, T_i <= breaks[k+1]}
#			- #{i: T_i <= C_i, C_i <= breaks[k]}
#		     = #{i: T_i <= C_i, C_i > breaks[k], T_i <= breaks[k+1]}
#	this ensures that numerator and denominator are 
#	comparable, u[k] <= cc[k] always.
#
	u <- cumsum(nco) - c(0,cumsum(ncc)[1:(n-1)])
	rs <- u/cc
#
#	Hence rs[k] = u[k]/cc[k] is an estimator of F(r) 
#	for r = breaks[k+1], i.e. for the right hand end of the interval.
#
        if(!show)
          return(rs)
        else
          return(list(rs=rs, numerator=u, denominator=cc))
}

"km.rs" <-
function(o, cc, d, breaks) {
#	o: censored lifetimes min(T_i,C_i)
#	cc: censoring times C_i
#	d: censoring indicators 1(T_i <= C_i)
#	breaks: histogram breakpoints (vector or 'breakpts' object)
#
        breaks <- as.breakpts(breaks)
# compile histograms
	obs <- hist( o,		breaks=breaks$val,plot=F,probability=F)$counts
	nco <- hist( o[d], 	breaks=breaks$val,plot=F,probability=F)$counts
	cen <- hist( cc,	breaks=breaks$val,plot=F,probability=F)$counts
	ncc <- hist( cc[d],	breaks=breaks$val,plot=F,probability=F)$counts
# go
	km <- kaplan.meier(obs, nco, breaks)
	rs <- reduced.sample(nco, cen, ncc)
#
	return(list(rs=rs, km=km$km, hazard=km$lambda,
                    r=breaks$val[-1], breaks=breaks$val))
}
#
#	$Revision: 4.6 $	$Date: 2001/12/11 15:29:17 $
#
#    mpl()
#          Fit a point process model to a two-dimensional point pattern
#          The model may include
#                  - trend (arbitrary);
#                  - dependence on covariates;  
#                   - arbitrary interaction structure
#			(with p-dimensional regular parameter)
#	   The window of observation may have arbitrary shape.
#
#
# -------------------------------------------------------------------
#	
"mpl" <- 
function(Q,
         trend = ~1,
	 interaction = NULL,
         data,
	 correction="border",
	 rbord = 0,
         use.gam=F
) {
#
# Extract quadrature scheme 
#
	if(verifyclass(Q, "ppp", fatal = F)) {
#		warning("using default quadrature scheme")
		Q <- quadscheme(Q)   
	} else if(!verifyclass(Q, "quad", fatal=F))
		stop("First argument Q should be a quadrature scheme")
#
# Data points
  X <- Q$data
#
# Data and dummy points together 
  P <- union.quad(Q)
#
#
# Interpret the call
want.trend <- !is.null(trend) && !identical.formulae(trend, ~1)
want.inter <- !is.null(interaction) && !is.null(interaction$family)

the.call <- deparse(sys.call())
the.version <- "Spatstat 1.0; mpl.S $Date: 2001/12/11 15:29:17 $"

if(use.gam && exists("is.R") && is.R())
  stop("gam() is not available in R; you have specified use.gam=T")
        
if(!want.trend && !want.inter) {
  # the model is the uniform Poisson process
  # The MPLE of its intensity is the MLE
  npts <- X$n
  volume <- area.owin(X$window) * markspace.integral(X)
  lambda <- npts/volume
  theta <- list("log(lambda)"=log(lambda))
  maxlogpl <- npts * (log(lambda) - 1)
  rslt <- list(
               theta       = theta,
               coef        = theta,
               trend       = NULL,
               interaction = NULL,
               Q           = Q,
               maxlogpl    = maxlogpl,
               internal    = list(),
	       correction  = correction,
               rbord       = rbord,
               call        = the.call,
               version     = the.version)
  class(rslt) <- "ppm"
  return(rslt)
}


################################################################
################ C o m p u t e     d a t a  ####################
################################################################


### Form the weights and the ``response variable''.

W <- w.quad(Q)
Z <- is.data(Q)
Y <- Z/W
n <- length(Y)

MARKS <- marks.quad(Q)    # is NULL for unmarked patterns
	
glmdata <- data.frame(W=W, Y=Y)
SUBSET <- rep(T, n)

zeroes <- attr(W, "zeroes")
if(!is.null(zeroes))
	SUBSET <-  !zeroes

####################### T r e n d ##############################

if(want.trend) {
  # Default explanatory variables for trend
  glmdata <- data.frame(glmdata, x=P$x, y=P$y)
  if(!is.null(MARKS))
    glmdata <- data.frame(glmdata, marks=MARKS)
  # 
  if(!missing(data)) {
# Append `data' to `glmdata'
# First validate 'data'
#   Check it has the right number of rows
    if(nrow(data) != nrow(glmdata))
      stop("Number of rows of \"data\" != number of points in \"Q\"")
#   Check any duplication of reserved names
    name.match <- outer(names(glmdata), names(data), "==")
    if(any(name.match)) {
      is.matched <- apply(name.match, 2, any)
      matched.names <- (names(data))[is.matched]
      if(sum(is.matched) == 1) {
        stop(paste("the variable name \"",
                   matched.names,
                   "\" in \'data\' is a reserved name - see help(mpl)"))
      } else {
        stop(paste("the variable names \"",
                   paste(matched.names, collapse=", "),
                   "\" in \'data\' are reserved names - see help(mpl)"))
      }
    }
#   OK, append 'data'    
    glmdata <- data.frame(glmdata,data)
  }
}

###################### I n t e r a c t i o n ####################

Vnames <- NULL

if(want.inter) {

  verifyclass(interaction, "interact")
  
  # Calculations require a matrix (data) x (data + dummy) indicating equality
  E <- equals.quad(Q)
  
  # Form the matrix of "regression variables" V.
  # The rows of V correspond to the rows of P (quadrature points)
  # while the column(s) of V are the regression variables (log-potentials)

  V <- interaction$family$eval(X, P, E,
                        interaction$pot,
                        interaction$par,
                        correction)

  if(!is.matrix(V))
    stop("interaction evaluator did not return a matrix")

  # Augment data frame by appending the regression variables for interactions.
  #
  # If there are no names provided for the columns of V,
  # call them "Interact.1", "Interact.2", ...

  if(is.null(dimnames(V)[[2]])) {
    # default names
    nc <- ncol(V)
    dimnames(V) <- list(dimnames(V)[[1]], 
      if(nc == 1) "Interaction" else paste("Interact.", 1:nc, sep=""))
  }
  Vnames <- dimnames(V)[[2]]
  
  glmdata <- data.frame(glmdata, V)   

# Keep only those quadrature points for which the
# conditional intensity is nonzero. 

#KEEP  <- apply(V != -Inf, 1, all)
KEEP  <- matrowall(V != -Inf)

SUBSET <- SUBSET & KEEP

if(any(Z & !KEEP)) {
        howmany <- sum(Z & !KEEP)
	warning(paste(howmany, "data point(s) are illegal (zero conditional intensity under the model)"))
#	browser()
}

}
######################################################################
################ F I T    M O D E L  #################################
######################################################################

# Determine the domain of integration for the pseudolikelihood.

if(correction == "border" && !missing(rbord)) {
	bd <- bdist.points(P)
	DOMAIN <- (bd >= rbord)
	SUBSET <- DOMAIN & SUBSET
}

glmdata <- data.frame(glmdata, SUBSET=SUBSET)

#################  F o r m u l a   ##################################


if(!want.trend) trend <- ~1 
trendpart <- paste(as.character(trend), collapse=" ")
rhs <- paste(c(trendpart, Vnames), collapse= "+")
fmla <- paste("Y ", rhs)
fmla <- as.formula(fmla)

################# F i t    i t   ####################################

# Fit the generalized linear/additive model.

if(want.trend && use.gam)
  FIT  <- gam(fmla, family=quasi(link=log, var=mu), weights=W,
              data=glmdata, subset=(SUBSET=="TRUE"), maxit=50)
else
  FIT  <- glm(fmla, family=quasi(link=log, var=mu), weights=W,
              data=glmdata, subset=(SUBSET=="TRUE"), maxit=50)
  
################  I n t e r p r e t    f i t   #######################

# Fitted coefficients

co <- FIT$coef
theta <- if(exists("is.R") && is.R()) NULL else dummy.coef(FIT)

  
# attained value of max log pseudolikelihood
maxlogpl <-  -(deviance(FIT)/2 + sum(log(W[Z & SUBSET])) + sum(Z & SUBSET))

######################################################################
# Clean up & return 

rslt <- list(
             theta        = theta,
             coef         = co,
             trend        = if(want.trend) trend       else NULL,
             interaction  = if(want.inter) interaction else NULL,
             Q            = Q,
             maxlogpl     = maxlogpl, 
             internal     = list(glmfit=FIT, glmdata=glmdata, Vnames=Vnames),
             correction   = correction,
             rbord        = rbord,
             call         = the.call,
             version      = the.version)
class(rslt) <- "ppm"
return(rslt)
}  

#
#
#    multipair.family.S
#
#    $Revision: 1.3 $	$Date: 2001/08/07 08:44:48 $
#
#    Pairwise interaction class for multitype point processes
#    i.e. marked point processes with a finite number of possible types
#
#    multipair.family:      object of class 'isf' defining
#                           pairwise interaction for multitype point processes
#	
# -------------------------------------------------------------------
#	

multipair.family <-
  list(
         name  = "multipair",
         print = function(self) {
                      cat("Multitype pairwise interaction family\n")
         },
         eval  = function(X,U,Equal,pairpot,potpars,correction) {
  #
  # multipair.family$eval
  #
  #  $Revision: 1.3 $  $Date: 2001/08/07 08:44:48 $
  #         
  # This auxiliary function is not meant to be called by the user.
  # It computes the distances between points,
  # evaluates the pair potential and applies edge corrections.
  #
  # Arguments:
  #   X           data point pattern (marked)             'ppp' object
  #   U           points at which to evaluate potential   list(x,y,marks)
  #   Equal       logical matrix X[i] == U[j]             matrix or NULL
  #                      (NB: equality only if marks are equal too)
  #                      NULL means all comparisons are FALSE
  #   pairpot     potential function (see above)          function()
  #   potpars     auxiliary parameters for pairpot        list(......)
  #   correction  edge correction type                    (string)
  #
  # Value:
  #    matrix of values of the total pair potential
  #    induced by the pattern X at each location given in U.
  #    The rows of this matrix correspond to the rows of U (the sample points);
  #    the k columns are the coordinates of the k-dimensional potential.
  #
  # Note:
  # The pair potential function 'pairpot' will be called as
  #    pairpot(M, V1, V2, potpars)
  # where M is a matrix of interpoint distances,
  #       V1 and V2 are vectors of marks for the rows and columns of M resp.
  # It must return a matrix with the same dimensions as M
  # or an array with its first two dimensions the same as the dimensions of M.
  ##########################################################################

# coercion should be unnecessary..
# X <- as.ppp(X)
# U <- as.ppp(U, X$window)   # i.e. X$window is DEFAULT window

x <- X$x
y <- X$y
m <- X$marks

xx <- U$x
yy <- U$y
mm <- U$marks

if(!any(correction == c("periodic", "border", "translate", "none")))
  stop(paste("Unrecognised edge correction \'", correction, "\'", sep=""))
          
#  
# Form the matrix of distances
	
sqdif <- function(u,v) {(u-v)^2}

MX <- outer(x,xx,sqdif)
MY <- outer(y,yy,sqdif)

if(correction=="periodic") {
	if(X$window$type != "rectangle")
          stop("Periodic edge correction can't be applied",
               "in an irregular window")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
	MX1 <- outer(x,xx-wide,sqdif)
	MX2 <- outer(x,xx+wide,sqdif)
	MX <- pmin(MX, MX1, MX2)
	MY1 <- outer(y,yy-high,sqdif)
	MY2 <- outer(y,yy+high,sqdif)
	MY <- pmin(MY, MY1, MY2)
}
M <- sqrt(MX + MY)

# Evaluate the pairwise potential 

POT <- pairpot(M, m, mm, potpars)
if(length(dim(POT)) == 1 || any(dim(POT)[1:2] != dim(M))) {
        whinge <- paste(
           "The pair potential function ",deparse(substitute(pairpot)),
           "must produce a matrix or array with its first two dimensions\n",
           "the same as the dimensions of its input.\n", sep="")
	stop(whinge)
}

# make it a 3D array
if(length(dim(POT))==2)
        POT <- array(POT, dim=c(dim(POT),1), dimnames=NULL)
                          
if(correction == "translate") {
	if(X$window$type != "rectangle")
		stop("sorry, translation correction is not yet implemented",
                     "for irregular windows")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
        DX <- abs(outer(x,xx,"-"))
        DY <- abs(outer(y,yy,"-"))
                                        # translation correction
        edgewt <- wide * high / ((wide - DX) * (high - DY))
        edgewt <- pmin(edgewt, 1000)    # arbitrary stabilisation
        POT <- c(edgewt) * POT
}

# No pair potential term between a point and itself
if(!is.null(Equal))
  POT[Equal] <- 0

# Sum the pairwise potentials 

V <- apply(POT, c(2,3), sum)

return(V)

}
######### end of function $eval                            
)
######### end of list

class(multipair.family) <- "isf"





###########    utilities for this family  #######################


MultiPair.checkmatrix <-
  function(mat, n, name) {
    if(!is.matrix(mat))
      stop(paste(name, "must be a matrix"))
    if(any(dim(mat) != rep(n,2)))
      stop(paste(name, "must be a square matrix,",
                 "of size", n, "x", n))
    isna <- is.na(mat)
    if(any(mat[!isna] <= 0))
      stop(paste("Entries of", name,
                 "must be positive numbers or NA"))
    if(any(isna != t(isna)) ||
       any(mat[!isna] != t(mat)[!isna]))
      stop(paste(name, "must be a symmetric matrix"))
  }

#
#
#    multistrauss.S
#
#    $Revision: 1.7 $	$Date: 2001/11/05 08:20:43 $
#
#    The multitype Strauss process
#
#    MultiStrauss()    create an instance of the multitype Strauss process
#                 [an object of class 'interact']
#	
# -------------------------------------------------------------------
#	

MultiStrauss <- function(types, radii) {
  if(length(types) == 1)
    stop("The \`types\' argument should be a vector of all possible types")
  if(!is.factor(types)) {
    # warning("\`types\' converted to a factor")
    types <- factor(types)
  }
  ct <- levels(types)
  dimnames(radii) <- list(ct, ct)
  out <- 
  list(
         name     = "Multitype Strauss process",
         family    = multipair.family,
         pot      = function(d, tx, tu, par) {
     # arguments:
     # d[i,j] distance between points X[i] and U[j]
     # tx[i]  type (mark) of point X[i]
     # tu[j]  type (mark) of point U[j]
     #
     # get matrix of interaction radii r[ , ]
     r <- par$radii
     #
     # get possible marks and validate
     if(!is.factor(tx) || !is.factor(tu))
	stop("marks of data and dummy points must be factor variables")
     lx <- levels(tx)
     lu <- levels(tu)
     if(length(lx) != length(lu) || any(lx != lu))
	stop("marks of data and dummy points do not have same possible levels")

     # list all UNORDERED pairs of types to be checked
     # (the interaction must be symmetric in type, and scored as such)
     uptri <- (row(r) <= col(r)) & !is.na(r)
     mark1 <- (lx[row(r)])[uptri]
     mark2 <- (lx[col(r)])[uptri]
     vname <- apply(cbind(mark1,mark2), 1, paste, collapse="x")
     vname <- paste("mark", vname, sep="")
     npairs <- length(vname)
     # list all ORDERED pairs of types to be checked
     # (to save writing the same code twice)
     different <- mark1 != mark2
     mark1o <- c(mark1, mark2[different])
     mark2o <- c(mark2, mark1[different])
     nordpairs <- length(mark1o)
     # unordered pair corresponding to each ordered pair
     ucode <- c(1:npairs, (1:npairs)[different])
     #
     # go....
     # assemble the relevant interaction distance for each pair of points
     rxu <- r[ tx, tu ]
     # apply relevant threshold to each pair of points
     str <- (d <= rxu)
     # create logical array for result
     z <- array(F, dim=c(dim(d), npairs),
                dimnames=list(character(0), character(0), vname))
     # assign str[i,j] -> z[i,j,k] where k is relevant interaction code
     for(i in 1:nordpairs) {
       # data points with mark m1
       Xsub <- (tx == mark1o[i])
       # quadrature points with mark m2
       Qsub <- (tu == mark2o[i])
       # assign
       z[Xsub, Qsub, ucode[i]] <- str[Xsub, Qsub]
     }
     return(z)
   },
     #### end of 'pot' function ####
     #       
         par      = list(types=types, radii = radii),
         parnames = c("possible types", "interaction distances"),
         init     = function(self) {
                      r <- self$par$radii
                      nt <- length(self$par$types)
                      MultiPair.checkmatrix(r, nt, "\`radii\'")
                    },
         update = NULL,  # default OK
         print = function(self) {
           print.isf(self$family)
           cat(paste("Interaction:\t", self$name, "\n"))
           cat(paste(length(self$par$types), "types of points\n"))
           cat("Possible types: \n")
           print(self$par$types)
           cat("Interaction radii:\n")
           print(self$par$radii)
           invisible()
         }
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#
#    multistrhard.S
#
#    $Revision: 1.6 $	$Date: 2001/11/05 08:20:43 $
#
#    The multitype Strauss/hardcore process
#
#    MultiStraussHard()
#                 create an instance of the multitype Strauss/ harcore
#                 point process
#                 [an object of class 'interact']
#	
# -------------------------------------------------------------------
#	

MultiStraussHard <- function(types, iradii, hradii) {
  if(length(types) == 1)
    stop("The \`types\' argument should be a vector of all possible types")
  if(!is.factor(types)) {
    # warning("\`types\' converted to a factor")
    types <- factor(types)
  }
  ct <- levels(types)
  dimnames(iradii) <- dimnames(hradii) <- list(ct, ct)
  out <- 
  list(
         name     = "Multitype Strauss Hardcore process",
         family    = multipair.family,
         pot      = function(d, tx, tu, par) {
     # arguments:
     # d[i,j] distance between points X[i] and U[j]
     # tx[i]  type (mark) of point X[i]
     # tu[i]  type (mark) of point U[j]
     #
     # get matrices of interaction radii
     r <- par$iradii
     h <- par$hradii

     # get possible marks and validate
     if(!is.factor(tx) || !is.factor(tu))
	stop("marks of data and dummy points must be factor variables")
     lx <- levels(tx)
     lu <- levels(tu)
     if(length(lx) != length(lu) || any(lx != lu))
	stop("marks of data and dummy points do not have same possible levels")

     # list all UNORDERED pairs of types to be checked
     # (the interaction must be symmetric in type, and scored as such)
     uptri <- (row(r) <= col(r)) & (!is.na(r) | !is.na(h))
     mark1 <- (lx[row(r)])[uptri]
     mark2 <- (lx[col(r)])[uptri]
     vname <- apply(cbind(mark1,mark2), 1, paste, collapse="x")
     vname <- paste("mark", vname, sep="")
     npairs <- length(vname)
     # list all ORDERED pairs of types to be checked
     # (to save writing the same code twice)
     different <- mark1 != mark2
     mark1o <- c(mark1, mark2[different])
     mark2o <- c(mark2, mark1[different])
     nordpairs <- length(mark1o)
     # unordered pair corresponding to each ordered pair
     ucode <- c(1:npairs, (1:npairs)[different])
     #
     # go....
     # apply the relevant interaction distance to each pair of points
     rxu <- r[ tx, tu ]
     str <- (d < rxu)
     str[is.na(str)] <- F
     # and the relevant hard core distance
     hxu <- h[ tx, tu ]
     forbid <- (d < hxu)
     forbid[is.na(forbid)] <- F
     # form the potential 
     value <- ifelse(forbid, -Inf, str)
     # create numeric array for result
     z <- array(0, dim=c(dim(d), npairs),
                dimnames=list(character(0), character(0), vname))
     # assign value[i,j] -> z[i,j,k] where k is relevant interaction code
     for(i in 1:nordpairs) {
       # data points with mark m1
       Xsub <- (tx == mark1o[i])
       # quadrature points with mark m2
       Qsub <- (tu == mark2o[i])
       # assign
       z[Xsub, Qsub, ucode[i]] <- value[Xsub, Qsub]
     }     
     return(z)
     },
     #### end of 'pot' function ####
     #       
         par      = list(types=types, iradii = iradii, hradii = hradii),
         parnames = c("possible types", "interaction distances", "hardcore distances"),
         init     = function(self) {

                      r <- self$par$iradii
                      h <- self$par$hradii
                      nt <- length(self$par$types)

                      MultiPair.checkmatrix(r, nt, "\`iradii\'")
                      MultiPair.checkmatrix(h, nt, "\`hradii\'")
                      
                    },
         update = NULL,  # default OK
         print = function(self) {
           print.isf(self$family)
           cat(paste("Interaction:\t", self$name, "\n"))
           cat(paste(length(self$par$types), "types of points\n"))
           cat("Possible types: \n")
           print(self$par$types)
           cat("Interaction radii:\n")
           print(self$par$iradii)
           cat("Hardcore radii:\n")
           print(self$par$hradii)
           invisible()
         }
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#
#    ord.S
#
#    $Revision: 1.2 $	$Date: 2001/08/07 11:52:17 $
#
#    Ord process with user-supplied potential
#
#    Ord()  create an instance of the Ord process
#                 [an object of class 'interact']
#                 with user-supplied potential
#	
#
# -------------------------------------------------------------------
#	

Ord <- function(pot, name) {
  if(missing(name))
    name <- "Ord process with user-defined potential"
  
  out <- 
  list(
         name     = name,
         family    = ord.family,
         pot      = pot,
         par      = NULL,
         parnames = NULL,
         init     = NULL,
         update   = NULL, 
         print = function(self) {
           cat(paste(self$name, "\n"))
           cat("Potential function:\n")
           print(self$pot)
           invisible()
         }
  )
  class(out) <- "interact"
  return(out)
}
#
#
#    ord.family.S
#
#    $Revision: 1.8 $	$Date: 2001/08/09 10:53:16 $
#
#    The Ord model (family of point process models)
#
#    ord.family:      object of class 'isf' defining Ord model structure
#	
#
# -------------------------------------------------------------------
#	

ord.family <-
  list(
         name  = "ord",
         print = function(self) {
                      cat("Ord model family\n")
         },
         eval  = function(X, U, Equal, pot, pars, ...) {
  #
  # This auxiliary function is not meant to be called by the user.
  # It computes the distances between points,
  # evaluates the pair potential and applies edge corrections.
  #
  # Arguments:
  #   X           data point pattern                      'ppp' object
  #   U           points at which to evaluate potential   list(x,y) suffices
  #   Equal       logical matrix X[i] == U[j]             matrix or NULL
  #   pot         potential function                      function(d, p)
  #   pars        auxiliary parameters for pot            list(......)
  #   ...         IGNORED                             
  #
  # Value:
  #    matrix of values of the potential
  #    induced by the pattern X at each location given in U.
  #    The rows of this matrix correspond to the rows of U (the sample points);
  #    the k columns are the coordinates of the k-dimensional potential.
  #
  # Note:
  # The potential function 'pot' will be called as
  #    pot(M, pars)   where M is a vector of tile areas.
  # It must return a vector of the same length as M
  # or a matrix with number of rows equal to the length of M
  ##########################################################################

nall <- length(U$x)       # number of data + dummy points

# determine which points in the combined list are data points
if(!is.null(Equal))           
  #is.data <- apply(Equal, 2, any)
  is.data <- matcolany(Equal)
else
  is.data <- rep(F, nall)

#############################################################################
# First compute Dirichlet tessellation of data
# and its total potential (which could be vector-valued)
#############################################################################

Wdata <- dirichlet.weights(X)   # sic - these are the tile areas.
Pdata <- pot(Wdata, pars)
summa <- function(P) {
  if(is.vector(P)) sum(P) else matrowsum(P)
}
total.data.potential <- summa(Pdata)

# Initialise V

dimpot <- dim(Pdata)[-1]  # dimension of each value of the potential function
                          # (NULL if scalar)
dimV <- c(nall, dimpot)
V <- array(0, dim=dimV)

rowV <- array(1:nall, dim=dimV)

#################### Next, evaluate V for the data points.  ###############
# For each data point, compute Dirichlet tessellation
# of the data with this point removed.
# Compute difference of total potential.
#############################################################################


Xminus <- X
for(j in seq(X$n)) {
  Xminus$x <- X$x[-j]
  Xminus$y <- X$y[-j]
        #  Dirichlet tessellation of data without point j
  Wminus <- dirichlet.weights(Xminus)
        #  regressor is the difference in total potential
  V[rowV == j] <- total.data.potential - summa(pot(Wminus, pars))
}


#################### Next, evaluate V for the dummy points   ################
# For each dummy point, compute Dirichlet tessellation
# of (data points together with this dummy point) only. 
# Take difference of total potential.
#############################################################################

Xplus <- X

for(j in seq(U$x)[!is.data]) {
		Xplus$x <- c(U$x[j], X$x)
		Xplus$y <- c(U$y[j], X$y)
                     #  compute Dirichlet tessellation (of these points only!)
                Wplus <- dirichlet.weights(Xplus)
                     #  regressor is difference in total potential
                V[rowV == j] <- summa(pot(Wplus)) - total.data.potential
}

return(V)

} ######### end of function $eval                            

) ######### end of list

class(ord.family) <- "isf"
#
#
#    ordthresh.S
#
#    $Revision: 1.1 $	$Date: 2000/07/11 10:49:01 $
#
#    Ord process with threshold potential
#
#    OrdThresh()  create an instance of the Ord process
#                 [an object of class 'interact']
#                 with threshold potential
#	
#
# -------------------------------------------------------------------
#	

OrdThresh <- function(r) {
  out <- 
  list(
         name     = "Ord process with threshold potential",
         family    = ord.family,
         pot      = function(d, par) {
                         ifelse(d <= par$r, 1, 0)
                    },
         par      = list(r = r),
         parnames = "threshold distance",
         init     = function(self) {
                      r <- self$par$r
                      if(!is.numeric(r) || length(r) != 1 || r <= 0)
                       stop("threshold distance r must be a positive number")
                    },
         update = NULL,  # default OK
         print = NULL    # default OK
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#
#    pairpiece.S
#
#    $Revision: 1.1 $	$Date: 2000/07/11 10:49:01 $
#
#    A pairwise interaction process with piecewise constant potential
#
#    PairPiece()   create an instance of the process
#                 [an object of class 'interact']
#	
#
# -------------------------------------------------------------------
#	

PairPiece <- function(r) {
  out <- 
  list(
         name     = "Piecewise constant pairwise interaction process",
         family    = pairwise.family,
         pot      = function(d, par) {
                         r <- par$r
                         nr <- length(r)
                         out <- array(F, dim=c(dim(d), nr))
                         out[,,1] <-  ifelse(d <= r[1], 1, 0)
                         for(i in 2:nr) 
                           out[,,i] <- ifelse((d > r[i-1]) & (d <= r[i]), 1, 0)
                         out
                    },
         par      = list(r = r),
         parnames = "interaction thresholds",
         init     = function(self) {
                      r <- self$par$r
                      if(!is.numeric(r) || length(r) < 2 || !all(r > 0))
                       stop("interaction thresholds r must be positive numbers")
                      if(!all(diff(r) > 0))
                        stop("interaction thresholds r must be strictly increasing")
                    },
         update = NULL,  # default OK
         print = NULL    # default OK
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#
#    pairsat.family.S
#
#    $Revision: 1.9 $	$Date: 2001/08/09 11:35:41 $
#
#    The saturated pairwise interaction family of point process models
#
#    (an extension of Geyer's saturation process to all pairwise interactions)
#
#    pairsat.family:         object of class 'isf'
#                     defining saturated pairwise interaction
#	
#
# -------------------------------------------------------------------
#	

pairsat.family <-
  list(
         name  = "saturated pairwise",
         print = function(self) {
                      cat("Saturated pairwise interaction family\n")
         },
         eval  = function(X,U,Equal,pairpot,potpars,correction) {
  #
  # This auxiliary function is not meant to be called by the user.
  # It computes the distances between points,
  # evaluates the pair potential and applies edge corrections.
  #
  # Arguments:
  #   X           data point pattern                      'ppp' object
  #   U           points at which to evaluate potential   list(x,y) suffices
  #   Equal       logical matrix X[i] == U[j]             matrix or NULL
  #   pairpot     potential function (see above)          function(d, p)
  #   potpars     auxiliary parameters for pairpot        list(......)
  #   correction  edge correction type                    (string)
  #
  #   Note the Geyer saturation threshold must be given in 'potpars$saturate'
  #
  # Value:
  #    matrix of values of the total pair potential
  #    induced by the pattern X at each location given in U.
  #    The rows of this matrix correspond to the rows of U (the sample points);
  #    the k columns are the coordinates of the k-dimensional potential.
  #
  # Note:
  # The pair potential function 'pairpot' will be called as
  #    pairpot(M, potpars)   where M is a matrix of interpoint distances.
  # It must return a matrix with the same dimensions as M
  # or an array with its first two dimensions the same as the dimensions of M.
  ##########################################################################

# coercion should be unnecessary, but this is useful for debugging
X <- as.ppp(X)
U <- as.ppp(U, X$window)   # i.e. X$window is DEFAULT window

# extract coordinates
x <- X$x
y <- X$y
xx <- U$x
yy <- U$y

ndata <- length(x)
nall <- length(xx)
ndummy <- nall - ndata

# saturation parameter
saturate <- potpars$saturate

# validate 'correction' argument
if(!any(correction == c("periodic", "border", "translate", "none")))
  stop(paste("Unrecognised edge correction \'", correction, "\'", sep=""))
          
#  
# Form the matrix of distances
	
sqdif <- function(u,v) {(u-v)^2}

MX <- outer(x,xx,sqdif)
MY <- outer(y,yy,sqdif)

if(correction=="periodic") {
	if(X$window$type != "rectangle")
		stop("Periodic edge correction can't be applied",
                     "in an irregular window")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
	MX1 <- outer(x,xx-wide,sqdif)
	MX2 <- outer(x,xx+wide,sqdif)
	MX <- pmin(MX, MX1, MX2)
	MY1 <- outer(y,yy-high,sqdif)
	MY2 <- outer(y,yy+high,sqdif)
	MY <- pmin(MY, MY1, MY2)
}
M <- sqrt(MX + MY)

# Evaluate the pairwise potential 

POT <- pairpot(M, potpars)
if(length(dim(POT)) == 1 || any(dim(POT)[1:2] != dim(M))) {
        whinge <- paste(
           "The pair potential function ",deparse(substitute(pairpot)),
           "must produce a matrix or array with its first two dimensions\n",
           "the same as the dimensions of its input.\n", sep="")
	stop(whinge)
}

# make it a 3D array
if(length(dim(POT))==2)
        POT <- array(POT, dim=c(dim(POT),1), dimnames=NULL)
                          
if(correction == "translate") {
	if(X$window$type != "rectangle")
		stop("sorry, translation correction is not yet implemented",
                     "for irregular windows")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
        DX <- abs(outer(x,xx,"-"))
        DY <- abs(outer(y,yy,"-"))
                                        # translation correction
        edgewt <- wide * high / ((wide - DX) * (high - DY))
        edgewt <- pmin(edgewt, 1000)    # arbitrary stabilisation
        POT <- c(edgewt) * POT
}

# No pair potential term between a point and itself
if(!is.null(Equal))
  POT[Equal] <- 0

# Sum the pairwise potentials (sum over the data points)
   V <- apply(POT, c(2,3), sum)

if(is.null(saturate))
  return(V)

#################################################################
################## saturation part ##############################
#################################################################

#
# (a) compute SATURATED potential sums
V.sat <- array(pmin(saturate, V), dim=dim(V))

#
# (b) compute effect of addition/deletion of dummy/data point j
# on the UNSATURATED potential sum of each data point i
#
# Identify data points
     # is.data <- apply(Equal, 2, any)
is.data <- matcolany(Equal)              # logical vector corresp. to V

# Extract potential sums for data points only
V.data <- V[is.data, , drop=F]

# replicate them so that V.dat.rep[i,j,k] = V.data[i, k]
V.dat.rep <- aperm(array(V.data, dim=c(dim(V.data), nall)), c(1,3,2))

# make a logical array   col.is.data[i,j,k] = is.data[j]
dip <- dim(POT)
mat <- matrix(, nrow=dip[1], ncol=dip[2])
izdat <- is.data[col(mat)] 
col.is.data <- array(izdat, dim=dip) # automatically replicates

# compute value of unsaturated potential sum for each data point i
# obtained after addition/deletion of each dummy/data point j
                                  
V.after <- V.dat.rep + ifelse(col.is.data, -POT, POT)
#
#
# (c) difference of SATURATED potential sums for each data point i
# before & after increment/decrement of each dummy/data point j
#
# saturated values after increment/decrement
V.after.sat <- array(pmin(saturate, V.after), dim=dim(V.after))
# saturated values before
V.dat.rep.sat <- array(pmin(saturate, V.dat.rep), dim=dim(V.dat.rep))
# difference
V.delta <- V.after.sat - V.dat.rep.sat
V.delta <- ifelse(col.is.data, V.delta, -V.delta)
#
# (d) Sum (c) over all data points i
V.delta.sum <- apply(V.delta, c(2,3), sum)
#
# (e) Result
V <- V.sat + V.delta.sum
  
return(V)

}     ######### end of function $eval                            

)     ######### end of list

class(pairsat.family) <- "isf"
#
#
#    pairwise.S
#
#    $Revision: 1.2 $	$Date: 2001/08/07 11:52:17 $
#
#    Pairwise()    create a user-defined pairwise interaction process
#                 [an object of class 'interact']
#	
# -------------------------------------------------------------------
#	

Pairwise <- function(pot, name) {
  if(missing(name))
    name <- "user-defined pairwise interaction process"

  out <- 
  list(
         name     = name,
         family   = pairwise.family,
         pot      = pot,
         par      = NULL,
         parnames = NULL,
         init     = NULL,
         update   = NULL,  
         print = function(self) {
           cat(paste(self$name, "\n"))
           cat("Potential function:\n")
           print(self$pot)
           invisible()
         }
  )
  class(out) <- "interact"
  return(out)
}
#
#
#    pairwise.family.S
#
#    $Revision: 1.4 $	$Date: 2001/08/07 08:44:48 $
#
#    The pairwise interaction family of point process models
#
#    pairwise.family:      object of class 'isf' defining pairwise interaction
#	
#
# -------------------------------------------------------------------
#	

pairwise.family <-
  list(
         name  = "pairwise",
         print = function(self) {
                      cat("Pairwise interaction family\n")
         },
         eval  = function(X,U,Equal,pairpot,potpars,correction) {
  #
  # This auxiliary function is not meant to be called by the user.
  # It computes the distances between points,
  # evaluates the pair potential and applies edge corrections.
  #
  # Arguments:
  #   X           data point pattern                      'ppp' object
  #   U           points at which to evaluate potential   list(x,y) suffices
  #   Equal        logical matrix X[i] == U[j]             matrix or NULL
  #   pairpot     potential function (see above)          function(d, p)
  #   potpars     auxiliary parameters for pairpot        list(......)
  #   correction  edge correction type                    (string)
  #
  # Value:
  #    matrix of values of the total pair potential
  #    induced by the pattern X at each location given in U.
  #    The rows of this matrix correspond to the rows of U (the sample points);
  #    the k columns are the coordinates of the k-dimensional potential.
  #
  # Note:
  # The pair potential function 'pairpot' will be called as
  #    pairpot(M, potpars)   where M is a matrix of interpoint distances.
  # It must return a matrix with the same dimensions as M
  # or an array with its first two dimensions the same as the dimensions of M.
  ##########################################################################

# coercion should be unnecessary..
# X <- as.ppp(X)
# U <- as.ppp(U, X$window)   # i.e. X$window is DEFAULT window

x <- X$x
y <- X$y
xx <- U$x
yy <- U$y
  
if(!any(correction == c("periodic", "border", "translate", "none")))
  stop(paste("Unrecognised edge correction \'", correction, "\'", sep=""))
          
#  
# Form the matrix of distances
	
sqdif <- function(u,v) {(u-v)^2}

MX <- outer(x,xx,sqdif)
MY <- outer(y,yy,sqdif)

if(correction=="periodic") {
	if(X$window$type != "rectangle")
		stop("Periodic edge correction can't be applied",
                     "in an irregular window")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
	MX1 <- outer(x,xx-wide,sqdif)
	MX2 <- outer(x,xx+wide,sqdif)
	MX <- pmin(MX, MX1, MX2)
	MY1 <- outer(y,yy-high,sqdif)
	MY2 <- outer(y,yy+high,sqdif)
	MY <- pmin(MY, MY1, MY2)
}
M <- sqrt(MX + MY)

# Evaluate the pairwise potential 

POT <- pairpot(M, potpars)
if(length(dim(POT)) == 1 || any(dim(POT)[1:2] != dim(M))) {
        whinge <- paste(
           "The pair potential function ",deparse(substitute(pairpot)),
           "must produce a matrix or array with its first two dimensions\n",
           "the same as the dimensions of its input.\n", sep="")
	stop(whinge)
}

# make it a 3D array
if(length(dim(POT))==2)
        POT <- array(POT, dim=c(dim(POT),1), dimnames=NULL)
                          
if(correction == "translate") {
	if(X$window$type != "rectangle")
		stop("sorry, translation correction is not yet implemented",
                     "for irregular windows")
	wide <- diff(X$window$xrange)
	high <- diff(X$window$yrange)
        DX <- abs(outer(x,xx,"-"))
        DY <- abs(outer(y,yy,"-"))
                                        # translation correction
        edgewt <- wide * high / ((wide - DX) * (high - DY))
        edgewt <- pmin(edgewt, 1000)    # arbitrary stabilisation
        POT <- c(edgewt) * POT
}

# No pair potential term between a point and itself
if(!is.null(Equal))
  POT[Equal] <- 0

# Sum the pairwise potentials 

V <- apply(POT, c(2,3), sum)

return(V)

}
######### end of function $eval                            
)
######### end of list

class(pairwise.family) <- "isf"
#
#   pcf.R
#
#   $Revision: 1.2 $   $Date: 2001/11/23 06:25:02 $
#
#
#   calculate pair correlation function
#   from estimate of K or Kcross
#
#

"pcf" <-
function(X, ..., method="c") { 
	require(modreg)

	if(verifyclass(X, "ppp", fatal=FALSE))
        # point pattern - estimate K and continue
		X <- Kest(X)
        else if(verifyclass(X, "fasp", fatal=FALSE)) {
          # function array - go to work on each function
          n <- length(X$fns)
          for(i in 1:n) {
            Xi <- X$fns[[i]]
            if(!is.data.frame(Xi))
              stop("Internal error - an entry in the function array is not a data frame")
            X$fns[[i]] <- pcf(Xi, ..., method=method)
          }
	return(X)
        }
	else if(!is.data.frame(X) || is.null(X$r) || is.null(X$border))
		stop("X should be either a point pattern or the value returned by Kest() or Kcross() or alltypes(..., \"K\")")
        
	# remove NA's
	ok <- !is.na(X$border) 
	X <- X[ok, ]
	r <- X$r
	K <- X$border
	switch(method,
		a = { 
			ss <- smooth.spline(r, K, ...)
			dK <- predict.smooth.spline(ss, r, deriv=1)$y
			g <- dK/(2 * pi * r)
		},
		b = {
			y <- K/(2 * pi * r)
			y[is.nan(y)] <- 0
			ss <- smooth.spline(r, y, ...)
			dy <- predict.smooth.spline(ss, r, deriv=1)$y
			g <- dy + y/r
		},
		c = {
			z <- K/(pi * r^2)
			z[is.nan(z)] <- 1
			ss <- smooth.spline(r, z, ...)
			dz <- predict.smooth.spline(ss, r, deriv=1)$y
			g <- (r/2) * dz + z
		},
		stop(paste("unrecognised method \"", method, "\""))
	)
	X$pcf <- g
	return(X)
}
#
#   plot.fasp.R
#
#   $Revision: 1.5 $   $Date: 2002/01/18 06:27:22 $
#
plot.fasp <- function(x,formula=NULL,subset=NULL,lty=NULL,
                      col=NULL,title=NULL,...) {

# If the formula is null, look for a default formula in x:
	if(is.null(formula)) {
		if(is.null(x$default.formula))
			stop("No formula supplied.\n")
		formula <- x$default.formula
	}
# The formula should be a single formula or a list of formulae.
# If it is a single formula, wrap it up in a list so that all
# formula arguments can be treated consistently.
        if(!is.list(formula)) formula <- list(formula)

# Check on the length of the formula argument.
nf <- length(formula)
if(nf > 1) {
	if(nf != length(x$fns))
		stop("Wrong number of entries in formula argument.\n")
	mfor <- T
} else mfor <- F

# Check on the length of the subset argument.
ns <- length(subset)
if(ns > 1) {
	if(ns != length(x$fns))
		stop("Wrong number of entries in subset argument.\n")
	msub <- T
} else msub <- F

# Set up the array of plotting regions.
	mfrow.save <- par("mfrow")
	oma.save   <- par("oma")
	on.exit(par(mfrow=mfrow.save,oma=oma.save))
	which <- x$which
	m  <- nrow(which)
	n  <- ncol(which)
	nm <- n * m
	par(mfrow=c(m,n))
        # decide whether panels require subtitles
        subtit <- (nm > 1) || !(is.null(x$titles[[1]]) || x$titles[[1]] == "")
	if(nm>1) par(oma=c(0,3,4,0))
        else if(subtit) par(oma=c(3,3,4,0))
        
# Run through the components of the structure x, plotting each
# in the appropriate region, according to the formula.
	k <- 0
	for(i in 1:m) {
		for(j in 1:n) {
# Now do the actual plotting.
			k <- which[i,j]
			if(is.na(k)) plot(0,0,type='n',xlim=c(0,1),
					  ylim=c(0,1),axes=F,xlab='',ylab='')
			else {
				dat.loc <- x$fns[[k]]
				sij <- if(msub) subset[[k]] else subset
				fij <- if(mfor) formula[[k]] else formula[[1]]
				conspire(dat.loc,fij,sij,lty,col)

# Add the (sub)title of each plot.
				if(!is.null(x$titles[[k]]))
					title(main=x$titles[[k]])
			}
		}
	}

# Add an overall title.
	if(!is.null(title)) overall <- title
	else if(!is.null(x$title)) overall <- x$title
	else {
		if(nm > 1)
			overall <- "Array of diagnostic functions"
		else
			overall <- "Diagnostic function"
		if(is.null(x$dataname)) overall <- paste(overall,".",sep="")
		else overall <- paste(overall," for ",x$dataname,".",sep="")
		
	}
	if(nm > 1 || subtit)
          mtext(side=3,outer=T,line=1,text=overall,cex=1.2)
	else title(main=overall)
	invisible()
}
#
#	plot.owin.S
#
#	The 'plot' method for observation windows (class "owin")
#
#	$Revision: 1.4 $	$Date: 2002/01/21 02:24:37 $
#
#
#

plot.owin <- function(x, main, ..., box=T, edge=0.04)
{
#
# Function plot.owin.  A method for plot.
#
# argument must be called 'x' for compatibility with plot()
        W <- x
        if(missing(main))
          main <- deparse(substitute(x))
# no, this cannot be inserted in the argument list!!!
	verifyclass(W, "owin")

#########        
        x <- W$xrange
        y <- W$yrange

# see 'help(par)' under 'xaxs': default is xaxs="r";
# we seize full control by setting xaxs="i",
# but mimic the extra 4% space allocated when xaxs="r"
        
        blowup <- function(v, s) { mean(v) + s * (v - mean(v)) }
        xlim <- blowup(x, 1+edge)
        ylim <- blowup(y, 1+edge)

# Constrain X scale = Y scale        
# I can't believe I have to program this
        pin <- par("pin")  #physical size of plot region
        xscale <- pin[1]/diff(xlim)
        yscale <- pin[2]/diff(ylim)
        sc <- min(xscale,yscale)
        if(xscale > sc) xlim <- blowup(xlim, xscale/sc)
        if(yscale > sc) ylim <- blowup(ylim, yscale/sc) 

# Commit scales
        plot(x, y, xlim=xlim, ylim=ylim, type="n", ...,
             main=main, axes=F, xlab="", ylab="", xaxs="i", yaxs="i")

# Draw window

        switch(W$type,
               rectangle = {
               },
               polygonal = {
                 p <- W$bdry
                 if(exists("is.R") && is.R()) {
                   for(i in seq(p))
                     polygon(p[[i]], ...)
                 } else {
                   for(i in seq(p))
                     polygon(p[[i]], density=0, ...)
                 }
               },
               mask = {
                 image(W$xcol, W$yrow, !t(W$m), add=T, ...)
               },
               stop("Don't know how to plot image type", W$type)
               )

# Draw surrounding box
        if(box)
          segments(x[c(1,2,2,1)],
                   y[c(1,1,2,2)],
                   x[c(2,2,1,1)],
                   y[c(1,2,2,1)])

	invisible()
}





#
#    plot.ppm.S
#
#    $Revision: 1.9 $    $Date: 2002/01/18 06:35:13 $
#
#    plot.ppm()
#         Plot a point process model fitted by mpl().
#        
#
#
plot.ppm <- function(x, nx = 40, ny = 40, 
		     superimpose = T,
                     trend = T, cif = T, pause = T,
                     how=c("persp","image", "contour"),  ...)
{
  
        model <- x
        
#       Plot a point process model fitted by mpl().
#
        verifyclass(model, "ppm")
#
#       find out what kind of model it is
#
        stationary <- is.stationary.ppm(model)
        poisson <- is.poisson.ppm(model)
        markeddata <- is.marked(model$Q$data)
        funnywindow <- (model$Q$data$window$type != "rectangle")

        data <- model$Q$data
        
        if(markeddata) {
          mrks <- data$marks
          if(!is.factor(mrks))
            stop("Marks are not a factor -- I don't know how to plot this")
          mrkvals <- levels(mrks)
        }
#
# Interpret options
#
#        
#                        plotting style
#        
        howmat <- outer(how, c("persp","image", "contour"), "==")
        howmatch <- apply(howmat, 1, any)
        if(any(!howmatch))
          stop(paste("unrecognised option", how[!howmatch]))

#                       whether to plot trend, cif
        
        if(!missing(trend) && !missing(cif) && !trend && !cif) {
          cat("Nothing plotted - both \'trend\' and \'cif\' are F\n")
          return(invisible(NULL))
        }
#                       suppress uninteresting plots
#                       unless explicitly instructed otherwise
        if(missing(trend))
          trend <- !stationary
	if(missing(cif))
          cif <- !poisson
        
        if(!trend && !cif) {
          cat("Nothing plotted -- all plots selected are flat surfaces.\n")
          return(invisible())
        }

#
#        
#
# Do the plotting --- first the trend, then the full cif
# (if both are requested).
#
	types <- c("trend","cif")[c(trend,cif)]
	for(ttt in types) {
          # compute the predictions
          xyz <-  predict.ppm(model, nx = nx, ny = ny, type=ttt)
          # 
          if(funnywindow) {
            # determine which prediction points are inside the window
            xy <- expand.grid(x=xyz$x, y=xyz$y)
            inside <- inside.owin(xy$x, xy$y, model$Q$data$window)
          }
          if(!markeddata) {
            # can use xyz directly
            #
            # first set predictions outside window to NA
            if(funnywindow)
              xyz$z[!inside] <- NA
            #
            for(style in how) {
              switch(style, 
                persp = persp(xyz, ...),
                image = {
                  plot(data$window, main="")
                  image(xyz, add=T, ...)
                  if(superimpose) plot(data, add=T)
                },
                contour = {
                  plot(data$window, main="")
                  contour(xyz, add=T, ...)
                  if(superimpose) plot(data, add=T)
                },
                {stop(paste("Unrecognised plot style", style))}
              )
              lastone <- (style == how[length(how)] &&
                          ttt == types[length(types)] )
              if(pause && !lastone) {
                cat('Next plot? ')
                readline()
              }
            }
          } else {
            # marked points - slice the prediction at each level
            slice <- xyz
            for(level in mrkvals) {
              slice$z <- xyz$z[,,level]
              # set predictions outside window to NA
              if(funnywindow)
                slice$z[!inside] <- NA
              #
              main <- paste("mark =",level)
              for(style in how) {
                switch(style, 
                  persp = {
                    persp(slice, ...)
                    title(main=main)
                  },
                  image = {
                    plot(data$window, main=main)
                    image(slice, add=T, ...)
                    if(superimpose) plot(data[data$marks == level], add=T)
                  },
                  contour = {
                    plot(data$window, main=main)
                    contour(slice, add=T, ...)
                    if(superimpose) plot(data[data$marks == level], add=T)
                  },
                  {
                    stop(paste("Unrecognised plot style", style))
                  }
                )
              lastone <- (style == how[length(how)] &&
                          ttt == types[length(types)] &&
                          level == mrkvals[length(mrkvals)])
              if(pause && !lastone) {
                cat('Next plot? ')
                readline()
              }
              }
            }
          }
        }
        invisible()
}      

#
#	plot.ppp.S
#
#	$Revision: 1.7 $	$Date: 2002/01/18 06:39:34 $
#
#
#--------------------------------------------------------------------------

plot.ppp <- function(x, main, ..., chars, use.marks=T, add=F)
{
#
# Function plot.ppp.
# A plot() method for the class 'ppp'
#
	if(missing(main))
		main <- deparse(substitute(x))

        x <- as.ppp(x)

        if(!add)
          plot.owin(x$window, ..., main=main)

        if(x$n == 0)
          return(invisible())
        
	if(is.null(x$marks) || !use.marks) {
		points(x$x, x$y, ...)
                return(invisible())
        } else {
	        um <- if(is.factor(x$marks))
                         levels(x$marks)
                      else
                         sort(unique(x$marks))

                if(missing(chars))
                  chars <- seq(um)
                
                for(i in seq(um)) {
                  relevant <- (x$marks == um[i])
                  if(any(relevant)) 
                    points(x$x[relevant], x$y[relevant], pch = chars[i])
                }
                names(chars) <- um
                return(chars)
	}
}
#
#
#    poisson.S
#
#    $Revision: 1.3 $	$Date: 2001/08/07 11:52:17 $
#
#    The Poisson process
#
#    Poisson()    create an object of class 'interact' describing
#                 the (null) interpoint interaction structure
#                 of the Poisson process.
#	
#
# -------------------------------------------------------------------
#	

Poisson <- function() {
  out <- 
  list(
         name     = "Poisson process",
         family   = NULL,
         pot      = NULL,
         par      = NULL,
         parnames = NULL,
         init     = function(...) { },
         update   = function(...) { },
         print    = function(self) {
           cat("Poisson process\n")
           invisible()
         }
  )
  class(out) <- "interact"
  return(out)
}
#
#	ppm.S
#
#	Class 'ppm' representing fitted point process models.
#
#
#	$Revision: 1.17 $	$Date: 2002/01/18 06:56:34 $
#
#       An object of class 'ppm' contains the following:
#
#            $coef             vector of fitted regular parameters
#                              as given by coef(glm(....))
#
#            $theta            vector of fitted regular parameters
#                              as given by dummy.coef(glm(....))
#
#            $trend            the trend formula
#                              or NULL 
#
#            $interaction      the interaction family 
#                              (an object of class 'interact') or NULL
#
#            $Q                the quadrature scheme used
#
#            $maxlogpl         the maximised value of log pseudolikelihood
#
#            $internal         list of internal calculation results
#
#            $correction       name of edge correction method used
#            $rbord            erosion distance for border correction (or NULL)
#
#            $the.call         the originating call to mpl()
#
#            $the.version      version of mpl() which yielded the fit
#
#
#------------------------------------------------------------------------

print.ppm <- function(x, ...) {
	verifyclass(x, "ppm")

        notrend <- no.trend.ppm(x)
	stationary <- is.stationary.ppm(x)
	poisson <- is.poisson.ppm(x)
        markeddata <- is.marked(x$Q$data)
        
        markedpoisson <- poisson && markeddata
        # i.e. if the data are marked and mpl() is called with
        # a Poisson interaction, the fitted model is a marked Poisson process.

        # names of interaction variables if any
        Vnames <- x$internal$Vnames

        # ----------- Print model type -------------------
        
	cat(paste(
		if(stationary) "Stationary " else "Nonstationary ",
                if(markedpoisson) "marked " else "", 
                if(poisson) "Poisson process" else x$interaction$name,
                "\n", sep=""))

        if(markeddata) {
          mrk <- x$Q$data$marks
          if(is.factor(mrk)) {
            cat("Possible marks: \n")
            print(factor(levels(mrk)))
          }
        }

        if(exists("is.R") && is.R()) 
          theta <- x$coef # result of coef(glm(...))
        else
          theta <- x$theta # result of dummy.coef(glm(....))
          
        # ----- trivial case: uniform Poisson --------
        
	if(poisson && notrend) {
          lambda <- exp(theta[[1]])
          cat(paste("Uniform intensity", lambda))
          if(markeddata) {
            cat("(for each mark level).\n")
            m <- markspace.integral(x$Q$data)
            cat(paste("Total point intensity", lambda * m))
          } 
          cat("\n")
          return(invisible(NULL))
	}

        # ----- trend --------------------------

        cat(paste("\n\n-----------",
                  if(poisson) "Intensity" else "First order term",
                  "-------------\n\n"))

	if(!notrend) {
		cat("Trend formula:\n")
		print(x$trend)
        }

        # process is at least one of: marked, nonstationary, non-poisson
        
        if(stationary) {
          if(!markeddata) {
            # it can't be poisson
            cat("First order term:\n")
            beta <- exp(theta[[1]])
            print(unlist(list(beta=beta)))
          } else {
            # marked and stationary
            cat( if(poisson)
                "Fitted intensities:\n"
                 else
                "Fitted first order terms:\n")
            lev <- factor(levels(mrk))
            betas <- predict(x, newdata=data.frame(marks=lev), type="trend")
            names(betas) <- paste("beta_", as.character(lev), sep="")
            print(betas)
          }
        } else {
          # not stationary 
          cat("Fitted coefficients for trend formula:\n")
          # extract trend terms without trying to understand them much
          if(is.null(Vnames)) 
            trendbits <- theta
          else {
            agree <- outer(names(theta), Vnames, "==")
            whichbits <- apply(!agree, 1, all)
            trendbits <- theta[whichbits]
          }
          # decide whether there are 'labels within labels'
          unlabelled <- unlist(lapply(trendbits,
                                    function(x) { is.null(names(x)) } ))
          if(all(unlabelled))
            print(unlist(trendbits))
          else 
            for(i in seq(trendbits))
              if(unlabelled[i])
                print(unlist(trendbits[i]))
              else
                print(trendbits[[i]])
        }
        
        # ---- Interaction ----------------------------

	if(!poisson) {
          cat("\n\n-----------Interaction------------------\n\n")
          print(x$interaction)
        
          cat("Fitted interaction terms:\n")
          print(exp(unlist(theta[x$internal$Vnames])))
        }

        cat("\n\n----------- gory details -----\n")
        
        cat("\nFitted regular parameters (theta): \n")
        print(theta)

        cat("\nFitted exp(theta): \n")
        print(exp(unlist(theta)))

	invisible(NULL)
}

no.trend.ppm <- function(x) {
  # the result is T only when the trend is ~1
  verifyclass(x, "ppm")
  length(termsinformula(x$trend)) == 0
}

is.stationary.ppm <- function(x) {
  verifyclass(x, "ppm")
  turms <- termsinformula(x$trend)
  (length(turms) == 0) ||
     ((length(turms) == 1) && turms[1] == "marks")
}

is.poisson.ppm <- function(x) {
  verifyclass(x, "ppm")
  is.null(x$interaction) || is.null(x$interaction$pot)
}

  
#
#	ppp.S
#
#	A class 'ppp' to define point patterns
#	observed in arbitrary windows in two dimensions.
#
#	$Revision: 4.5 $	$Date: 2002/01/18 06:53:42 $
#
#	A point pattern contains the following entries:	
#
#		$window:	an object of class 'owin'
#				defining the observation window
#
#		$n:	the number of points (for efficiency)
#	
#		$x:	
#		$y:	vectors of length n giving the Cartesian
#			coordinates of the points.
#
#	It may also contain the entry:	
#
#		$marks:	a vector of length n
#			whose entries are interpreted as the
#			'marks' attached to the corresponding points.	
#	
#--------------------------------------------------------------------------
ppp <- function(x, y, ..., window, marks ) {
	# Constructs an object of class 'ppp'
	#
        if(!missing(window))
          verifyclass(window, "owin")
        else
          window <- owin(...)
          
	n <- length(x)
	if(length(y) != n)
		stop("coordinate vectors x and y are not of equal length")
	pp <- list(window=window, n=n, x=x, y=y)
	if(!missing(marks) && !is.null(marks)) {
                if(is.matrix(marks) || is.data.frame(marks))
                  stop(paste("Attempted to create point pattern with",
                             ncol(marks), "columns of mark data;",
                             "multidimensional marks are not yet implemented"))
		if(length(marks) != n)
			stop("length of marks vector != length of x and y")
		pp$marks <- marks
	}
	class(pp) <- "ppp"
	pp
}
#
#--------------------------------------------------------------------------
#

as.ppp <- function(X, W = NULL) {
	# tries to coerce data X to a point pattern
	# X may be:
	#	1. an object of class 'ppp'
	#	2. a structure with entries x, y, xl, xu, yl, yu
	#	3. a two-column matrix
	#	4. a structure with entries x, y
        #       5. a quadrature scheme (object of class 'quad')
	# In cases 3 and 4, we need the second argument W
	# which is coerced to an object of class 'owin' by the 
	# function "as.owin" in window.S
        # In cases 2 and 4, if X also has an entry X$marks
        # then this will be interpreted as the marks vector for the pattern.
	#
	if(verifyclass(X, "ppp", fatal=F))
		return(X)
        else if(verifyclass(X, "quad", fatal=F))
                return(union.quad(X))
	else if(checkfields(X, 	c("x", "y", "xl", "xu", "yl", "yu"))) {
		xrange <- c(X$xl, X$xu)
		yrange <- c(X$yl, X$yu)
		if(is.null(X$marks))
			Z <- ppp(X$x, X$y, xrange, yrange)
		else
			Z <- ppp(X$x, X$y, xrange, yrange, 
				marks=X$marks)
		return(Z)
	} else if(is.matrix(X) && is.numeric(X)) {
		if(is.null(W))
			stop("x,y coords given but no window specified")
		win <- as.owin(W)
		Z <- ppp(X[,1], X[,2], window = win)
		return(Z)
	} else if(checkfields(X, c("x", "y"))) {
		if(is.null(W))
			stop("x,y coords given but no window specified")
		win <- as.owin(W)
		if(is.null(X$marks))
                  Z <- ppp(X$x, X$y, window=win)
                else
                  Z <- ppp(X$x, X$y, window=win, marks=X$marks)
                return(Z)
	} else
		stop("Can't interpret X as a point pattern")
}

# --------------------------------------------------------------

"[.ppp" <-
"subset.ppp" <-
  function(x, subset, window, drop, ...) {

        verifyclass(x, "ppp")

        trim <- !missing(window)
        thin <- !missing(subset)
        if(!thin && !trim)
          stop("Please specify a subset (to thin the pattern) or a window (to trim it)")

        # thin first, according to 'subset'
        if(!thin)
          Y <- x
        else
          Y <- ppp(x$x[subset],
                   x$y[subset],
                   window=x$window,
                   marks=if(is.null(x$marks)) NULL else x$marks[subset])

        # now trim to window 
        if(trim) {
          ok <- inside.owin(Y$x, Y$y, window)
          Y <- ppp(Y$x[ok], Y$y[ok],
                   window=window,  # SIC
                   marks=if(is.null(Y$marks)) NULL else Y$marks[ok])
        }
        
        return(Y)
}

# ------------------------------------------------------------------
#
#
scanpp <- function(filename, window, header=T, dir="", multitype=F) {
  filename <- paste(dir, filename, sep="")
  df <- read.table(filename, header=header)
  if(header) {
    x <- df$x
    y <- df$y
    colnames <- dimnames(df)[[2]]
    xycolumns <- match(colnames, c("x","y"), 0)
  } else {
    # assume x, y given in columns 1, 2 respectively
    x <- df[,1]
    y <- df[,2]
    xycolumns <- c(1,2)
  }
  if(ncol(df) == 2) 
      X <- ppp(x, y, window=window)
  else {
    marks <- df[ , -xycolumns]
    if(multitype) 
      marks <- factor(marks)
    X <- ppp(x, y, window=window, marks = marks)
  }
  X
}


"superimpose" <-
  function(...)
{
  # superimpose any number of point patterns
  # ASSUMED TO BE IN THE SAME WINDOW
  # WITH THE SAME MARK SPACE if relevant
  
  # concatenate lists of (x,y) coordinates
  XY <- concatxy(...)
  # concatenate vectors of marks
  M <- unlist(lapply(list(...), function(x) {x$marks}))
  if(length(M) > 0 && length(M) != length(XY$x)) {
    warning("marks not present in all patterns -- ignored them.")
    M <- NULL
  }
  # determine window
  P <- ..1
  verifyclass(P, "ppp")
  win <- P$window
  #
  # determine type of marks
  if(!is.null(M) && is.factor(P$marks)) {
    M <- factor(M)
    levels(M) <- levels(P$marks)
  }
  #
  ppp(XY$x, XY$y, window=win, marks=M)
}


"is.marked.ppp" <-
function(X, na.action="warn") {
    verifyclass(X, "ppp")
    if(is.null(X$marks))
      return(F)
    if(any(is.na(X$marks)))
      switch(na.action,
             warn = {
               warning(paste("some mark values are NA in the point pattern",
                    deparse(substitute(X))))
             },
             fatal = {
               return(F)
             },
             ignore = {
               return(T)
             }
      )
    return(T)
}

"is.marked" <-
function(X, ...) {
  UseMethod("is.marked")
}

"is.marked.default" <-
  function(...) { return(F) }

"unmark" <-
function(X) {
  verifyclass(X, "ppp")
  X$marks <- NULL
  X
}

"markspace.integral" <-
  function(X) {
  verifyclass(X, "ppp")
  if(!is.marked(X))
    return(1)
  if(is.factor(X$marks))
    return(length(levels(X$marks)))
  else
    stop("Don't know how to compute total mass of mark space")
}
#
#    predictppm.S
#
#	$Revision: 1.11 $	$Date: 2002/01/18 06:57:40 $
#
#    predict.ppm()
#	   From fitted model obtained by mpl(),	
#	   evaluate the fitted trend or conditional intensity 
#	   at a grid/list of other locations 
#
#
# -------------------------------------------------------------------

predict.ppm <-
function(object, newdata, nx = 40, ny = NULL, type="trend", ...) {
#
#	'object' is the output of mpl()
#
  model <- object
  verifyclass(model, "ppm")
#
#       find out what kind of model it is
#
  stationary <- is.stationary.ppm(model)
  poisson <- is.poisson.ppm(model)
  markeddata <- is.marked(model$Q$data)
  notrend <- no.trend.ppm(model)
  trivial <- poisson && notrend

  if(markeddata) {
      mrks <- model$Q$data$marks
  }    
#
#  Create data frame to pass to predict.glm()
#		
  want.grid <- missing(newdata)

  if(want.grid) {
  # create a rectangular grid of nx x ny points in window
    window <- model$Q$data$window
    if(is.null(ny))
      ny <- nx
    xr <- window$xrange
    yr <- window$yrange
  # gam with lo() will not allow extrapolation beyond the range of x,y
  # values actually used for the fit. Check this:
    tums <- termsinformula(model$trend)
    if(any(
           tums == "lo(x)" |
           tums == "lo(y)" |
           tums == "lo(x,y)" |
           tums == "lo(y,x)")
    ) {
      gg <- model$internal$glmdata
      gxr <- range(gg$x[gg$SUBSET])
      gyr <- range(gg$y[gg$SUBSET])
      intersect.ranges <- function(r, s) {
        c(max(r[1],s[1]), min(r[2],s[2]))
      }
      xr <- intersect.ranges(xr, gxr)
      yr <- intersect.ranges(yr, gyr)
    }
  # determine x and y vectors for nx * ny grid
    dx <- diff(xr)/nx
    dy <- diff(yr)/ny
    xvals <- seq(xr[1] + dx/2, xr[2]-dx/2, length=nx)
    yvals <- seq(yr[1] + dy/2, yr[2]-dy/2, length=ny)
    if(!markeddata) {
      # create nx * ny grid
      newdata <- expand.grid(x=xvals, y = yvals)
      # initialise output list
      out <- list(x=xvals, y=yvals)
    } else {
      if(!is.factor(mrks))
        stop("marks are not a factor")
      mrkvals <- levels(mrks)
      # create (nx * ny * nmarks) marked points
      newdata <- expand.grid(x=xvals, y=yvals, marks=mrkvals)
      # initialise output list
      out <- list(x=xvals, y=yvals, marks = mrkvals)
    }
  }

#
######## Set up prediction variables ################################
#
#
# Provide SUBSET variable
#
        if(is.null(newdata$SUBSET))
          newdata$SUBSET <- rep(T, nrow(newdata))
#
# Dig out information used in the original call to mpl().
#        Vnames:     the names for the ``interaction variables''
#        glmdata:    the data frame used for the glm fit
# and set the contrasts according to the convention followed in mpl().
#
  if(!trivial) {
	Vnames <- model$internal$Vnames

        if(exists("is.R") && is.R())
	    glmdata <- model$internal$glmdata
        else        
            assign("glmdata", model$internal$glmdata, f = 1)
  }
  
############  COMPUTE PREDICTION ##############################
#
#   Compute the predicted value z[i] for each row of 'newdata'
#   Store in a vector z and reshape it later
#

###############################################################  
  if(trivial) {
#############  COMPUTE CONSTANT INTENSITY #####################

    lambda <- exp(model$theta[[1]])
    z <- rep(lambda, nrow(newdata))
    
################################################################
  } else if(type == "trend" || poisson) {
#
#############  COMPUTE TREND ###################################
#	
#   set explanatory variables to zero
#	
    zeroes <- rep(0, nrow(newdata))    
    for(vn in Vnames)    
      newdata[[vn]] <- zeroes
#
#   invoke predict.glm()
#  
    z <- predict(model$internal$glmfit, newdata, type="response")

##############################################################  
  } else if(type == "cif" || type =="lambda") {
######### COMPUTE FITTED CONDITIONAL INTENSITY ################
#
# 	
  # set up arguments
    inter <- model$interaction
    X <- model$Q$data
    U <- list(x=newdata$x, y=newdata$y)
    Equal <- outer(X$x, U$x, "==") & outer(X$y, U$y, "==")
    if(markeddata) {
      U$marks <- newdata$marks
      Equal <- Equal & outer(X$marks, U$marks, "==")
    }
  # compute values of potential at the new sample points
    Vnew <- inter$family$eval(X, U, Equal,
                              inter$pot, inter$par, model$correction)
    if(!is.matrix(Vnew))
      stop("internal error: eval.pair.inter() did not return a matrix")
  
  # Insert the potential into the relevant column(s) of `newdata'
    if(ncol(Vnew) == 1)
      # Potential is real valued (Vnew is a column vector)
      # Assign values to a column of the same name in newdata
      newdata[[Vnames]] <- as.vector(Vnew)
      #
    else if(is.null(dimnames(Vnew)[[2]])) {
      # Potential is vector-valued (Vnew is a matrix)
      # with unnamed components.
      # Assign the components, in order of their appearance,
      # to the columns of newdata labelled Vnames[1], Vnames[2],... 
      for(i in seq(Vnames))
        newdata[[Vnames[i] ]] <- Vnew[,i]
      #
    } else {
      # Potential is vector-valued (Vnew is a matrix)
      # with named components.
      # Match variables by name
      for(vn in Vnames)    
        newdata[[vn]] <- Vnew[,vn]
      #
    }
  # invoke predict.glm
  z <- predict(model$internal$glmfit, newdata, type="response")

#################################################################    
  } else
     stop(paste("Unrecognised type \'", type, "\'\n", sep=""))

#################################################################
#
# reshape the result
#
    if(!want.grid) 
      out <- as.vector(z)
    else {
      out$z <-
        if(!markeddata)
          matrix(z, nx, ny)
        else
          array(z, dim=c(nx,ny,length(mrkvals)),
                dimnames=list(NULL, NULL, as.character(mrkvals)))
    }

####################################################################
#
#  
  invisible(out)
}

#
#	quadclass.S
#
#	Class 'quad' to define quadrature schemes
#	in (rectangular) windows in two dimensions.
#
#	$Revision: 4.2 $	$Date: 2002/01/17 09:44:54 $
#
# An object of class 'quad' contains the following entries:
#
#	$data:	an object of class 'ppp'
#		defining the OBSERVATION window, 
#		giving the locations (& marks) of the data points.
#
#	$dummy:	object of class 'ppp'
#		defining the QUADRATURE window, 
#		giving the locations (& marks) of the dummy points.
#	
#	$w: 	vector giving the nonnegative weights for the
#		data and dummy points (data first, followed by dummy)
#
#		w may also have an attribute attr(w, "zeroes")
#               equivalent to (w == 0). If this is absent
#               then all points are known to have positive weights.
#
#       The combined (data+dummy) vectors of x, y coordinates of the points, 
#       and their weights, are extracted using standard functions 
#       x.quad(), y.quad(), w.quad() etc.
#
#-------------------------------------------------------------

quad <- function(data, dummy, w) {
  
  data <- as.ppp(data)
  dummy <- as.ppp(dummy)

  n <- data$n + dummy$n
	
  if(missing(w))
    w <- rep(1, n)
  else if(length(w) != n)
    stop("length of weights vector w is not equal to total number of points")

  if(is.null(attr(w, "zeroes")) && any( w == 0))
	attr(w, "zeroes") <- (w == 0)

  Q <- list(data=data, dummy=dummy, w=w)
  class(Q) <- "quad"

  invisible(Q)
}

# ------------------ extractor functions ----------------------

x.quad <- function(Q) {
  verifyclass(Q, "quad")
  c(Q$data$x, Q$dummy$x)
}

y.quad <- function(Q) {
  verifyclass(Q, "quad")
  c(Q$data$y, Q$dummy$y)
}

w.quad <- function(Q) {
  verifyclass(Q, "quad")
  Q$w
}

n.quad <- function(Q) {
  verifyclass(Q, "quad")
  Q$data$n + Q$dummy$n
}

marks.quad <- function(Q) {
  verifyclass(Q, "quad")
  mdat <- Q$data$marks
  mdum <- Q$dummy$marks
  if(is.null(mdat) && is.null(mdum))
    return(NULL)
  if(is.null(mdat))
    mdat <- rep(NA, Q$data$n)
  if(is.null(mdum))
    mdum <- rep(NA, Q$dummy$n)
  mall <- c(mdat, mdum)
  if(is.factor(mdat) && is.factor(mdum) && all(levels(mdat) == levels(mdum))) {
    mall <- factor(mall)
    levels(mall) <- levels(mdat)
  }
  return(mall)
}

is.data <- function(Q) {
  verifyclass(Q, "quad")
  return(c(rep(T, Q$data$n),
	   rep(F, Q$dummy$n)))
}

equals.quad <- function(Q) {
    # return matrix E such that E[i,j] = (X[i] == U[j])
    # where X = Q$data and U = union.quad(Q)
    n <- Q$data$n
    m <- Q$dummy$n
    E <- matrix(F, nrow=n, ncol=n+m)
    diag(E) <- T
    E
}

  
union.quad <- function(Q) {
  verifyclass(Q, "quad")
  ppp(x= c(Q$data$x, Q$dummy$x),
      y= c(Q$data$y, Q$dummy$y),
      window=Q$dummy$window,
      marks=marks.quad(Q))
}
	
#
#
#      quadscheme.S
#
#      $Revision: 4.1 $    $Date: 2001/08/07 10:14:29 $
#
#      quadscheme()    generate a quadrature scheme from 
#		       data and dummy point patterns.
#
#      quadscheme.spatial()    case where both patterns are unmarked
#
#      quadscheme.replicated() case where data are multitype
#
#
#---------------------------------------------------------------------

quadscheme <- function(data, dummy=default.dummy(data), ...) {
        #
	# generate a quadrature scheme from data and dummy patterns.
	#
	# Other arguments control how the quadrature weights are computed
        #
  mX <- is.marked(data)
  mQ <- is.marked(dummy)
  
  if(!mX && !mQ)
    quadscheme.spatial(data, dummy, ...)
  else if(mX && !mQ)
    quadscheme.replicated(data, dummy, ...)
  else if(!mX && mQ)
    stop("dummy points are marked but data are unmarked")
  else
    stop("marked data and marked dummy points -- sorry, this case is not implemented")
}

quadscheme.spatial <-
  function(data, dummy=default.dummy(data), method="grid", ...) {
        #
	# generate a quadrature scheme from data and dummy patterns.
	#
	# The 'method' may be "grid" or "dirichlet"
	#
	# '...' are passed to gridweights() or dirichlet.weights()
        #
        # quadscheme.spatial:
        #       for unmarked point patterns.
        #
        #       weights are determined only by spatial locations
        #       (i.e. weight computations ignore any marks)
	#
        # No two points should have the same spatial location
        # 

	data <- as.ppp(data)
	dummy <- as.ppp(dummy, data$window)
		# note data$window is the DEFAULT quadrature window
		# unless otherwise specified in 'dummy'

        if(is.marked(data))
          warning("marks in data pattern - ignored")
        if(is.marked(dummy))
          warning("marks in dummy pattern - ignored")
        
	both <- as.ppp(concatxy(data, dummy), dummy$window)
	switch(method,
		grid={
			w <- gridweights(both, window= dummy$window, ...)
		},
		dirichlet = {
			w <- dirichlet.weights(both, window=dummy$window, ...)
		},
		{ 
			stop(paste("unrecognised method \'", method, "\'")) 
		}
	)
	Q <- quad(data, dummy, w)
	invisible(Q)
}

"quadscheme.replicated" <-
  function(data, dummy=default.dummy(data), method="grid", ...) {
        #
	# generate a quadrature scheme from data and dummy patterns.
	#
	# The 'method' may be "grid" or "dirichlet"
	#
	# '...' are passed to gridweights() or dirichlet.weights()
        #
        # quadscheme.replicated:
        #       for multitype point patterns.
        #
        # No two points in 'data'+'dummy' should have the same spatial location

	data <- as.ppp(data)
	dummy <- as.ppp(dummy, data$window)
		# note data$window is the DEFAULT quadrature window
		# unless otherwise specified in 'dummy'

        if(!is.marked(data))
          stop("data pattern does not have marks")
        if(is.marked(dummy))
          warning("dummy points have marks --- ignored")

        # first, ignore marks and compute spatial weights
        P <- quadscheme.spatial(unmark(data), dummy, method, ...)
        W <- w.quad(P)
        iz <- is.data(P)
        Wdat <- W[iz]
        Wdum <- W[!iz]

        # find the set of all possible marks

        if(!is.factor(data$marks))
          stop("data$marks is not a factor")
        markset <- levels(data$marks)
        nmarks <- length(markset)
        
        # replicate dummy points, one copy for each possible mark
        # -> dummy x {1,..,K}
        
        dumdum <- replicate(dummy, markset)
        Wdumdum <- rep(Wdum, nmarks)
        
        # also make dummy marked points at same locations as data points
        # but with different marks

        dumdat <- replicate(unmark(data), markset)
        Wdumdat <- rep(Wdat, nmarks)
        Mdumdat <- dumdat$marks
        
        Mrepdat <- rep(data$marks, nmarks)

        ok <- (Mdumdat != Mrepdat)
        dumdat <- dumdat[ok,]
        Wdumdat <- Wdumdat[ok]

        # combine the two dummy patterns

        dumb <- superimpose(dumdum, dumdat)
        Wdumb <- c(Wdumdum, Wdumdat)

        # wrap up

	Q <- quad(data, dumb, c(Wdat, Wdumb))
	invisible(Q)
}


"replicate" <-
function(pp, markset, fac=T) {
  # given an unmarked point pattern 'pp'
  # and a finite set of marks,
  # create the marked point pattern which is
  # the Cartesian product, consisting of all pairs (u,k)
  # where u is a point of 'pp' and k is a mark in 'markset'
  nmarks <- length(markset)
  result <- ppp(
                rep(pp$x, nmarks),
                rep(pp$y, nmarks),
                window=pp$window,
                marks=rep(markset, rep(pp$n, nmarks))
  )
  if(fac)
    result$marks <- factor(result$marks, levels=markset)
  result
}
#
#    random.S
#
#    Functions for generating random point patterns
#
#    $Revision: 4.4 $   $Date: 2001/11/26 09:51:16 $
#
#
#    runifpoint()      n i.i.d. uniform random points ("binomial process")
#
#    runifpoispp()     uniform Poisson point process
#
#    rpoispp()         general Poisson point process (rejection method)
#
#    rMaternI()        Mat'ern model I 
#    rMaternII()       Mat'ern model II
#    rSSI()            Simple Sequential Inhibition process
#
#    rNeymanScott()    Neyman-Scott process (generic)
#    rMatClust()       Mat'ern cluster process
#    rThomas()         Thomas process
#
#
#
#    Examples:
#          u01 <- owin(0:1,0:1)
#          plot(runifpoispp(100, u01))
#          X <- rpoispp(function(x,y) {100 * (1-x/2)}, 100, u01)
#          X <- rpoispp(function(x,y) {ifelse(x < 0.5, 100, 20)}, 100)
#          plot(X)
#          plot(rMaternI(100, 0.02))
#          plot(rMaternII(100, 0.05))
#

"runifrect" <-
  function(n, win=owin(c(0,1),c(0,1)))
{
  # no checking
      x <- runif(n, min=win$xrange[1], max=win$xrange[2])
      y <- runif(n, min=win$yrange[1], max=win$yrange[2])  
      return(ppp(x, y, window=win))
}

"runifdisc" <-
  function(n, r=1, x=0, y=0)
{
  # i.i.d. uniform points in the disc of radius r and centre (x,y)
  theta <- runif(n, min=0, max= 2 * pi)
  s <- sqrt(runif(n, min=0, max=r^2))
  return(list(x = x + s * cos(theta), y = y + s * sin(theta)))
}


"runifpoint" <-
  function(n, win=owin(c(0,1),c(0,1)), giveup=1000)
{
    win <- as.owin(win)

    if(win$type == "rectangle")
      return(runifrect(n, win))

    # otherwise - window is irregular
    
    # rejection method
    # initialise empty pattern
    x <- numeric(0)
    y <- numeric(0)
    X <- ppp(x, y, window=win)
    #
    # rectangle in which trial points will be generated
    box <- bounding.box(win)
    # 
    ntries <- 0
    repeat {
      ntries <- ntries + 1
      # generate trial points in batches of n
      qq <- runifrect(n, box) 
      # retain those which are inside 'win'
      qq <- qq[, win]
      # add them to result
      browser()
      X <- superimpose(X, qq)
      # if we have enough points, exit
      if(X$n > n) 
        return(X[1:n])
      else if(X$n == n)
        return(X)
      # otherwise get bored eventually
      else if(ntries >= giveup)
        stop(paste("Gave up after", giveup * n, "trials,",
                   np, "points accepted"))
    }
}

"runifpoispp" <-
function(lambda, win = owin(c(0,1),c(0,1))) {
    win <- as.owin(win)
    if(!is.numeric(lambda) || length(lambda) > 1 || lambda <= 0)
      stop("Intensity lambda must be a single number > 0")

    # generate Poisson process in enclosing rectangle 
    box <- bounding.box(win)
    mean <- lambda * area.owin(box)
    n <- rpois(1, mean)
    X <- runifpoint(n, box)

    # trim to window
    if(win$type != "rectangle")
      X <- X[, win]  

    return(X)
}

"rpoispp" <-
  function(lambda, max, win = owin(c(0,1),c(0,1))) {
    # arguments:
    #     win     observation window (of class 'owin')
    #     lambda  intensity - constant or function(x,y)
    #     max     maximum possible value of lambda(x,y)
    win <- as.owin(win)
    
    if(is.numeric(lambda))
      # uniform Poisson
      return(runifpoispp(lambda, win))
    # inhomogeneous Poisson - use rejection filtering 
    X <- runifpoispp(max, win)  # includes sanity checks on `max'
    prob <- lambda(X$x, X$y)/max
    u <- runif(X$n)
    retain <- (u <= prob)
    X <- X[retain, ]
    return(X)
}
    
"rMaternI" <-
  function(lambda, r, win = owin(c(0,1),c(0,1)))
{
    win <- as.owin(win)
    X <- runifpoispp(lambda, win)
    d <- nndist(X$x, X$y)
    qq <- X[d > r]
    return(qq)
}
    
"rMaternII" <-
  function(lambda, r, win = owin(c(0,1),c(0,1)))
{
    win <- as.owin(win)

    X <- runifpoispp(lambda, win)

    # matrix of pairwise distances
    d <- pairdist(X$x, X$y)
    close <- (d <= r)

    # random order 1:n
    age <- sample(seq(X$n), X$n, replace=F)
    earlier <- outer(age, age, ">")

    conflict <- close & earlier
    # delete <- apply(conflict, 1, any)
    delete <- matrowany(conflict)
    
    qq <- X[ !delete]
    return(qq)
}
  
"rSSI" <-
  function(r, n, win = owin(c(0,1),c(0,1)), giveup = 1000)
{
     # Simple Sequential Inhibition process
     # fixed number of points
     # Naive implementation, proposals are uniform
     win <- as.owin(win)
     X <- ppp(numeric(0),numeric(0), window=win)
     r2 <- r^2
     if(n * pi * r2/4  > area.owin(win))
       stop(paste("Window is too small to fit", n, "points",
                  "at minimum separation", r))
     ntries <- 0
     while(ntries < giveup) {
       ntries <- ntries + 1
       qq <- runifpoint(1, win)
       x <- qq$x[1]
       y <- qq$y[1]
       if(X$n == 0 || all(((x - X$x)^2 + (y - X$y)^2) > r2))
         X <- superimpose(X, qq)
       if(X$n == n)
         return(X)
     }
     warning(paste("Gave up after", giveup,
                "attempts with only", X$n, "points placed"))
     return(X)
}

"rNeymanScott" <-
  function(lambda, rmax, rcluster, win = owin(c(0,1),c(0,1)), ...)
{
  # Generic Neyman-Scott process
  # Implementation for bounded cluster radius
  #
  # 'rcluster' is a function(x,y) that takes the coordinates
  # (x,y) of the parent point and generates a list(x,y) of offspring
  #
  # "..." are arguments to be passed to 'rcluster()'
  #
  # Generate parents in dilated window
  frame <- bounding.box(win)
  dilated <- owin(frame$xrange + c(-rmax, rmax),
                  frame$yrange + c(-rmax, rmax))
  parents <- runifpoispp(lambda, win=dilated)
  #
  result <- ppp(numeric(0), numeric(0), window = win)
  
  if(parents$n == 0)
    return(result)
  
  for(i in seq(parents$n)) {
    # generate random offspring of i-th parent point
    cluster <- rcluster(parents$x[i], parents$y[i], ...)
    cluster <- ppp(cluster$x, cluster$y, window=frame)
    # trim to window
    cluster <- cluster[,win]
    # add to pattern
    result <- superimpose(result, cluster)
  }
  return(result)
}  

"rMatClust" <-
  function(lambda, r, mu, win = owin(c(0,1),c(0,1)))
{
  # Matern Cluster Process with Poisson (mu) offspring distribution
  #
  poisclus <-  function(x0, y0, radius, mu) {
                           n <- rpois(1, mu)
                           return(runifdisc(n, radius, x0, y0))
                         }
  result <- rNeymanScott(lambda, r, poisclus, win, radius=r, mu=mu)
  return(result)
}
    
"rThomas" <-
  function(lambda, sigma, mu, win = owin(c(0,1),c(0,1)))
{
  # Thomas process with Poisson(mu) number of offspring
  # at isotropic Normal(0,sigma^2) displacements from parent
  #
  thomclus <-  function(x0, y0, sigma, mu) {
                           n <- rpois(1, mu)
                           x <- rnorm(n, mean=x0, sd=sigma)
                           y <- rnorm(n, mean=y0, sd=sigma)
                           return(list(x=x, y=y))
                         }
  result <- rNeymanScott(lambda, 4 * sigma, thomclus, win, sigma=sigma, mu=mu)
  return(result)
}
  
#
#
#    rmh.R
#
#    $Revision: 1.4 $     $Date: 2001/12/12 14:35:06 $
#
#
#

rmh <- function(cif,par,w,ntypes=0,ptypes=NULL,tpar=NULL,n.start,
                expand=NULL,periodic=F,nrep=1e6,p=0.9,q=0.5,
                iseed=NULL,nverb=0) {
#
# Function rmh.  To simulate realizations of 2-dimensional point
# patterns, given the conditional intensity function of the 
# underlying process, via the Metropolis-Hastings algorithm
#

# Check that cif is available:
if(!is.loaded(symbol.For(cif)))
	stop(paste("Unrecognized cif: ",cif,".\n",sep=""))

# Turn the name of the cif into a number
name.list <- c('strauss','straush','sftcr','straussm','straushm',
               'dig1','dig2','geyer')
nmbr <- match(cif,name.list)
if(is.na(nmbr)) stop("Name of cif not in name list.\n")

# Check for compatibility of ntypes and length of ptypes.
if(ntypes <= 1) ptypes <- 1
else {
	if(is.null(ptypes)) ptypes <- rep(1/ntypes,ntypes)
	if(length(ptypes) != ntypes | sum(ptypes) != 1)
		stop("Arguments ntypes and/or ptypes do not make sense.\n")
}

# Turn w into an object of class "owin" if necessary; then
# turn it (back) into a vector of length 4 determining the
# enclosing box (which may simply be the window if the
# original window was rectangular).
w.save <- as.owin(w)
rw <- c(w.save$xrange,w.save$yrange)

# Set the ``period''; if periodic make sure expand is 1.
if(periodic) {
	if(w.save$type != "rectangle") {
		whinge <- paste("\"periodic\" only makes sense",
                                "for rectanglar windows.\n")
		stop(whinge)
	}
	if(is.null(expand)) expand <- 1
	else if(expand > 1) stop("If periodic, expand must be 1.\n")
	period <- c(rw[2] - rw[1], rw[4] - rw[3])
} else {
	if(is.null(expand)) expand <- 2
	period <- c(-1,-1)
}

# Adjust n.start; it is/should be given as the ``expected'' number
# of points in the actual window.  It should be magnified to the
# expected number of points in the bounding box ``rw''.
a1 <- area.owin(w.save)
a2 <- area.owin(as.owin(rw))
n.start <- ceiling(a2*n.start/a1)

# Set need.aux to F (it gets set to T iff cif == 'geyer').
need.aux <- F

# Do some rudimentary pre-processing of the par and tpar arguments.
# Save the supplied values of par and tpar first.
par.save  <- par
tpar.save <- tpar

# 1. Strauss.
if(cif=="strauss") {
	if(length(par) != 3) {
		cat("For strauss cif, par should be a vector of\n")
		cat("length 3, consisting of beta, gamma and r.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative parameters.\n")
	if(par[2] > 1)
		stop("For Strauss processes, gamma must be <= 1.\n")
}

# 2. Strauss with hardcore.
if(cif=="straush") {
	if(length(par) != 4) {
		cat("For straush cif, par should be a vector of\n")
		cat("length 4, consisting of beta, gamma, r, and r_0.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative parameters.\n")
}

# 3. Softcore.
if(cif=="sftcr") {
	if(length(par) != 3) {
		cat("For sftcr cif, par should be a vector of\n")
		cat("length 3, consisting of beta, sigma and kappa.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative  parameters.\n")
	if(par[3] > 1)
		stop("For Softcore processes, kappa must be <= 1.\n")
}

# 4. Marked Strauss.
if(cif=="straussm") {
	if(ntypes<=1)
		stop("Argument ntypes must be at least 2 for straussm.\n")
	if(!is.atomic(par)) {
		if(!is.list(par)) {
			cat("For straussm cif, par must be either a\n")
			cat("vector, or a list with components beta,\n")
			cat("gamma, and r.\n")
			stop("Bailing out.\n")
		}
		beta <- par$beta
		if(length(beta) != ntypes)
			stop("Component beta of par is of wrong length.\n")
		gamma <- par$gamma
		if(!is.matrix(gamma) | length(gamma) != ntypes^2)
			stop("Component gamma of par is of wrong shape.\n")
		r <- par$r
		if(!is.matrix(r) | length(r) != ntypes^2)
			stop("Component r of par is the wrong shape.\n")
		gamma <- t(gamma)[row(gamma)>=col(gamma)]
		r <- t(r)[row(r)>=col(r)]
		par <- c(beta,gamma,r)
	}
	if(any(par[!is.na(par)]<0))
		stop("Negative  parameters.\n")
}

# 5. Marked Strauss with hardcore.
if(cif=="straushm") {
	if(ntypes<=1)
		stop("Argument ntypes must be at least 2 for straushm.\n")
	if(!is.atomic(par)) {
		if(!is.list(par)) {
			cat("For straushm cif, par must be either a\n")
			cat("vector, or a list with components beta,\n")
			cat("gamma, r, and rhc.\n")
			stop("Bailing out.\n")
		}
		beta <- par$beta
		if(length(beta) != ntypes)
			stop("Component beta of par is of wrong length.\n")

		gamma <- par$gamma
		if(!is.matrix(gamma) | length(gamma) != ntypes^2)
			stop("Component gamma of par is of wrong shape.\n")

		r <- par$r
		if(!is.matrix(r) | length(r) != ntypes^2)
			stop("Component r of par is the wrong shape.\n")

		rhc <- par$rhc
		if(!is.matrix(rhc) | length(rhc) != ntypes^2)
			stop("Component rhc of par is the wrong shape.\n")

		gamma <- t(gamma)[row(gamma)>=col(gamma)]
		r <- t(r)[row(r)>=col(r)]
		rhc <- t(rhc)[row(rhc)>=col(rhc)]

		par <- c(beta,gamma,r,rhc)
	}
	if(any(par[!is.na(par)]<0))
		stop("Negative  parameters.\n")
}

# 6. Using interaction function number 1 from Diggle, Gates,
#    and Stibbard.

if(cif=="dig1") {
	if(length(par) != 2) {
		cat("For dig1 cif, par should be a vector of\n")
		cat("length 2, consisting of beta and rho.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative parameters.\n")
}

# 7. Using interaction function number 2 from Diggle, Gates,
#    and Stibbard.

if(cif=="dig2") {
	if(length(par) != 4) {
		cat("For dig2 cif, par should be a vector of\n")
		cat("length 4, consisting of beta, kappa, delta\n")
		cat("and rho.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative parameters.\n")
	if(par[3] >= par[4])
		stop("Radius delta must be less than radius rho.\n")
}

# 8. The Geyer conditional intensity function.

if(cif=="geyer") {
	if(length(par) != 4) {
		cat("For the geyer cif, par should be a vector of\n")
		cat("length 4, consisting of beta, gamma, r\n")
		cat("and s.\n")
		stop("Bailing out.\n")
	}
	if(any(par<0))
		stop("Negative parameters.\n")
	if(par[4] > .Machine$integer.max-100)
		par[4] <- .Machine$integer.max-100
	need.aux <- T
}

# Calculate the degree(s) of the polynomial(s) in the log polynomial
# trend(s), if any, and tack it/them onto tpar.  Note that
# if tpar is NULL, this will make tpar into a scalar equal to 0.
if(ntypes > 1 & (!is.null(tpar))) {
	if(!is.atomic(tpar)) {
		if(!is.list(tpar) || length(tpar) != ntypes)
			stop("Argument tpar has wrong form.\n")
		tmp <- NULL
		for(cc in tpar) {
			nc <- length(cc)
			nd <- (sqrt(9+8*nc) - 3)/2
			if(nd%%1 > .Machine$double.eps) {
				stpms <- "Length of a component"
				stpms <- paste(stpms,"does not make sense.\n")
				stop(stpms)
			}
			tmp <- c(tmp,nd)
		}
		tpar <- c(ntypes,tmp,unlist(tpar))
	}
}
else {
	nt <- length(tpar)
	nd <- (sqrt(9+8*nt) - 3)/2
	if(nd%%1 > .Machine$double.eps)
	stop("Length of tpar does not make sense.\n")
	tpar <- c(nd,tpar)
}

# Set the 3-vector of integer seeds needed by the subroutine arand:
if(is.null(iseed)) iseed <- sample(1:1000000,3)
iseed.save <- iseed

# Prepare vectors x and y (and perhaps marks) to hold the generated
# process; note that we are guessing at how big they will need to be.
# We start off with twice the length of the ``initial state'',
# and structure things so that the storage space may be incremented
# without losing the ``state'' which has already been generated.

npts  <- if(expand > 1) ceiling(expand*n.start) else n.start
n2    <- 2*npts
x     <- numeric(n2)
y     <- numeric(n2)
marks <- if(ntypes > 1) numeric(n2) else 0
aux   <- if(need.aux) numeric(n2) else 0
npmax <- 0
mrep  <- 0 # Indicates starting out; subroutine methas will
           # reset mrep to 1 after generating ``initial state''.

# Build the expanded window within which to suspend the actual
# window of interest (in order to approximate the simulation of a
# windowed process, rather than a process existing only in the given
# window.  If expand == 1, then we are simulating the latter.  The
# larger ``expand'' is, the better we approximate the former.  Note
# that any value of ``expand'' smaller than 1 is treated as if it
# were 1.
if(expand>1) {
	xdim <- rw[2]-rw[1]
	ydim <- rw[4]-rw[3]
	fff  <- (sqrt(expand)-1)*0.5
	erw  <- c(rw[1] - fff*xdim,rw[2] + fff*xdim,
                  rw[3] - fff*ydim,rw[4] + fff*ydim)
}
else erw <- rw

# The repetion is to allow the storage space to be incremented if
# necessary.
repeat {
	npmax <- npmax + n2
# Call the Metropolis-Hastings simulator:
	rslt <- .Fortran(
			"methas",
			nmbr=as.integer(nmbr),
			rw=as.double(erw),
			par=as.double(par),
			period=as.double(period),
			tpar=as.double(tpar),
			ntypes=as.integer(ntypes),
			ptypes=as.double(ptypes),
			iseed=as.integer(iseed),
			nrep=as.integer(nrep),
			mrep=as.integer(mrep),
			p=as.double(p),
			q=as.double(q),
			npmax=as.integer(npmax),
			nverb=as.integer(nverb),
			x=as.double(x),
			y=as.double(y),
			marks=as.integer(marks),
			aux=as.integer(aux),
			npts=as.integer(npts)
		)

# If npts > npmax we've run out of storage space.  Tack some space
# onto the end of the ``state'' already generated, increase npmax
# correspondingly, and re-call the hasmet subroutine.  Note that
# mrep is the number of the repetion on which things stopped due
# to lack of storage; so we start again at the ***beginning*** of
# the mrep repetion.
	npts <- rslt$npts
	if(npts <= npmax) break
	cat('Number of points greater than ',npmax,';\n',sep='')
	cat('increasing storage space and continuing.\n')
	x     <- c(rslt$x,numeric(n2))
	y     <- c(rslt$y,numeric(n2))
	marks <- if(ntypes>1) c(rslt$marks,numeric(n2)) else 0
	aux   <- if(need.aux) c(rslt$aux,numeric(n2)) else 0
	mrep  <- rslt$mrep
	iseed <- rslt$iseed
	npts  <- npts-1
}

x <- rslt$x[1:npts]
y <- rslt$y[1:npts]
if(ntypes>1) marks <- rslt$marks[1:npts]
tmp <- as.ppp(list(x=x,y=y),W=as.owin(erw))
if(ntypes>1) tmp$marks <- factor(marks)

# Now window the returned process by the original window:
tmp <- tmp[,w.save]

# Append to the result information about how it was generated.
tmp$info <- list(cif=cif,par=par.save,tpar=tpar.save,n.start=n.start,
                  nrep=nrep,p=p,q=q,expand=expand,periodic=periodic,
                  iseed=iseed.save)
class(tmp) <- "ppp"
tmp
}
#
#	rotate.S
#
#	$Revision: 1.1 $	$Date: 2001/11/21 09:07:22 $
#

rotxy <- function(X, angle=pi/2) {
  co <- cos(angle)
  si <- sin(angle)
  list(x = co * X$x - si * X$y,
       y = si * X$x + co * X$y)
}

"rotate.owin" <- function(W, angle=pi/2) {
  verifyclass(W, "owin")
  switch(W$type,
         rectangle={
           # convert rectangle to polygon
           P <- owin(W$xrange, W$yrange, poly=
                     list(x=W$xrange[c(1,2,2,1)],
                          y=W$yrange[c(1,1,2,2)]))
           # call polygonal case
           return(rotate.owin(P, angle))
         },
         polygonal={
           # First rotate the polygonal boundaries
           bdry <- lapply(W$bdry, rotxy, angle=angle)
           # Compute bounding box of new polygons
           xr <- range(unlist(lapply(bdry, function(a) a$x)))
           yr <- range(unlist(lapply(bdry, function(a) a$y)))
           # wrap up
           return(owin(xr, yr, poly=bdry))
         },
         mask={
           stop("Sorry, \'rotate.owin\' is not yet implemented for masks")
         },
         stop("Unrecognised window type")
         )
}

"rotate.ppp" <- function(X, angle=pi/2) {
  verifyclass(X, "ppp")
  r <- rotxy(X, angle)
  w <- rotate.owin(X$window, angle)
  return(ppp(r$x, r$y, window=w, marks=X$marks))
}


"rotate" <- function(X, ...) {
  UseMethod("rotate")
}

  
#
#
#    saturated.S
#
#    $Revision: 1.2 $	$Date: 2001/08/07 11:52:17 $
#
#    Saturated pairwise process with user-supplied potential
#
#    Saturated()  create a saturated pairwise process
#                 [an object of class 'interact']
#                 with user-supplied potential
#	
#
# -------------------------------------------------------------------
#	

Saturated <- function(pot, name) {
  if(missing(name))
    name <- "Saturated process with user-defined potential"
  
  out <- 
  list(
         name     = name,
         family    = pairsat.family,
         pot      = pot,
         par      = NULL,
         parnames = NULL,
         init     = NULL,
         update   = NULL, 
         print = function(self) {
           cat(paste(self$name, "\n"))
           cat("Potential function:\n")
           print(self$pot)
           invisible()
         }
  )
  class(out) <- "interact"
  return(out)
}
#
#
#    softcore.S
#
#    $Revision: 1.2 $   $Date: 2000/07/11 10:52:01 $
#
#    Soft core processes.
#
#    Softcore()    create an instance of a soft core process
#                 [an object of class 'interact']
#
#
# -------------------------------------------------------------------
#

Softcore <- function(kappa) {
  out <- 
  list(
         name     = "Soft core process",
         family   = pairwise.family,
         pot      = function(d, par) {
                        -d^(-2/par$kappa)
                    },
         par      = list(kappa = kappa),
         parnames = "Exponent kappa",
         init     = function(self) {
                      kappa <- self$par$kappa
                      if(!is.numeric(kappa) || length(kappa) != 1 ||
                         kappa <= 0 || kappa >= 1)
                       stop("Exponent kappa must be a positive \
number less than 1")
                    },
         update = NULL,  # default OK
         print = NULL    # default OK
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}

#
#
#    strauss.S
#
#    $Revision: 1.5 $	$Date: 2001/07/29 07:25:54 $
#
#    The Strauss process
#
#    Strauss()    create an instance of the Strauss process
#                 [an object of class 'interact']
#	
#
# -------------------------------------------------------------------
#	

Strauss <- function(r) {
  out <- 
  list(
         name     = "Strauss process",
         family    = pairwise.family,
         pot      = function(d, par) {
                         (d <= par$r)
                    },
         par      = list(r = r),
         parnames = "interaction distance",
         init     = function(self) {
                      r <- self$par$r
                      if(!is.numeric(r) || length(r) != 1 || r <= 0)
                       stop("interaction distance r must be a positive number")
                    },
         update = NULL,  # default OK
         print = NULL    # default OK
  )
  class(out) <- "interact"
  out$init(out)
  return(out)
}
#
#
#    strausshard.S
#
#    $Revision: 1.5 $	$Date: 2001/08/10 09:13:59 $
#
#    The Strauss/hard core process
#
#    StraussHard()     create an instance of the Strauss-hardcore process
#                      [an object of class 'interact']
#	
#
# -------------------------------------------------------------------
#	

StraussHard <- function(r, hc) {
  out <- 
  list(
         name   = "Strauss - hard core process",
         family  = pairwise.family,
         pot    = function(d, par) {
           v <- ifelse(d <= par$r, 1, 0)
           v[ d <= par$hc ] <-  (-Inf)
           v
         },
         par    = list(r = r, hc = hc),
         parnames = c("interaction distance",
                      "hard core distance"), 
         init   = function(self) {
           r <- self$par$r
           hc <- self$par$hc
           if(!is.numeric(hc) || length(hc) != 1 || hc <= 0)
             stop("hard core distance hc must be a positive number")
           if(!is.numeric(r) || length(r) != 1 || r <= hc)
             stop("interaction distance r must be a number greater than hc")
         },
         update = NULL,       # default OK
         print = NULL         # default OK
  )
  class(out) <- "interact"
  (out$init)(out)
  return(out)
}
#
#	tryFGJKest.S
#
#	Test the F, G, J and K estimation routines
#
#	$Revision: 4.2 $ $Date: 2001/11/13 04:31:36 $
#
#   Note: this is not the most efficient way to calculate F, G, J and K
#         if all four are required. See the source for 'allstats'
#
################################################################
#
"try.FGJKest"<-
function(niter = 20, lambda = 25, r = seq(0, sqrt(2), 0.02), eps=0.01, slow=F)
{
	Frs <- 
	Fkm <- 
	Grs <- 
	Gkm <- 
	Jrs <- 
	Jkm <- 
	Kbord <- matrix(0, nrow = niter, ncol = length(r))

	cat("computing realisation ")
	for(i in 1:niter) {
		cat(paste(i))

		pp <- rpoispp(lambda)
                cat("(")
                cat("F")
                FF <- Fest(pp, eps, r)
		Fkm[i,  ] <- FF$km
		Frs[i,  ] <- FF$rs
                cat("G")
                G <- Gest(pp, r)
		Gkm[i,  ] <- G$km
		Grs[i,  ] <- G$rs
                cat("J")
                J <- Jest(pp, eps, r)
		Jkm[i,  ] <- J$km
		Jrs[i,  ] <- J$rs
                cat("K") ;
                K <- Kest(pp, r, slow=slow)
		Kbord[i,  ] <- K$border
                cat("), ")
	}
	cat("Done.\n")

        oldpar <- par(ask=T)
	plotteststuff <- function(r, mat, correction, symb, trueval) {
          	rrange <- range(r)
                mrange <- range(c(mat, trueval), na.rm=T)
                plot(rrange, mrange, xlab="r", ylab=symb, sub=correction,
                     main=paste("estimated", symb), type="n")
                niter <- nrow(mat)
                for(i in 1:niter) {
                  lines(r, mat[i,  ])
                }
                lines(r, trueval, lty=2)
                dev <- mat - t(matrix(trueval, ncol=nrow(mat), nrow=ncol(mat)))
                derange <- range(c(dev,0), na.rm=T)
                plot(rrange, derange, xlab="r", ylab="Deviation",
                     sub=correction,
                     main=paste("deviation of estimated", symb), type="n")
                for(i in 1:niter) {
                  lines(r, dev[i,  ])
                }
                abline(0,0,lty=2)
                bias <- apply(mat, 2, mean, na.rm=T) - trueval
                brange <- range(c(bias, 0), na.rm=T)
                plot(r, bias, xlab="r", ylab="Bias",
                     sub=correction, ylim=brange,
                     main=paste("Bias of estimated", symb), type="l")
                abline(0,0,lty=2)
                
                varnaok <- function(x) {
                  nbg <- is.na(x)
                  if(all(nbg))
                    NA
                  else
                    var(x[!nbg])   # should work in all dialects
                }
                sd <- sqrt(apply(mat, 2, varnaok))
# Splus 5.1:    sd <- sqrt(apply(mat, 2, var, na.method="available"))
# R:            sd <- sqrt(apply(mat, 2, var, na.rm=T))
                
                plot(r, sd, xlab="r", ylab="SD",
                     sub=correction,
                     main=paste("Standard deviation of estimated", symb),
                     type="l")
                invisible(NULL)
              }

        trueK  <- pi * r^2
        trueFG <- 1 - exp( - lambda * pi * r^2)
        trueJ <- rep(1, length(r))

        plotteststuff(r, Fkm, "Kaplan-Meier", "F", trueFG)
        plotteststuff(r, Frs, "reduced sample", "F", trueFG)
        plotteststuff(r, Gkm, "Kaplan-Meier", "G", trueFG)
        plotteststuff(r, Grs, "reduced sample", "G", trueFG)
        plotteststuff(r, Jkm, "Kaplan-Meier", "J", trueJ)
        plotteststuff(r, Jrs, "reduced sample", "J", trueJ)
        plotteststuff(r, Kbord, "border method", "K", trueK)

        par(oldpar)
	invisible(NULL)
}
#
#	tryKcross.S
#
#	Test the routine Kcross()
#
#	$Revision: 4.2 $ $Date: 2001/11/13 04:31:36 $
#
################################################################
#
"try.Kcross"<-
function(niter = 20, lambda1 = 25, lambda2 = 25, r = seq(0, 1, 0.02), R=0.2)
{
	lambda <- lambda1 + lambda2
	probs <- c(lambda1,lambda2)/lambda
	
	k <- matrix(0, nrow = niter, ncol = length(r))
	k2 <- k
	cat("computing realisation ")
	for(i in 1:niter) {
		cat(paste(i,", ", sep=""))

                X <- rpoispp(lambda)
		X$marks <- factor(sample(1:2, X$n, prob=probs, replace=T))
		out <- Kcross(X, "1", "2", r)
		k[i,  ] <- out$border
		k2[i,  ] <- out$bord.modif
	}
	cat("\n")

        # restrict to r in [0,R]
        ok <- (r <= R)
        r <- r[ok]
        k <- k[, ok]
        k2 <- k2[, ok]
        
	truek <- pi * r^2
	rrange <- range(r)
	krange <- range(c(k, k2, truek), na.rm=T)
	
	plot(rrange, krange, xlab = "r", ylab = "Kcross", 
				main = "border method", type = "n")
	for(i in 1:niter) {
		lines(r, k[i,  ])
	}
	plot(krange, krange, xlab = "true Kcross", ylab = "estimated Kcross", 
				main = "border method", type = "n")
	for(i in 1:niter) {
		lines(truek, k[i,  ])
	}
	plot(rrange, krange, xlab = "r", ylab = "Kcross", 
				main = "border (modified)", type = "n")
	for(i in 1:niter) {
		lines(r, k2[i,  ])
	}
	plot(krange, krange, xlab = "true Kcross", ylab = "estimated Kcross", 
				main = "border (modified)", type = "n")
	for(i in 1:niter) {
		lines(truek, k2[i,  ])
	}
	invisible(NULL)
}
#
#    util.S    miscellaneous utilities
#
#    $Revision: 1.3 $    $Date: 2001/07/31 02:45:19 $
#
#  (a) for matrices only:
#
#    matrowany(X) is equivalent to apply(X, 1, any)
#    matrowall(X) "   "  " "  "  " apply(X, 1, all)
#    matcolany(X) "   "  " "  "  " apply(X, 2, any)
#    matcolall(X) "   "  " "  "  " apply(X, 2, all)
#
#  (b) for 3D arrays only:
#    apply23sum(X)  "  "   "  " apply(X, c(2,3), sum)
#

matrowsum <- function(x) {
  x %*% rep(1, ncol(x))
}

matcolsum <- function(x) {
  rep(1, nrow(x)) %*% x
}
  
matrowany <- function(x) {
  (matrowsum(x) > 0)
}

matrowall <- function(x) {
  (matrowsum(x) == ncol(x))
}

matcolany <- function(x) {
  (matcolsum(x) > 0)
}

matcolall <- function(x) {
  (matcolsum(x) == nrow(x))
}

########
    # hm, this is SLOWER

apply23sum <- function(x) {
  dimx <- dim(x)
  if(length(dimx) != 3)
    stop("x is not a 3D array")
  result <- array(0, dimx[-1])

  nz <- dimx[3]
  for(k in 1:nz) {
    result[,k] <- matcolsum(x[,,k])
  }
  result
}
    
#
#	weights.S
#
#	Utilities for computing quadrature weights
#
#	$Revision: 4.2 $	$Date: 2001/12/11 16:06:27 $
#
#
# Main functions:
		
#	gridweights()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. Given an arbitrary
#			    pattern of (data + dummy) points derive the
#			    'counting weights'.
#
#	dirichlet.weights() Compute the areas of the tiles of the
#			    Dirichlet tessellation generated by the 
#			    given pattern of (data+dummy) points,
#			    restricted to the window.
#	
# Auxiliary functions:	
#			
#       countingweights()   compute the counting weights
#                           for a GENERIC tiling scheme and an arbitrary
#			    pattern of (data + dummy) points,
#			    given the tile areas and the information
#			    that point number k belongs to tile number id[k]. 
#
#
#	gridindex()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. 
#			    Compute tile membership for arbitrary x,y.
#				    
#       discretise()        1-dimensional analogue of gridindex()
#
#
#-------------------------------------------------------------------
	
countingweights <- function(id, areas, check=T) {
	#
	# id:        cell indices of n points
	#                     (length n, values in 1:k)
	#
	# areas:     areas of k cells 
	#                     (length k)
	#
    id <- factor(id, levels=seq(areas))
    counts <- table(id)
    w <- areas[id] / counts[id]     # ensures denominator > 0
#	
# that's it; but check for funny business
#
    zerocount <- (counts == 0)
    zeroarea <- (areas == 0)
    if(any(!zeroarea & zerocount))
	warning("some tiles with positive area do not contain any points")
    if(any(!zerocount & zeroarea)) {
	warning("Some tiles with zero area contain points")
	warning("Some weights are zero")
	attr(w, "zeroes") <- zeroarea[id]
    }
#
    names(w) <- NULL
    w
}

gridindex <- function(x, y, xrange, yrange, nx, ny) {
	#
	# The box with dimensions xrange, yrange is divided
	# into nx * ny cells.
	#
	# For each point (x[i], y[i]) compute the index (ix, iy)
	# of the cell containing the point.
	# 
	ix <- discretise(x, xrange, nx)
	iy <- discretise(y, yrange, ny)
	#
	return(list(ix=ix, iy=iy, index=(iy-1) * nx + ix))
}

discretise <- function(x, xrange, nx) {
	i <- ceiling( nx * (x - xrange[1])/diff(xrange))
	i <- pmax(1, i)
	i <- pmin(i, nx)
	i
}

gridweights <- function(X, nx, ny, window=NULL) {
	#
	# Compute counting weights based on a regular tessellation of the
	# window frame into nx * ny rectangular tiles.
	#
	# Arguments X and (optionally) 'window' are interpreted as a
	# point pattern.
	#
	# The window frame is divided into a regular nx * ny grid
	# of rectangular tiles. The counting weights based on this tessellation
	# are computed for the points (x, y) of the pattern.
	#
	
	X <- as.ppp(X, window)
	x <- X$x
	y <- X$y
	win <- X$window

        if(missing(nx))
          nx <- default.ngrid(X)
        if(missing(ny))
          ny <- nx

	# classify each point according	to its tile
	
	id <- gridindex(x, y, win$xrange, win$yrange, nx, ny)$index

	# compute tile areas
	if(win$type == "rectangle") {

		tilearea <- area.owin(win)/(nx * ny)
		areas <- rep(tilearea, nx * ny)

	} else {
                # convert to mask
                win <- as.mask(win)

                # extract pixel coordinates inside window
		xx <- as.vector(raster.x(win)[win$m])
		yy <- as.vector(raster.y(win)[win$m])
                                
		# classify all pixels into tiles
		pixelid <- gridindex(xx, yy, 
				win$xrange, win$yrange, nx, nx)$index
                pixelid <- factor(pixelid, levels=seq(nx * ny))
                                
		# compute digital areas of tiles
		tilepixels <- table(pixelid)
		pixelarea <- win$xstep * win$ystep
		areas <- tilepixels * pixelarea

	} 

	# compute counting weights 
	w <- countingweights(id, areas)

	w
}


dirichlet.weights <- function(X, window = NULL, exact=T) {
	#
	# Compute weights based on Dirichlet tessellation of the window 
	# induced by the point pattern X. 
	# The weights are just the tile areas.
	#
	# NOTE:	X should contain both data and dummy points,
	# if you need these weights for the B-T-B method.
	#
	# Arguments X and (optionally) 'window' are interpreted as a
	# point pattern.
	#
	# If the window is a rectangle, we invoke Rolf Turner's "deldir"
	# package to compute the areas of the tiles of the Dirichlet
	# tessellation of the window frame induced by the points.
	# [NOTE: the functionality of deldir to create dummy points
	# is NOT used. ]
	#	if exact=T	compute the exact areas, using "deldir"
	#	if exact=F      compute the digital areas using exactdt()
	# 
	# If the window is a mask, we compute the digital area of
	# each tile of the Dirichlet tessellation by counting pixels.
	#
	#
	# 
	#
	
	X <- as.ppp(X, window)
	x <- X$x
	y <- X$y
	win <- X$window

        if(exact && !exists("deldir")) {
          warning("\'deldir\' package not found; using discrete approximation")
          exact <- F
        }

	if(exact && (win$type == "rectangle")) {
		rw <- c(win$xrange, win$yrange)
	        # invoke deldir() with NO DUMMY POINTS
		tessellation <- deldir(x, y, dpl=NULL, rw=rw)
	        # extract tile areas
	        w <- tessellation$summary[, 'dir.area']
		return(w)
	} else {
		# Compute digital areas of Dirichlet tiles.
                win <- as.mask(win)
                X$window <- win
		#
                # Nearest data point to each pixel:
                tileid <- exactdt(X)$i
                # 
		if(win$type == "mask") 
			# Restrict to window (result is a vector - OK)
			tileid <- tileid[win$m]
		# Count pixels in each tile
		id <- factor(tileid, levels=seq(X$n))
		counts <- table(id)
                # turn off the christmas lights
                class(counts) <- NULL
                names(counts) <- NULL
                dimnames(counts) <- NULL
		# Convert to digital area
		pixelarea <- win$xstep * win$ystep
		w <- pixelarea * counts
		# Check for zero pixel counts
		zeroes <- (counts == 0)
		if(any(zeroes)) {
			warning("some Dirichlet tiles have zero digital area")
			attr(w, "zeroes") <- zeroes
		}
		return(w)
	} 
}


#
#	window.S
#
#	A class 'owin' to define the "observation window"
#
#	$Revision: 4.10 $	$Date: 2002/01/17 10:15:06 $
#
#
#	A window may be either
#
#		- rectangular:
#                       a rectangle in R^2
#                       (with sides parallel to the coordinate axes)
#
#		- polygonal:
#			delineated by one or more non-self-intersecting
#                       polygons, possibly including polygonal holes.
#	
#		- digital mask:
#			defined by a binary image
#			whose pixel values are T wherever the pixel
#                       is inside the window
#
#	Any window is an object of class 'owin', 
#       containing at least the following entries:	
#
#		$type:	a string ("rectangle", "polygonal" or  "mask")
#
#		$xrange   
#		$yrange
#			vectors of length 2 giving the real dimensions 
#			of the enclosing box.
#
#	The 'rectangle' type has only these entries.
#
#       The 'polygonal' type has an additional entry
#
#               $bdry
#                       a list of polygons.
#                       Each entry bdry[[i]] determines a closed polygon.
#
#                       bdry[[i]] has components $x and $y which are
#                       the cartesian coordinates of the vertices of
#                       the i-th boundary polygon (without repetition of
#                       the first vertex, i.e. same convention as in the
#                       plotting function polygon().)
#
#
#	The 'mask' type has entries
#
#		$m		logical matrix
#		$dim		its dimension array
#		$xstep,ystep	x and y dimensions of a pixel
#		$xcol	        vector of x values for each column
#               $yrow           vector of y values for each row
#	
#	(the row index corresponds to increasing y coordinate; 
#	 the column index "   "     "   "  "  "  x "   "    ".)
#
#
#-----------------------------------------------------------------------------
#
owin <- function(xrange=c(0,1), yrange=c(0,1), poly=NULL, mask=NULL) {

  ## Exterminate ambiguities
  if(!missing(poly) && !is.null(poly) && !missing(mask) && !is.null(mask))
     stop("Ambiguous -- both polygonal boundary and digital mask supplied")
     
  if(missing(xrange) != missing(yrange))
    stop("If one of xrange, yrange is specified then both must be.")

  if(missing(poly) && missing(mask)) {
    ######### rectangle #################
    if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] <= xrange[1])
      stop("xrange should be a vector of length 2 giving (xmin, xmax)")
    if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] <= yrange[1])
      stop("yrange should be a vector of length 2 giving (ymin, ymax)")
    w <- list(type="rectangle", xrange=xrange, yrange=yrange)
    class(w) <- "owin"
    return(w)
  } else if(!missing(poly)) {
    ######### polygonal boundary ########
    #
    # test whether it's a single polygon or multiple polygons
    if(verify.xypolygon(poly, fatal=F))
      psingle <- T
    else if(all(unlist(lapply(poly, verify.xypolygon, fatal=F))))
      psingle <- F
    else
      stop("poly must be either a list(x,y) or a list of list(x,y)")
                  
    if(psingle) {
      # single boundary polygon
      if(area.xypolygon(poly) < 0)
        stop("Area of polygon is negative - maybe traversed in wrong direction?")
      bdry <- list(poly)
    } else {
      # multiple boundary polygons
      bdry <- poly
      if(sum(unlist(lapply(poly, area.xypolygon))) < 0)
        stop(paste("Area of window is negative;\n",
             "check that all polygons were traversed in the right direction"))
    }

    actual.xrange <- range(unlist(lapply(bdry, function(a) a$x)))
    if(missing(xrange))
      xrange <- actual.xrange
    else {
      if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] <= xrange[1])
        stop("xrange should be a vector of length 2 giving (xmin, xmax)")
      if(!all(xrange == range(c(xrange, actual.xrange))))
        stop("polygon's x coordinates outside xrange")
    }
    
    actual.yrange <- range(unlist(lapply(bdry, function(a) a$y)))
    if(missing(yrange))
      yrange <- actual.yrange
    else {
      if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] <= yrange[1])
        stop("yrange should be a vector of length 2 giving (ymin, ymax)")
      if(!all(yrange == range(c(yrange, actual.yrange))))
      stop("polygon's y coordinates outside yrange")
    }

    w <- list(type="polygonal", xrange=xrange, yrange=yrange, bdry=bdry)
    class(w) <- "owin"
    return(w)
    
  } else if(!missing(mask)) {
    ######### digital mask #####################
    
    if(!is.matrix(mask))
      stop("\`mask\' must be a matrix")
    if(!is.logical(mask))
      stop("The entries of \`mask\' must be logical")
    
    nc <- ncol(mask)
    nr <- nrow(mask)

    if(missing(xrange) && missing(yrange)) {
      # take pixels to be 1 x 1 unit
      xrange <- c(0,nc)
      yrange <- c(0,nr)
    } else {
      if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] <= xrange[1])
        stop("xrange should be a vector of length 2 giving (xmin, xmax)")
      if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] <= yrange[1])
        stop("yrange should be a vector of length 2 giving (ymin, ymax)")
    }

    xstep <- diff(xrange)/nc
    ystep <- diff(yrange)/nr
  
    out <- list(type     = "mask",
                xrange   = xrange,
                yrange   = yrange,
                dim      = c(nr, nc),
                xstep    = xstep,
                ystep    = ystep,
                warnings = c(
"Row index corresponds to increasing y coordinate; column to increasing x",
"Transpose matrices to get the standard presentation in S",
"Example: image(result$xcol,result$yrow,t(result$d))"
                ),
                xcol    = seq(xrange[1]+xstep/2, xrange[2]-xstep/2, length=nc),
                yrow    = seq(yrange[1]+ystep/2, yrange[2]-ystep/2, length=nr),
                m       = mask)
    class(out) <- "owin"
    return(out)
  }
  # never reached
  NULL
}

#
#-----------------------------------------------------------------------------
#

as.owin <- function(W) {
	# Tries to interpret data as an object of class 'window'
	# W may be
	#	an object of class 'window'
	#	a structure with entries xrange, yrange
	#	a four-element vector (interpreted xmin, xmax, ymin, ymax)
	#	a structure with entries xl, xu, yl, yu
	#	an object of class 'ppp'

	if(verifyclass(W, "owin", fatal=F))
		return(W)
	else if(checkfields(W, c("xrange", "yrange"))) {
		Z <- owin(W$xrange, W$yrange)
		return(Z)
	} else if(is.vector(W) && is.numeric(W) && length(W) == 4) {
		Z <- owin(W[1:2], W[3:4])
		return(Z)
	} else if(checkfields(W, c("xl", "xu", "yl", "yu"))) {
		Z <- owin(c(X$xl, X$xu),c(X$yl, X$yu))
		return(Z)
	} else if(verifyclass(W, "ppp", fatal=F))
		return(W$window)
	else
		stop("Can't interpret W as a window")
}		

#
#-----------------------------------------------------------------------------
#
#
as.rectangle <- function(...) {
        w <- as.owin(...)
        return(owin(w$xrange, w$yrange))
}

#
#-----------------------------------------------------------------------------
#
as.mask <- function(w, eps=NULL, dimyx=NULL) {
#	eps:		   grid mesh (pixel) size
#	dimyx:		   dimensions of pixel raster
#
#  NOTE: rows <=> y coordinate

        verifyclass(w, "owin")
        if(w$type == "mask")
          return(w)

        # determine row & column dimensions of output raster
        if(!is.null(dimyx)) {
          nr <- dimyx[1]
          nc <- dimyx[2]
        } else {
          # use pixel size 'eps'
          if(!is.null(eps)) {
            nr <- ceiling(diff(w$yrange)/eps)
            nc <- ceiling(diff(w$xrange)/eps)
          } else {
            # warning("raster resolution defaults to 1/100 window width")
            nr <- 100
            nc <- 100
          }
        }

        # mask with all entries True
        out <- owin(w$xrange, w$yrange,
                    mask=matrix(T, nrow=nr, ncol=nc))
        
        switch(w$type,
               rectangle={
                 return(out)
               },
               polygonal={
                 # test every pixel
                 x <- as.vector(raster.x(out))
                 y <- as.vector(raster.y(out))
                 value <- inside.owin(x, y, w)
                 return(owin(w$xrange, w$yrange,
                             mask=matrix(value, nrow=nr, ncol=nc)))
               },
               stop("Unrecognised window type \`", w$type, "\'")
               )
}
#
#-----------------------------------------------------------------------------
#
as.polygonal <- function(W) {
  # the only sensible use of this function is to
  # convert a rectangle to a polygon
  verifyclass(W, "owin")
  switch(W$type,
         rectangle = {
           xr <- W$xrange
           yr <- W$yrange
           return(owin(xr, yr, poly=list(x=xr[c(1,2,2,1)],y=yr[c(1,1,2,2)])))
         },
         polygonal = {
           return(W)
         },
         mask = {
           stop("A mask cannot be converted to a polygon")
         }
         )
}

#
# ----------------------------------------------------------------------

validate.mask <- function(w, fatal=T) {
  verifyclass(w, "owin", fatal=fatal)
  if(w$type == "mask")
    return(T)
  if(fatal)
      stop(deparse(substitute(w)), "is not a binary mask")
  else {
      warning(deparse(substitute(w)), "is not a binary mask")
      return(F)
  }
}
             
raster.x <- function(w) {
	validate.mask(w)
        di <- w$dim
        nr <- di[1]
        nc <- di[2]
	m <- matrix(0, nrow=nr, ncol=nc)
	matrix(w$xcol[col(m)], nrow=nr, ncol=nc)
}

raster.y <- function(w) {
	validate.mask(w)
        di <- w$dim
        nr <- di[1]
        nc <- di[2]
	m <- matrix(0, nrow=nr, ncol=nc)
	matrix(w$yrow[row(m)], nrow=nr, ncol=nc)
}

nearest.raster.point <- function(x,y,w) {
	validate.mask(w)
	nr <- w$dim[1]
	nc <- w$dim[2]
	cc <- round(0.5 + (x - w$xrange[1])/w$xstep)
	rr <- round(0.5 + (y - w$yrange[1])/w$ystep)
	cc <- pmax(1,pmin(cc, nc))
	rr <- pmax(1,pmin(rr, nr))
	return(list(row=rr, col=cc))
}

#------------------------------------------------------------------
		
bounding.box <- function(w) {
        # determine a tight bounding box for the window w
        verifyclass(w, "owin")

        switch(w$type,
               rectangle = {
                 return(w)
               },
               polygonal = {
                 bdry <- w$bdry
                 xr <- range(unlist(lapply(bdry, function(a) a$x)))
                 yr <- range(unlist(lapply(bdry, function(a) a$y)))
                 return(owin(xr, yr))
               },
               mask = {
                 m <- w$m
                 x <- raster.x(w)
                 y <- raster.y(w)
                 xr <- range(x[m])
                 yr <- range(y[m])
                 return(owin(xr, yr))
               },
               stop("unrecognised window type", w$type)
               )
}
  
complement.owin <- function(w) {
	verifyclass(w, "owin")
        switch(w$type,
               mask = {
                 w$m <- !(w$m)
               },
               polygonal = {

                 bdry <- w$bdry
                 
                 # bounding box, in anticlockwise order
                 box <- list(x=w$xrange[c(1,2,2,1)],
                             y=w$yrange[c(1,1,2,2)])
                 boxarea <- area.xypolygon(box)
                 
                 # first check whether one of the current boundary polygons
                 # is the bounding box itself (with + sign)
                 nvert <- unlist(lapply(bdry, function(a) { length(a$x) }))
                 area <- unlist(lapply(bdry, area.xypolygon, test01=F))
                 boxarea.mineps <- boxarea * (1 - .Machine$single.eps)
                 is.box <- (nvert == 4 & area >= boxarea.mineps)
                 if(sum(is.box) > 1)
                   stop("Internal error: multiple copies of bounding box")
                 
                 # if box is present (with + sign), remove it
                 if(any(is.box))
                   bdry <- bdry[!is.box]
                 
                 # reverse the direction of each polygon
                 bdry <- lapply(bdry, reverse.xypolygon)
                 
                 # if box was absent, add it
                 if(!any(is.box))
                   bdry <- c(bdry, list(box))   # sic
                 
                 # put back into w
                 w$bdry <- bdry
               },
               rectangle = {
                 stop("window is a rectangle - its complement is empty")
               },
               stop("unrecognised window type", w$type)
               )
	return(w)
}

#-----------------------------------------------------------

inside.owin <- function(x, y, w) {
  # test whether (x,y) is inside window w
  # x, y may be vectors 
  
  verifyclass(w, "owin")

  # test whether inside bounding rectangle
  xr <- w$xrange
  yr <- w$yrange
  frameok <- (xr[1] <= x) & (x <= xr[2]) & (yr[1] <= y) & (y <= yr[2])

  if(all(!frameok))  # all points OUTSIDE window - no further work needed
    return(frameok)

  ok <- frameok
  switch(w$type,
         rectangle = {
           return(ok)
         },
         polygonal = {
           xy <- list(x=x,y=y)
           bdry <- w$bdry
           total <- rep(0, length(x))
           on.bdry <- rep(F, length(x))
           for(i in seq(bdry)) {
             score <- inside.xypolygon(xy, bdry[[i]], test01=F)
             total <- total + score
             on.bdry <- on.bdry | attr(score, "on.boundary")
           }
           # any points identified as belonging to the boundary get score 1
           total[on.bdry] <- 1
           # check for sanity now..
           if(any(total * (1-total) != 0)) 
             stop("internal error: some total scores are neither 0 nor 1")
           return(ok & (total != 0))
         },
         mask = {
           # consider only those points which are inside the frame
           xf <- x[frameok]
           yf <- y[frameok]
           # map locations to raster (row,col) coordinates
           loc <- nearest.raster.point(xf,yf,w)
           # look up mask values
           mas <- w$m
           nf <- sum(frameok)
           okf <- logical(nf)
           for(i in 1:nf) 
             okf[i] <- mas[loc$row[i],loc$col[i]]
           # insert into 'ok' vector
           ok[frameok] <- okf
           return(ok)
         },
         stop("unrecognised window type", w$type)
         )
}
#
#	wingeom.S	Various geometrical computations in windows
#
#
#	$Revision: 4.5 $	$Date: 2002/01/17 10:43:13 $
#
#
#
#
#-------------------------------------
area.owin <- function(w) {
	verifyclass(w, "owin")
        switch(w$type,
               rectangle = {
		width <- abs(diff(w$xrange))
		height <- abs(diff(w$yrange))
		area <- width * height
               },
               polygonal = {
                 area <- sum(unlist(lapply(w$bdry, area.xypolygon)))
               },
               mask = {
                 pixelarea <- abs(w$xstep * w$ystep)
                 npixels <- sum(w$m)
                 area <- pixelarea * npixels
               },
               stop("Unrecognised window type")
        )
        return(area)
}

eroded.areas <- function(w, r) {
	verifyclass(w, "owin")
	
	switch(w$type,
               rectangle = {
                 width <- abs(diff(w$xrange))
                 height <- abs(diff(w$yrange))
                 areas <- pmax(width - 2 * r, 0) * pmax(height - 2 * r, 0)
               },
               polygonal = {
                 # warning("Approximating polygonal window by digital image")
                 w <- as.mask(w)
                 areas <- eroded.areas(w, r)
               },
               mask = {
                 # distances from each pixel to window boundary
                 b <- bdist.pixels(w, coords=F)
                 # histogram breaks to satisfy hist()
                 Bmax <- max(b, r)
                 breaks <- c(-1,r,Bmax+1)
                 # histogram of boundary distances
                 h <- hist(b, breaks=breaks, plot=F, probability=F)$counts
                 # reverse cumulative histogram
                 H <- rev(cumsum(rev(h)))
                 # drop first entry corresponding to r=-1
                 H <- H[-1]
                 # convert count to area
                 pixarea <- w$xstep * w$ystep
                 areas <- pixarea * H
               },
 	       stop("unrecognised window type")
               )
	areas
}	

diameter <- function(w) {
	verifyclass(w, "owin")
	
        width <- abs(diff(w$xrange))
        height <- abs(diff(w$yrange))
        
        sqrt(width^2 + height^2)
}

even.breaks.owin <- function(w) {
	verifyclass(w, "owin")
        Rmax <- diameter(w)
        make.even.breaks(Rmax, Rmax/(100 * sqrt(2)))
}

unit.square <- function() { owin(c(0,1),c(0,1)) }

overlap.owin <- function(A, B) {
  # compute the area of overlap between two windows
  At <- A$type
  Bt <- B$type
  if(At=="rectangle" && Bt=="rectangle") {
    xmin <- max(A$xrange[1],B$xrange[1])
    xmax <- min(A$xrange[2],B$xrange[2])
    if(xmax <= xmin) return(0)
    ymin <- max(A$yrange[1],B$yrange[1])
    ymax <- min(A$yrange[2],B$yrange[2])
    if(ymax <= ymin) return(0)
    return((xmax-xmin) * (ymax-ymin))
  }
  if((At=="rectangle" && Bt=="polygonal")
     || (At=="polygonal" && Bt=="rectangle")
     || (At=="polygonal" && Bt=="polygonal"))
  {
    AA <- as.polygonal(A)$bdry
    BB <- as.polygonal(B)$bdry
    area <- 0
    for(i in seq(AA))
      for(j in seq(BB))
        area <- area + overlap.xypolygon(AA[[i]], BB[[j]])
    return(area)
  }
  if(At=="mask") {
    # count pixels in A that belong to B
    pixelarea <- abs(A$xstep * A$ystep)
    x <- as.vector(raster.x(A))
    y <- as.vector(raster.y(A))
    ok <- inside.owin(x, y, B)
    return(pixelarea * sum(ok))
  }
  if(Bt== "mask") {
    # count pixels in B that belong to A
    pixelarea <- abs(B$xstep * B$ystep)
    x <- as.vector(raster.x(B))
    y <- as.vector(raster.y(B))
    ok <- inside.owin(x, y, A)
    return(pixelarea * sum(ok))
  }
  stop("Internal error")
}
  


  
#
#    xypolygon.S
#
#    $Revision: 1.3 $    $Date: 2002/01/17 10:01:49 $
#
#    low-level functions defined for polygons in list(x,y) format
#
verify.xypolygon <- function(p, fatal=T) {
  whinge <- NULL
  if(!is.list(p) || length(p) != 2 ||
     length(names(p)) != 2 || any(sort(names(p)) != c("x","y")))
    whinge <- "polygon must be a list with two components x and y"
  else if(is.null(p$x) || is.null(p$y) || !is.numeric(p$x) || !is.numeric(p$y))
    whinge <- "components x and y must be numeric vectors"
  else if(length(p$x) != length(p$y))
    whinge <- "lengths of x and y vectors unequal"
  ok <- is.null(whinge)
  if(!ok && fatal)
    stop(whinge)
  return(ok)
}

inside.xypolygon <- function(pts, polly, test01=T) {
  # pts:  list(x,y) points to be tested
  # polly: list(x,y) vertices of a single polygon (n joins to 1)
  # test01: logical - if TRUE, test whether all values in output are 0 or 1
  verify.xypolygon(pts)
  verify.xypolygon(polly)
  
  x <- pts$x
  y <- pts$y
  xp <- polly$x
  yp <- polly$y

  npts <- length(x)
  nedges <- length(xp)   # sic

  score <- rep(0, npts)
  on.boundary <- rep(F, npts)

  for(i in 1:nedges) {
    x0 <- xp[i]
    y0 <- yp[i]
    x1 <- if(i == nedges) xp[1] else xp[i+1]
    y1 <- if(i == nedges) yp[1] else yp[i+1]
    dx <- x1 - x0
    dy <- y1 - y0
    if(dx < 0) {
      # upper edge
      xcriterion <- (x - x0) * (x - x1)
      consider <- (xcriterion <= 0)
      if(any(consider)) {
        ycriterion <- y[consider] * dx - x[consider] * dy + (x0 * dy - y0 * dx)
        # closed inequality
        contrib <- (ycriterion >= 0) * ifelse(xcriterion[consider] == 0, 1/2, 1)
        # positive edge sign
        score[consider] <- score[consider] + contrib
        # detect whether any point lies on this segment
        on.boundary[consider] <- on.boundary[consider] | (ycriterion == 0)
      }
    } else if(dx > 0) {
      # lower edge
      xcriterion <- (x - x0) * (x - x1)
      consider <- (xcriterion <= 0)
      if(any(consider)) {
        ycriterion <- y[consider] * dx - x[consider] * dy + (x0 * dy - y0 * dx)
        # open inequality
        contrib <- (ycriterion < 0) * ifelse(xcriterion[consider] == 0, 1/2, 1)
        # negative edge sign
        score[consider] <- score[consider] - contrib
        # detect whether any point lies on this segment
        on.boundary[consider] <- on.boundary[consider] | (ycriterion == 0)
      }
    } else {
      # vertical edge
      consider <- (x == x0)
      if(any(consider)) {
        # zero score
        # detect whether any point lies on this segment
        yconsider <- y[consider]
        ycriterion <- (yconsider - y0) * (yconsider - y1)
        on.boundary[consider] <- on.boundary[consider] | (ycriterion <= 0)
      }
    }
  }
  
  # any point recognised as lying on the boundary gets score 1.
  score[on.boundary] <- 1

  if(test01) {
    # check sanity
    if(!all((score == 0) | (score == 1)))
      warning("internal error: some scores are not equal to 0 or 1")
  }

  attr(score, "on.boundary") <- on.boundary
  
  return(score)
}

area.xypolygon <- function(polly) {
  #
  # polly: list(x,y) vertices of a single polygon (n joins to 1)
  #
  verify.xypolygon(polly)
  
  xp <- polly$x
  yp <- polly$y
  
  nedges <- length(xp)   # sic
  
  # place x axis below polygon
  yp <- yp - min(yp) 

  # join vertex n to vertex 1
  nxt <- c(2:nedges, 1)

  # x step, WITH sign
  dx <- xp[nxt] - xp

  # average height 
  ym <- (yp + yp[nxt])/2
  
  -sum(dx * ym)
}

bdrylength.xypolygon <- function(polly) {
  verify.xypolygon(polly)
  xp <- polly$x
  yp <- polly$y
  nedges <- length(xp)
  nxt <- c(2:nedges, 1)
  dx <- xp[nxt] - xp
  dy <- yp[nxt] - yp
  sum(sqrt(dx^2 + dy^2))
}

reverse.xypolygon <- function(p) {
  # reverse the order of vertices
  # (=> change sign of polygon)

  verify.xypolygon(p)
  
  n <- length(p$x)

  return(list(x=p$x[n:1], y=p$y[n:1]))
}

overlap.xypolygon <- function(P, Q) {
  # compute area of overlap of two simple closed polygons 
  verify.xypolygon(P)
  verify.xypolygon(Q)
  
  xp <- P$x
  yp <- P$y
  np <- length(xp)
  nextp <- c(2:np, 1)

  xq <- Q$x
  yq <- Q$y
  nq <- length(xq)
  nextq <- c(2:nq, 1)

  # adjust y coordinates so all are nonnegative
  ylow <- min(c(yp,yq))
  yp <- yp - ylow
  yq <- yq - ylow

  area <- 0
  for(i in 1:np) {
    ii <- c(i, nextp[i])
    xpii <- xp[ii]
    ypii <- yp[ii]
    for(j in 1:nq) {
      jj <- c(j, nextq[j])
      area <- area +
        overlap.trapezium(xpii, ypii, xq[jj], yq[jj])
    }
  }
  return(area)
}

overlap.trapezium <- function(xa, ya, xb, yb, verb=F) {
  # compute area of overlap of two trapezia
  # which have same baseline y = 0
  #
  # first trapezium has vertices
  # (xa[1], 0), (xa[1], ya[1]), (xa[2], ya[2]), (xa[2], 0).
  # Similarly for second trapezium
  
  # Test for vertical edges
  dxa <- diff(xa)
  dxb <- diff(xb)
  if(dxa == 0 || dxb == 0)
    return(0)

  # Order x coordinates, x0 < x1
  if(dxa > 0) {
    signa <- 1
    lefta <- 1
    righta <- 2
    if(verb) cat("A is positive\n")
  } else {
    signa <- -1
    lefta <- 2
    righta <- 1
    if(verb) cat("A is negative\n")
  }
  if(dxb > 0) {
    signb <- 1
    leftb <- 1
    rightb <- 2
    if(verb) cat("B is positive\n")
  } else {
    signb <- -1
    leftb <- 2
    rightb <- 1
    if(verb) cat("B is negative\n")
  }
  signfactor <- signa * signb # actually (-signa) * (-signb)
  if(verb) cat(paste("sign factor =", signfactor, "\n"))

  # Intersect x ranges
  x0 <- max(xa[lefta], xb[leftb])
  x1 <- min(xa[righta], xb[rightb])
  if(x0 >= x1)
    return(0)
  if(verb) {
    cat(paste("Intersection of x ranges: [", x0, ",", x1, "]\n"))
    abline(v=x0, lty=3)
    abline(v=x1, lty=3)
  }

  # Compute associated y coordinates
  slopea <- diff(ya)/diff(xa)
  y0a <- ya[lefta] + slopea * (x0-xa[lefta])
  y1a <- ya[lefta] + slopea * (x1-xa[lefta])
  slopeb <- diff(yb)/diff(xb)
  y0b <- yb[leftb] + slopeb * (x0-xb[leftb])
  y1b <- yb[leftb] + slopeb * (x1-xb[leftb])
  
  # Determine whether upper edges intersect
  # if not, intersection is a single trapezium
  # if so, intersection is a union of two trapezia

  yd0 <- y0b - y0a
  yd1 <- y1b - y1a
  if(yd0 * yd1 >= 0) {
    # edges do not intersect
    areaT <- (x1 - x0) * (min(y1a,y1b) + min(y0a,y0b))/2
    if(verb) cat(paste("Edges do not intersect\n"))
  } else {
    # edges do intersect
    # find intersection
    xint <- x0 + (x1-x0) * abs(yd0/(yd1 - yd0))
    yint <- y0a + slopea * (xint - x0)
    if(verb) {
      cat(paste("Edges intersect at (", xint, ",", yint, ")\n"))
      points(xint, yint, cex=2, pch="O")
    }
    # evaluate left trapezium
    left <- (xint - x0) * (min(y0a, y0b) + yint)/2
    # evaluate right trapezium
    right <- (x1 - xint) * (min(y1a, y1b) + yint)/2
    areaT <- left + right
    if(verb)
      cat(paste("Left area = ", left, ", right=", right, "\n"))    
  }

  # return area of intersection multiplied by signs 
  return(signfactor * areaT)
}

#
#      xysegment.S
#
#     $Revision: 1.2 $    $Date: 2001/11/19 05:08:06 $
#
# Low level utilities for analytic geometry for line segments
#
# author: Adrian Baddeley 2001
#         from an original by Rob Foxall 1997
#
# distpl(p, l) 
#       distance from a single point p  = (xp, yp)
#       to a single line segment l = (x1, y1, x2, y2)
#
# distppl(p, l) 
#       distances from each of a list of points p[i,]
#       to a single line segment l = (x1, y1, x2, y2)
#       [uses only vector parallel ops]
#
# distppll(p, l) 
#       distances from each of a list of points p[i,]
#       to each of a list of line segments l[i,] 
#       [uses large matrices and 'outer()']
#

distpl <- function(p, l) {
  xp <- p[1]
  yp <- p[2]
  dx <- l[3]-l[1]
  dy <- l[4]-l[2]
  leng <- sqrt(dx^2 + dy^2)
  # vector from 1st endpoint to p
  xpl <- xp - l[1]
  ypl <- yp - l[2]
  # distance from p to 1st & 2nd endpoints
  d1 <- sqrt(xpl^2 + ypl^2)
  d2 <- sqrt((xp-l[3])^2 + (yp-l[4])^2)
  dmin <- min(d1,d2)
  # rotation sine & cosine
  co <- dx/leng
  si <- dy/leng
  # back-rotated coords of p
  xpr <- co * xpl + si * ypl
  ypr <-  - si * xpl + co * ypl
  # test
  if(xpr >= 0 && xpr <= leng)
    dmin <- min(dmin, abs(ypr))
  return(dmin)
}

distppl <- function(p, l) {
  xp <- p[,1]
  yp <- p[,2]
  dx <- l[3]-l[1]
  dy <- l[4]-l[2]
  leng <- sqrt(dx^2 + dy^2)
  # vector from 1st endpoint to p
  xpl <- xp - l[1]
  ypl <- yp - l[2]
  # distance from p to 1st & 2nd endpoints
  d1 <- sqrt(xpl^2 + ypl^2)
  d2 <- sqrt((xp-l[3])^2 + (yp-l[4])^2)
  dmin <- pmin(d1,d2)
  # rotation sine & cosine
  co <- dx/leng
  si <- dy/leng
  # back-rotated coords of p
  xpr <- co * xpl + si * ypl
  ypr <-  - si * xpl + co * ypl
  # ypr is perpendicular distance to infinite line
  # Applies only when xp, yp in the middle
  middle <- (xpr >= 0 & xpr <= leng)
  if(any(middle))
    dmin[middle] <- pmin(dmin[middle], abs(ypr[middle]))
  
  return(dmin)
}

distppll <- function(p, l) {
  np <- nrow(p)
  nl <- nrow(l)
  xp <- p[,1]
  yp <- p[,2]
  dx <- l[,3]-l[,1]
  dy <- l[,4]-l[,2]
  # segment lengths
  leng <- sqrt(dx^2 + dy^2)
  # rotation sines & cosines
  co <- dx/leng
  si <- dy/leng
  co <- matrix(co, nrow=np, ncol=nl, byrow=T)
  si <- matrix(si, nrow=np, ncol=nl, byrow=T)
  # matrix of squared distances from p[i] to 1st endpoint of segment j
  xp.x1 <- outer(xp, l[,1], "-")
  yp.y1 <- outer(yp, l[,2], "-")
  d1 <- xp.x1^2 + yp.y1^2
  # ditto for 2nd endpoint
  xp.x2 <- outer(xp, l[,3], "-")
  yp.y2 <- outer(yp, l[,4], "-")
  d2 <- xp.x2^2 + yp.y2^2
  # for each (i,j) rotate p[i] around 1st endpoint of segment j
  # so that line segment coincides with x axis
  xpr <- xp.x1 * co + yp.y1 * si
  ypr <-  - xp.x1 * si + yp.y1 * co
  d3 <- ypr^2
  # test
  lenf <- matrix(leng, nrow=np, ncol=nl, byrow=T)
  outside <- (xpr < 0 | xpr > lenf)
  if(any(outside))
    d3[outside] <- Inf

  dmin <- matrix(pmin(d1, d2, d3),nrow=np, ncol=nl)
  return(sqrt(dmin))
}



