.packageName <- "spdep"
# Copyright 2002-2003 by Roger Bivand and Marilia Carvalho
#

EBImoran <- function (z, listw, nn, S0, zero.policy = FALSE) 
{
    zm <- mean(z)
    zz <- sum((z - zm)^2)
    lz <- lag.listw(listw, z, zero.policy = zero.policy)
    EBI <- (nn/S0) * ((t(z) %*% lz)/zz)
    res <- EBI
    res
}
EBImoran.mc <- function (n, x, listw, nsim, zero.policy = FALSE,
 alternative = "greater", spChk = NULL) 
{
    alternative <- match.arg(alternative, c("greater", "less"))
    if (!inherits(listw, "listw")) 
        stop(paste(deparse(substitute(listw)), "is not a listw object"))
    if (missing(nsim)) 
        stop("nsim must be given")
    m <- length(listw$neighbours)
    if (m != length(x)) 
        stop("objects of different length")
    if (m != length(n)) 
        stop("objects of different length")
    if (is.null(spChk)) 
        spChk <- get.spChkOption()
    if (spChk && !chkIDs(x, listw)) 
        stop("Check of data and weights ID integrity failed")
    if (spChk && !chkIDs(n, listw)) 
        stop("Check of data and weights ID integrity failed")
    if (nsim > gamma(m + 1)) 
        stop("nsim too large for this number of observations")
    S0 <- Szero(listw)
    EB <- EBest(n, x)
    p <- EB$raw
    b <- attr(EB, "parameters")$b
    a <- attr(EB, "parameters")$a
    v <- a + (b/x)
    v[v < 0] <- b/x
    z <- (p - b)/sqrt(v)
    res <- numeric(length = nsim + 1)
    for (i in 1:nsim) res[i] <- EBImoran(sample(z), listw, m, 
        S0, zero.policy)
    res[nsim + 1] <- EBImoran(z, listw, m, S0, zero.policy)
    rankres <- rank(res)
    zrank <- rankres[length(res)]
    diff <- nsim - zrank
    diff <- ifelse(diff > 0, diff, 0)
    if (alternative == "less") 
        pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
    else if (alternative == "greater") 
        pval <- punif((diff + 1)/(nsim + 1))
    if (pval < 0 || pval > 1) 
	warning("Out-of-range p-value: reconsider test arguments")
    statistic <- res[nsim + 1]
    names(statistic) <- "statistic"
    parameter <- zrank
    names(parameter) <- "observed rank"
    method <- "Monte-Carlo simulation of Empirical Bayes Index"
    data.name <- paste("cases: ", deparse(substitute(n)), ", risk population: ", 
        deparse(substitute(x)), "\nweights: ", deparse(substitute(listw)), 
        "\nnumber of simulations + 1: ", nsim + 1, "\n", sep = "")
    lres <- list(statistic = statistic, parameter = parameter, 
        p.value = pval, alternative = alternative, method = method, 
        data.name = data.name, res = res, z = z)
    class(lres) <- c("htest", "mc.sim")
    lres
}

probmap <- function(n, x, row.names=NULL) {
    if (!is.numeric(x)) 
        stop(paste(deparse(substitute(x)), "is not a numeric vector"))
    if (!is.numeric(n)) 
        stop(paste(deparse(substitute(n)), "is not a numeric vector"))
    if (any(is.na(x))) 
        stop("NA in at risk population")
    if (any(is.na(n))) 
        stop("NA in cases")
    if (any(x < 0)) 
        stop("negative risk population")
    if (any(n < 0)) 
        stop("negative number of cases")
    p <- n/x
    nsum <- sum(n)
    xsum <- sum(x)
    b <- nsum/xsum
    expCount <- x*b
    relRisk <- 100*(n/expCount)
    pmap <- ppois(n, expCount)
    if (is.null(row.names)) 
	res <- data.frame(raw=p, expCount=expCount, relRisk=relRisk, 
	pmap=pmap)
    else
    	res <- data.frame(raw=p, expCount=expCount, relRisk=relRisk, 
	pmap=pmap, row.names=row.names)
    res
}


EBest <- function(n, x) {
    if (!is.numeric(x)) 
        stop(paste(deparse(substitute(x)), "is not a numeric vector"))
    if (!is.numeric(n)) 
        stop(paste(deparse(substitute(n)), "is not a numeric vector"))
    if (any(is.na(x))) 
        stop("NA in at risk population")
    if (any(is.na(n))) 
        stop("NA in cases")
    if (any(x < .Machine$double.eps)) 
        stop("non-positive risk population")
    if (any(n < 0)) 
        stop("negative number of cases")
    m <- length(n)
    p <- n/x
    nsum <- sum(n)
    xsum <- sum(x)
    b <- nsum/xsum
    s2 <- sum(x * (((p - b)^2)/xsum))
    a <- s2 - (b/(xsum/m))
    if (a < 0) a <- 0
    est <- b + (a*(p - b)) / (a + (b/x))
    res <- data.frame(raw=p, estmm=est)
    attr(res, "parameters") <- list(a=a, b=b)
    res
}

EBlocal <- function(ri, ni, nb, zero.policy = FALSE,
    spChk = NULL, geoda = FALSE) {
    if (class(nb) != "nb") 
        stop(paste(deparse(substitute(nb)), "is not an nb object"))
    lnb <- length(nb)
    if (lnb != length(ri)) 
        stop("objects of different length")
    if (lnb != length(ni)) 
        stop("objects of different length")
    if (is.null(spChk)) 
        spChk <- get.spChkOption()
    if (spChk && !chkIDs(ni, nb)) 
        stop("Check of data and neighbour ID integrity failed")
    if (spChk && !chkIDs(ri, nb)) 
        stop("Check of data and neighbour ID integrity failed")
    if (!is.numeric(ni)) 
        stop(paste(deparse(substitute(ni)), "is not a numeric vector"))
    if (!is.numeric(ri)) 
        stop(paste(deparse(substitute(ri)), "is not a numeric vector"))
    if (any(is.na(ni))) 
        stop("NA in at risk population")
    if (any(is.na(ri))) 
        stop("NA in cases")
    if (any(ni < 0)) 
        stop("negative risk population")
    if (any(ri < 0)) 
        stop("negative number of cases")
    lw <- nb2listw(include.self(nb), style="B", zero.policy=zero.policy)
    xi <- ri/ni
    r.i <- lag.listw(lw, ri, zero.policy = zero.policy)
    n.i <- lag.listw(lw, ni, zero.policy = zero.policy)
    nbar.i <- lag.listw(nb2listw(include.self(nb), style="W",
        zero.policy=zero.policy), ni, zero.policy = zero.policy)

    m.i <- r.i/n.i
    if (geoda) {
	work <- vector(mode="list", length=lnb)
	cnb <- card(lw$neighbours)
	for (i in 1:lnb) {
		work[[i]] <- rep(m.i[i], cnb[i])
		work[[i]] <- xi[lw$neighbours[[i]]] - work[[i]]
		work[[i]] <- work[[i]]^2
		work[[i]] <- ni[lw$neighbours[[i]]] * work[[i]]
	}
	C.i <- sapply(work, sum)
    } else {
        C.i <- lag.listw(lw, (ni * (xi - m.i)^2), zero.policy = zero.policy)
    }
    a.i <- (C.i/n.i) - (m.i/nbar.i)
    a.i[a.i < 0] <- 0
    est <- m.i + (xi - m.i) * (a.i / (a.i + (m.i/ni)))

    res <- data.frame(raw=xi, est=est)
    attr(res, "parameters") <- list(a=a.i, m=m.i)
    res
}
# Copyright 2002-2003 by Roger Bivand 
#

LR.sarlm <- function(x, y)
{
	if (!inherits(x, "logLik")) LLx <- logLik(x)
	else LLx <- x
	if (!inherits(y, "logLik")) LLy <- logLik(y)
	else LLy <- y
	statistic <- 2*(LLx - LLy)
	attr(statistic, "names") <- "Likelihood ratio"
	parameter <- abs(attr(LLx, "df") - attr(LLy, "df"))
	if (parameter < 1) 
		stop("non-positive degrees of freedom: no test possible")
	attr(parameter, "names") <- "df"
	p.value <- 1 - pchisq(abs(statistic), parameter)
	estimate <- c(LLx, LLy)
	attr(estimate, "names") <- c(paste("Log likelihood of",
		deparse(substitute(x))), paste("Log likelihood of",
		deparse(substitute(y))))
	method <- "Likelihood ratio for spatial linear models"
	res <- list(statistic=statistic, parameter=parameter,
		p.value=p.value, estimate=estimate, method=method)
	class(res) <- "htest"
	res
}

logLik.sarlm <- function(object, ...) {
	LL <- c(object$LL)
	class(LL) <- "logLik"
	N <- length(residuals(object))
	attr(LL, "nall") <- N
	attr(LL, "nobs") <- N
	attr(LL, "df") <- object$parameters
	LL
}

#residuals.sarlm <- function(object, ...) return(object$residuals)

# Copyright 2001 by Roger Bivand 
#

airdist <- function(ann=FALSE) {
	usr <- diff(par("usr"))[c(1,3)]
	plt <- diff(par("plt"))[c(1,3)]
	if (abs(diff(plt/usr)) > 0.005)
		warning("plot x and y scales may differ: use plot(..., asp=1)")
	coords <- locator(2)
	res <- sqrt(diff(coords$x)^2 + diff(coords$y)^2)
	if (ann) {
		lines(coords)
		text(mean(coords$x), mean(coords$y), format(res, digits=3),
			pos=4, offset=0.2, cex=0.7)
	}
	if (.Platform$OS.type == "windows") bringToTop(-1)
	print(res)
	invisible(list(dist=res, coords=coords))
}
# Copyright 2003 by Roger Bivand 
#

anova.sarlm <- function(object, ...) {
    if (length(list(object, ...)) > 1) {
        getResponseFormula <- function (object) 
        {
            form <- formula(object)
            if (!(inherits(form, "formula") && (length(form) == 3))) {
                stop("\"Form\" must be a two sided formula")
            }
            eval(parse(text = paste("~", deparse(form[[2]]))))
        }

        object <- list(object, ...)
        ancall <- sys.call()
        nmodels <- length(object)
        if (nmodels == 1) return(anova(object))
        termsClass <- unlist(lapply(object, data.class))
        if (!all(match(termsClass, c("lm", "sarlm"), 0))) {
            stop(paste("Objects must inherit from classes \"sarlm\" or \"lm\""))
        }
        resp <- unlist(lapply(object, 
	    function(el) deparse(getResponseFormula(el)[[2]])))
        subs <- as.logical(match(resp, resp[1], FALSE))
        if (!all(subs)) 
            warning(paste("Some fitted objects deleted because", 
                "response differs from the first model"))
        if (sum(subs) == 1) 
            stop("First model has a different response from the rest")
        object <- object[subs]
        aux <- lapply(object, logLik)
        if (length(unique(unlist(lapply(aux, 
	    function(el) length(residuals(el)))))) > 1) {
            stop("All fitted objects must use the same number of observations")
        }
        dfModel <- unlist(lapply(aux, function(el) attr(el, "df")))
        logLik <- unlist(lapply(aux, function(el) c(el)))
        AIC <- unlist(lapply(aux, AIC))
        aod <- data.frame(Model = (1:nmodels), df = dfModel, 
        AIC = AIC, logLik = logLik, check.names = FALSE)
        ddf <- diff(dfModel)
        if (sum(abs(ddf)) > 0) {
	   effects <- rep("", nmodels)
	   for (i in 2:nmodels) {
                if (ddf[i - 1] != 0) {
            	    effects[i] <- paste(i - 1, i, sep = " vs ")
                }
            }
            pval <- rep(NA, nmodels - 1)
            ldf <- as.logical(ddf)
            lratio <- 2 * abs(diff(logLik))
            lratio[!ldf] <- NA
            pval[ldf] <- 1 - pchisq(lratio[ldf], abs(ddf[ldf]))
	    aod <- data.frame(aod, Test = effects, L.Ratio = c(NA, 
                lratio), "p-value" = c(NA, pval), check.names = FALSE)
    	}
    	row.names(aod) <- unlist(lapply(as.list(ancall[-1]), 
            deparse))
    	attr(aod, "nmodels") <- nmodels
    	class(aod) <- c("anova.sarlm", "data.frame")
    	return(aod)

    } else {
    	if (!inherits(object, "sarlm")) 
            stop("object not a fitted simultaneous autoregressive model")
        LL <- logLik(object)
        AIC <- AIC(LL)
        res <- data.frame("AIC"=AIC, "Log likelihood"=LL, "df"=attr(LL, "df"),
	    row.names=deparse(substitute(object)))
        return(res)
    }
}


# Copyright 2004 by Roger Bivand (original taken from bptest() in the lmtest
# package, Copyright (C) 2001 Torsten Hothorn and Achim Zeileis and released
# under GNU General Public License, Version 2.
#

bptest.sarlm <- function (object, studentize = TRUE) 
{
    if(!inherits(object, "sarlm")) stop("not sarlm object")
    Z <- model.matrix(object$lm.target)
    k <- ncol(Z)
    n <- nrow(Z)
    resi <- residuals(object$lm.target)
    sigma2 <- sum(resi^2)/n
    if (studentize) {
        w <- resi^2 - sigma2
        fv <- lm.fit(Z, w)$fitted
        bp <- n * sum(fv^2)/sum(w^2)
        method <- "studentized Breusch-Pagan test"
    }
    else {
        f <- resi^2/sigma2 - 1
        fv <- lm.fit(Z, f)$fitted
        bp <- 0.5 * sum(fv^2)
        method <- "Breusch-Pagan test"
    }
    names(bp) <- "BP"
    df <- k - 1
    names(df) <- "df"
    RVAL <- list(statistic = bp, parameter = df, method = method, 
        p.value = 1 - pchisq(bp, df))
    class(RVAL) <- "htest"
    return(RVAL)
}

# Copyright 2001 by Roger Bivand 
#

card <- function(nb) {
    if (!inherits(nb, "nb")) stop("not a neighbours list")
    z <- .Call("card", nb, PACKAGE="spdep")
    invisible(z)
}
# Copyright 2001 by Roger Bivand 
#

rookcell <- function(rowcol, nrow, ncol, torus=FALSE, rmin=1, cmin=1) {
	if (is.null(dim(rowcol))) rowcol <- t(as.matrix(rowcol))
        if(nrow(rowcol) != 1) stop("only single grid cell handled")
	row <- rowcol[1]
	col <- rowcol[2]
	if (torus) {
		y <- c(ifelse(col-1 < cmin, ncol, col-1), col, col,
			ifelse(col+1 > (ncol+(cmin-1)), cmin, col+1))
		x <- c(row, ifelse(row-1 < rmin, nrow, row-1),
			ifelse(row+1 > (nrow+(rmin-1)), rmin, row+1), row)
	} else {
		y <- c(ifelse(col-1 < cmin, NA, col-1), col, col,
			ifelse(col+1 > (ncol+(cmin-1)), NA, col+1))
		x <- c(row, ifelse(row-1 < rmin, NA, row-1),
			ifelse(row+1 > (nrow+(rmin-1)), NA, row+1), row)
	}
	res <- as.data.frame(list(row=x, col=y))
	res <- na.omit(res)
	res <- as.matrix(res)
	rownames(res) <- NULL
	attr(res, "coords") <- c(col, row)
	res
}

queencell <- function(rowcol, nrow, ncol, torus=FALSE, rmin=1, cmin=1) {
	if (is.null(dim(rowcol))) rowcol <- t(as.matrix(rowcol))
        if(nrow(rowcol) != 1) stop("only single grid cell handled")
 	row <- rowcol[1]
	col <- rowcol[2]
	if (torus) {
		y <- c(rep(ifelse(col-1 < cmin, ncol, col-1), 3), col, col,
			rep(ifelse(col+1 > (ncol+(cmin-1)), cmin, col+1), 3))
		x <- integer(8)
		x[c(1,4,6)] <- rep(ifelse(row+1 > (nrow+(rmin-1)),
			rmin, row+1), 3)
		x[c(2,7)] <- rep(row, 2)
		x[c(3,5,8)] <- rep(ifelse(row-1 < rmin, nrow, row-1), 3)
	} else {
		y <- c(rep(ifelse(col-1 < cmin, NA, col-1), 3), col, col,
			rep(ifelse(col+1 > (ncol+(cmin-1)), NA, col+1), 3))
		x <- integer(8)
		x[c(1,4,6)] <- rep(ifelse(row+1 > (nrow+(rmin-1)),
			NA, row+1), 3)
		x[c(2,7)] <- rep(row, 2)
		x[c(3,5,8)] <- rep(ifelse(row-1 < rmin, NA, row-1), 3)
	}
	res <- as.data.frame(list(row=x, col=y))
	res <- na.omit(res)
	res <- as.matrix(res)
	rownames(res) <- NULL
	attr(res, "coords") <- c(col, row)
	res
}

mrc2vi <- function(rowcol, nrow, ncol) {
	i <- ((rowcol[,2]-1) * nrow) + rowcol[,1]
	if (i > nrow*ncol || i < 1) stop("row or column out of range")
	as.integer(i)
}

vi2mrc <- function(i, nrow, ncol) {
	col <- ceiling(i/nrow)
	tmp <- i %% nrow
	row <- ifelse(tmp == 0, nrow, tmp)
	if (row < 1 || row > nrow) stop("i out of range")
	if (col < 1 || col > ncol) stop("i out of range")
	res <- cbind(row, col)
	res
}

cell2nb <- function(nrow, ncol, type="rook", torus=FALSE) {
	nrow <- as.integer(nrow)
	if (nrow < 1) stop("nrow nonpositive")
	ncol <- as.integer(ncol)
	if (ncol < 1) stop("nrow nonpositive")
	xcell <- NULL
	if (type == "rook") xcell <- rookcell
	if (type == "queen") xcell <- queencell
	if (is.null(xcell))
		stop(paste(type, ": no such cell function", sep=""))
	n <- nrow * ncol
	res <- vector(mode="list", length=n)
	rownames <- character(n)
	for (i in 1:n) {
		res[[i]] <- sort(mrc2vi(xcell(vi2mrc(i, nrow, ncol),
			nrow, ncol, torus), nrow, ncol))
		rownames[i] <- paste(vi2mrc(i, nrow, ncol), collapse=":")
	}
	class(res) <- "nb"
	attr(res, "call") <- match.call()
	attr(res, "region.id") <- rownames
	attr(res, "cell") <- TRUE
	attr(res, type) <- TRUE
	if (torus) attr(res, "torus") <- TRUE
	res <- sym.attr.nb(res)
	invisible(res)
}


# Copyright 2004 by Roger Bivand 
#

choynowski <- function(n, x, row.names=NULL, tol=.Machine$double.eps^0.5) {
  len <- length(n)
  res <- numeric(len)
  nsum <- sum(n)
  xsum <- sum(x)
  b <- nsum/xsum
  E <- x*b
  type <- (n < E)
  for (i in 1:len) {
    if(type[i]) {
      for (j in 0:n[i]) {
        xx <- (E[i]^j*exp(-E[i])) / gamma(j + 1)
        res[i] <- res[i] + xx
      }
    } else {
      xx <- 1
      x <- n[i]
      while (xx > tol) {
        xx <- (E[i]^x*exp(-E[i])) / gamma(x + 1)
        res[i] <- res[i] + xx
        x <- x + 1
      }
    }
  }
  if (is.null(row.names)) 
    res <- data.frame(pmap=res, type=type)
  else
    res <- data.frame(pmap=res, type=type, row.names=row.names)
  res
}
# Copyright 2001 by Nicholas Lewin-Koh 
#


n.comp.nb <- function(nb.obj){
  if(!inherits(nb.obj,"nb"))stop("not a neighbours list")
  nb.obj <- make.sym.nb(nb.obj)
  comp <- rep(0,length(nb.obj))
  comp <- .Call("g_components", nb.obj, as.integer(comp), PACKAGE="spdep")
  answ <- list(nc=length(unique(comp)), comp.id=comp)
  answ
}

# Copyright 2001-3 by Roger Bivand 
#


diffnb <- function(x, y, verbose=TRUE) {
	if (!inherits(x, "nb")) stop("not a neighbours list")
	if (!inherits(y, "nb")) stop("not a neighbours list")
	n <- length(x)
	if(n != length(y)) stop("lengths differ")
	if (any(attr(x, "region.id") != attr(y, "region.id")))
		warning("region.id differ; using ids of first list")
	ids <- attr(x, "region.id")
	res <- vector(mode="list", length=n)
	for (i in 1:n) {
		xi <- x[[i]]
		yi <- y[[i]]
		xt <- xi %in% yi
		yt <- yi %in% xi
		if (!(all(xt) && all(yt))) {
			res[[i]] <- sort(unique(c(xi[which(!xt)],
				yi[which(!yt)])))
			if(verbose && (res[[i]] != 0))
				cat("Neighbour difference for region id:",
				ids[i], "in relation to id:", ids[res[[i]]], "\n")
		}
	}
	class(res) <- "nb"
	attr(res, "region.id") <- attr(x, "region.id")
	attr(res, "call") <- match.call()
	invisible(res)
}	
	
# Copyright 2000-3 by Roger S. Bivand. 
#

dnearneigh <- function(x, d1, d2, row.names=NULL, lonlat=FALSE) {
    if (!is.numeric(x)) stop("Data non-numeric")
    if (!is.matrix(x)) stop("Data not in matrix form")
    if (any(is.na(x))) stop("Data include NAs")
    if (!is.double(x)) storage.mode(x) <- "double"
    np <- nrow(x)
    if (!is.null(row.names)) {
	if(length(row.names) != np)
            stop("row.names wrong length")
	if (length(unique(row.names)) != length(row.names))
	    stop("non-unique row.names given")
    }
    if (is.null(row.names)) row.names <- as.character(1:np)
    dimension <- ncol(x)
    if (dimension > 2) stop("Only 2D data accepted")
    md <- 0
    if (d1 < 0) d1 <- 0.0
    if (!lonlat) {
	for (i in 1:dimension) md <- sum(md, (diff(range(x[,i]))^2))
    	if (d2 > sqrt(md)) d2 <- sqrt(md)
    }
    z <- .Call("dnearneigh", as.double(d1), as.double(d2), as.integer(np),
        as.integer(dimension), as.double(x), as.integer(lonlat), 
	PACKAGE="spdep")
    attr(z[[1]], "region.id") <- row.names
    attr(z[[1]], "call") <- match.call()
    attr(z[[1]], "dnn") <- c(d1, d2)
    z[[1]] <- sym.attr.nb(z[[1]])
    invisible(z[[1]])
}
# Copyright 2001 by Roger Bivand 
#

droplinks <- function(nb, drop, sym=TRUE) {
	if (class(nb) != "nb") stop("Not a neighbours list")
	n <- length(nb)
	if (is.logical(drop)) {
		if(length(drop) != n) stop("Argument lengths differ")
		idrop <- which(drop == TRUE)
	} else if(is.character(drop)) {
		row.names <- as.character(attr(nb, "region.id"))
		idrop <- match(drop, row.names)
		if(any(is.na(idrop))) stop("Region to drop not found")
	} else {
		idrop <- match(drop, 1:n)
		if(any(is.na(idrop))) stop("Region to drop not found")
	}
	if((attr(nb, "sym") == FALSE) && (sym == TRUE)) {
		warning("setting sym to FALSE")
		sym <- FALSE
	}
	for (i in idrop) {
		if (sym) {
			for (j in nb[[i]])
				nb[[j]] <- nb[[j]][nb[[j]] != i]
		}
		nb[[i]] <- 0
	}
	invisible(nb)
}

# Copyright 2001 by Roger Bivand and Nicholas Lewin-Koh
#

edit.nb <- function(name, coords, polys=NULL, ...) {
  nb <- name
  cnb <- card(nb)
  if (class(nb) != "nb") stop("Not a neighbours list")
  x <- coords[,1]
  y <- coords[,2]
  n <- length(nb)
  row.names <- attr(nb, "region.id")
  if (is.null(row.names)) row.names <- as.character(1:n)
  xlim <- range(x)
  ylim <- range(y)
  plot.new()
  plot.window(xlim = xlim, ylim = ylim, "", asp=1)
  if (!is.null(polys))
    plot(polys, border="grey", add=TRUE)
  for (i in 1:n) {
    #arrows(x[i],y[i],x[nb[[i]]],y[nb[[i]]],lenght=.08, angle=.15)
###
    if (cnb[i] > 0) segments(x[i],y[i],x[nb[[i]]],y[nb[[i]]])
###
    #inb <- nb[[i]]
    #for (j in inb)
    #lines(c(x[i], x[j]), c(y[i], y[j]), col="black")
  }
  points(x, y)
  finished <- "n"
  deletions <- NULL
  additions <- NULL
###
  edit.segs<-list()
  e.seg.stat<-NULL
  enum<-0
  erase.col<-par()$bg
###
  while (finished == "n") {
    cat("Identifying contiguity for deletion ...\n")
    cand <- identify(x, y, n=2)
    lines(x[cand], y[cand], col="red")

    if (.Platform$OS.type == "windows") bringToTop(-1)
    if ((cand[2] %in% nb[[cand[1]]]) && (cand[1] %in% nb[[cand[2]]])) {
      delete <- readline("Delete this line (y/n) ")
      if (delete != "y") delete <- "n"
      else {
        deletions <- c(deletions, paste(cand, collapse="-"))
        nb[[cand[1]]] <- nb[[cand[1]]][nb[[cand[1]]] != cand[2]]
				if(length(nb[[cand[1]]]) == 0) {
                                  nb[[cand[1]]] <- 0
                                  cat(cand[1], "is now an island\n")
				}
        nb[[cand[2]]] <- nb[[cand[2]]][nb[[cand[2]]] != cand[1]]
        if(length(nb[[cand[2]]]) == 0) {
          nb[[cand[2]]] <- 0
          cat(cand[2], "is now an island\n")
        }
###
        lines(x[cand], y[cand], col=erase.col)
        lines(x[cand], y[cand], col='brown',lty=4)
        enum<-enum+1
        edit.segs[[enum]]<-cand
        e.seg.stat<-c(e.seg.stat,0)
###
        cat("deleted contiguity between point", cand[1], "and", cand[2], "\n")
      }

      #plot.new()
      #plot.window(xlim = xlim, ylim = ylim, "", asp=1)
      #if (!is.null(polys))
      #  plot(polys, border="grey", add=TRUE)
      #for (i in 1:n) {
      #  inb <- nb[[i]]
      #  for (j in inb)
      #    lines(c(x[i], x[j]), c(y[i], y[j]),
      #          col="black")
      #}
      #points(x, y)
    }
      else {
        if (length(cand == 2)) {
          cat("No contiguity between chosen points\n")
          addcont <- readline("Add contiguity? (y/n) ")
          if (addcont != "y") addcont <- "n"
          if (addcont == "y") {
            nb[[cand[1]]] <-
              sort(unique(c(nb[[cand[1]]], cand[2])))
            nb[[cand[2]]] <-
              sort(unique(c(nb[[cand[2]]], cand[1])))
            cat("added contiguity between point",
                cand[1], "and", cand[2], "\n")
            additions <- c(additions, paste(cand, collapse="-"))
###
            enum<-enum+1
            edit.segs[[enum]]<-cand
            e.seg.stat<-c(e.seg.stat,1)
            lines(x[cand], y[cand], col='yellow')
###
          }
#          plot.new()
#          plot.window(xlim = xlim, ylim = ylim, "", asp=1)
#          if (!is.null(polys))
#            plot(polys, border="grey", add=TRUE)
#          for (i in 1:n) {
#            inb <- nb[[i]]
#            for (j in inb)
#              lines(c(x[i], x[j]),
#                    c(y[i], y[j]),
#                    col="black")
#				}
#          points(x, y)
        }
      }
#    finished <- readline("Finished yet? (y/n) ")
### 
    finished <- readline("Options: quit[q] refresh[r] continue[c] ")
    if (finished == "r") {
      plot.new()
      plot.window(xlim = xlim, ylim = ylim, "", asp=1)
      if (!is.null(polys))
        plot(polys, border="grey", add=TRUE)
      for (i in 1:n) {
        if(nb[[i]][1]!=0 & length(nb[[i]])>0)
          segments(x[i],y[i],x[nb[[i]]],y[nb[[i]]])
      }
      if(enum>1){
        for(i in 1:enum){
          if(e.seg.stat[i]==0){
            lines(x[edit.segs[[i]]], y[edit.segs[[i]]], col=erase.col)
            lines(x[edit.segs[[i]]], y[edit.segs[[i]]], col='brown',lty=4)
          }
          else lines(x[edit.segs[[i]]], y[edit.segs[[i]]], col='yellow')
        }
      }
      points(x, y)
      finished <- readline("Options: quit[q] continue[c]")
    }
    if (finished != "q") finished <- "n"
####
  }
  
  attributes(nb) <- list(deleted=deletions)
  attr(nb, "added") <- additions
  attr(nb, "region.id") <- row.names
  class(nb) <- "nb"
  nb <- sym.attr.nb(nb)
  invisible(nb)
}

# Copyright 2002-3 by Roger Bivand 
#

eigenw <- function(listw, quiet=TRUE)
{
	if(!inherits(listw, "listw")) stop("not a listw object")
	w <- listw2mat(listw)
	e <- eigen(w, only.values=TRUE)$values
	if (!quiet) {
		cat("Largest eigenvalue:", 
		if(is.complex(e)) max(Re(e)) else max(e),
		"Sum of eigenvalues:", sum(e), "\n")
	}
	e
}

# Copyright 1998-2004 by Roger Bivand (non-W styles Rein Halbersma)
#

errorsarlm <- function(formula, data = list(), listw, na.action=na.fail, 
	method="eigen", quiet=TRUE, zero.policy=FALSE, tol.solve=1.0e-10, 
        tol.opt=.Machine$double.eps^0.5, sparsedebug=FALSE) {
	mt <- terms(formula, data = data)
	mf <- lm(formula, data, na.action=na.action, method="model.frame")
	na.act <- attr(mf, "na.action")
	if (!inherits(listw, "listw")) stop("No neighbourhood list")
	if (!is.null(na.act)) {
	    subset <- !(1:length(listw$neighbours) %in% na.act)
	    listw <- subset(listw, subset, zero.policy=zero.policy)
	}
	if (!quiet) cat(paste("\nSpatial autoregressive error model\n", 
		"Jacobian calculated using "))
	switch(method,
		eigen = if (!quiet) cat("neighbourhood matrix eigenvalues\n"),
		sparse = if (!quiet) cat("sparse matrix techniques\n"),
		stop("...\n\nUnknown method\n"))
	y <- model.response(mf, "numeric")
	if (any(is.na(y))) stop("NAs in dependent variable")
	x <- model.matrix(mt, mf)
	if (any(is.na(x))) stop("NAs in independent variable")
	if (NROW(x) != length(listw$neighbours))
	    stop("Input data and neighbourhood list have different dimensions")
	n <- NROW(x)
	m <- NCOL(x)
	xcolnames <- colnames(x)
	K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1)
	wy <- lag.listw(listw, y, zero.policy=zero.policy)
	if (any(is.na(wy)))
	    stop("NAs in lagged dependent variable")
	if (m > 1) {
	    WX <- matrix(nrow=n,ncol=(m-(K-1)))
	    for (k in K:m) {
		wx <- lag.listw(listw, x[,k], zero.policy=zero.policy)
		if (any(is.na(wx)))
		    stop("NAs in lagged independent variable")
		WX[,(k-(K-1))] <- wx
	    }
	}
	if (K == 2) {
# modified to meet other styles, email from Rein Halbersma
		wx1 <- as.double(rep(1, n))
		wx <- lag.listw(listw, wx1, zero.policy=zero.policy)
		if (m > 1) WX <- cbind(wx, WX)
		else WX <- matrix(wx, nrow=n, ncol=1)
	}
	colnames(WX) <- xcolnames
	rm(wx)
	if (method == "eigen") {
		if (!quiet) cat("Computing eigenvalues ...\n")
		eig <- eigenw(listw)
		if (!quiet) cat("\n")
# range inverted 031031, email from Salvati Nicola (and Rein Halbersma)
		if (is.complex(eig)) eig.range <- 1/range(Re(eig))
		else eig.range <- 1/range(eig)
		opt <- optimize(sar.error.f, interval=eig.range, maximum=TRUE,
			tol=tol.opt, eig=eig,
			y=y, wy=wy, x=x, WX=WX, n=n, quiet=quiet)
	} else {
		sn <- listw2sn(listw)
		opt <- optimize(sar.error.f.s, interval=c(-1,1), maximum=TRUE,
			tol=tol.opt, sn=sn, y=y, wy=wy, x=x, WX=WX, 
			n=n, quiet=quiet, sparsedebug=sparsedebug)
	}
	lambda <- opt$maximum
	names(lambda) <- "lambda"
	LL <- opt$objective
	lm.target <- lm(I(y - lambda*wy) ~ I(x - lambda*WX) - 1)
	r <- as.vector(residuals(lm.target))
	fit <- as.vector(y - r)
	p <- lm.target$rank
	SSE <- deviance(lm.target)
	s2 <- SSE/n
	rest.se <- (summary(lm.target)$coefficients[,2])*sqrt((n-p)/n)
	coef.lambda <- coefficients(lm.target)
	names(coef.lambda) <- xcolnames
	lm.model <- lm(formula, data)
	ase <- FALSE
	lambda.se <- NULL
	LMtest <- NULL
	asyvar1 <- FALSE
	if (method == "eigen") {
		tr <- function(A) sum(diag(A))
		W <- listw2mat(listw)
		A <- solve(diag(n) - lambda*W, tol=tol.solve)
		WA <- W %*% A
		asyvar <- matrix(0, nrow=2+p, ncol=2+p)
		asyvar[1,1] <- n / (2*(s2^2))
		asyvar[2,1] <- asyvar[1,2] <- tr(WA) / s2
		asyvar[2,2] <- tr(WA %*% WA) + tr(t(WA) %*% WA)
		asyvar[3:(p+2),3:(p+2)] <- s2*(t(x - lambda*WX) %*% 
			(x - lambda*WX))
		asyvar1 <- solve(asyvar, tol=tol.solve)
		rownames(asyvar1) <- colnames(asyvar1) <- 
			c("sigma", "lambda", xcolnames)
		
		lambda.se <- sqrt(asyvar1[2,2])
		ase <- TRUE
	}
	call <- match.call()
	names(r) <- names(y)
	names(fit) <- names(y)
	ret <- structure(list(type="error", lambda=lambda, 
		coefficients=coef.lambda, rest.se=rest.se, 
		LL=LL, s2=s2, SSE=SSE, parameters=(m+2), lm.model=lm.model, 
		method=method, call=call, residuals=r, lm.target=lm.target,
		fitted.values=fit, ase=ase, formula=formula,
		se.fit=NULL, resvar=asyvar1,
		lambda.se=lambda.se, LMtest=LMtest, zero.policy=zero.policy), 
		class=c("sarlm"))
	if (zero.policy) {
		zero.regs <- attr(listw$neighbours, 
			"region.id")[which(card(listw$neighbours) == 0)]
		if (length(zero.regs) > 0)
			attr(ret, "zero.regs") <- zero.regs
	}
	if (!is.null(na.act))
		ret$na.action <- na.act
	ret
}

sar.error.f <- function(lambda, eig, y, wy, x, WX, n, quiet)
{
	yl <- y - lambda*wy
	xl <- x - lambda*WX
	xl.q <- qr.Q(qr(xl))
	xl.q.yl <- t(xl.q) %*% yl
	SSE <- t(yl) %*% yl - t(xl.q.yl) %*% xl.q.yl
	s2 <- SSE/n
	if (is.complex(eig)) det <- Re(prod(1 - lambda*eig)) 
	else det <- prod(1 - lambda*eig)
	ret <- (log(det) - ((n/2)*log(2*pi)) - (n/2)*log(s2) - (1/(2*(s2)))*SSE)
	if (!quiet) cat("Lambda:\t", lambda, "\tfunction value:\t", ret, "\n")
	ret
}

sar.error.f.s <- function(lambda, sn, y, wy, x, WX, n, quiet, sparsedebug)
{
	yl <- y - lambda*wy
	xl <- x - lambda*WX
	xl.q <- qr.Q(qr(xl))
	xl.q.yl <- t(xl.q) %*% yl
	SSE <- t(yl) %*% yl - t(xl.q.yl) %*% xl.q.yl
	s2 <- SSE/n
	ret <- (logSpwdet(sparseweights=sn, rho=lambda, debug=sparsedebug) - 
		((n/2)*log(2*pi)) - (n/2)*log(s2) - (1/(2*(s2)))*SSE)
	if (!quiet) cat("Lambda:\t", lambda, "\tfunction value:\t", ret, "\n")
	ret
}

# Copyright 2001 by Nicholas Lewin-Koh 
#

gabrielneigh <- function(coords) {
    x <- coords
    if (!is.matrix(x)) stop("Data not in matrix form")
    if (any(is.na(x))) stop("Data cannot include NAs")
    np <- nrow(x)
    if(ncol(x)!=2) stop("Planar graphs only work in 2d")
    g1<-g2<-rep(0,np*3)
    nogab <- 0
    z <- .C("compute_gabriel", np=as.integer(np), from=as.integer(g1),
             to=as.integer(g2), nedges=as.integer(nogab),
             x=as.double(x[,1]), y=as.double(x[,2]), PACKAGE="spdep")
    z$from<-z$from[1:z$nedges]
    z$to<-z$to[1:z$nedges]
    attr(z, "call") <- match.call()
    class(z)<-c("Graph","Gabriel")
    invisible(z)
}

plot.Gabriel<-function(x, show.points=FALSE, add=FALSE,
                       linecol=par(col), ...)
{
  if(!add) plot(gab$x,gab$y,type='n',...)
  segments(gab$x[gab$from], gab$y[gab$from],
           gab$x[gab$to], gab$y[gab$to], col=linecol)
  if(show.points) points(gab$x,gab$y)
}

# Copyright 2001-3 by Roger Bivand 
#


geary <- function(x, listw, n, n1, S0, zero.policy=FALSE) {
	z <- scale(x, scale=FALSE)
	zz <- sum(z^2)
	K <- (n*sum(z^4))/(zz^2)
	res <- geary.intern(x, listw, n, zero.policy, type="geary")
	C <- (n1 / (2*S0)) * (sum(res) / zz)
	res <- list(C=C, K=K)
	res
}

geary.intern <- function(x, listw, n, zero.policy, type="geary") {
	cardnb <- card(listw$neighbours)
	if (type == "geary") ft <- TRUE
	else if (type == "sokal") ft <- FALSE
	else stop("type unknown")
	res <- .Call("gearyw", listw$neighbours, listw$weights,
		as.numeric(x), as.integer(cardnb),
		as.logical(zero.policy), as.logical(ft), PACKAGE="spdep")
	if (any(is.na(res))) warning("NAs in lagged values")
	invisible(res)
}

geary.test <- function(x, listw, randomisation=TRUE, zero.policy=FALSE,
    alternative="less", spChk=NULL) {
	alternative <- match.arg(alternative, c("less", "greater", "two.sided"))
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.numeric(x)) stop(paste(deparse(substitute(x)),
		"is not a numeric vector"))
	if (any(is.na(x))) stop("NA in X")
	n <- length(listw$neighbours)
	if (n != length(x)) stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
	
	wc <- spweights.constants(listw, zero.policy)
	S02 <- wc$S0*wc$S0
	res <- geary(x, listw, wc$n, wc$n1, wc$S0, zero.policy)
	C <- res$C
	if (is.na(C)) stop("NAs generated in geary - check zero.policy")
	K <- res$K
	EC <- 1
	if(randomisation) {
		VC <- (wc$n1*wc$S1*(wc$nn - 3*n + 3 - K*wc$n1))
		VC <- VC - ((1/4) * (wc$n1*wc$S2*(wc$nn + 3*n - 6 - 
			K*(wc$nn - n + 2))))
		VC <- VC + (S02*(wc$nn - 3 - K*(wc$n1^2)))
		VC <- VC / (n*wc$n2*wc$n3*S02)
	} else {
		VC <- ((2*wc$S1 + wc$S2)*wc$n1 - 4*S02) / (2*(n + 1)*S02)
	}
	ZC <- (C - EC) / sqrt(VC)
	statistic <- ZC
	names(statistic) <- "Geary C statistic standard deviate"
	PrC <- NA
	if (is.finite(ZC)) {
        	if (alternative == "two.sided") PrC <- 2 * pnorm(abs(ZC), 
			lower.tail=FALSE)
        	else if (alternative == "greater")
            	PrC <- pnorm(ZC, lower.tail=FALSE)
        	else PrC <- pnorm(ZC)
		if (PrC < 0 || PrC > 1) 
		    warning("Out-of-range p-value: reconsider test arguments")
	}
	vec <- c(C, EC, VC)
	names(vec) <- c("Geary C statistic", "Expectation", "Variance")
	method <- paste("Geary's C test under", ifelse(randomisation,
	    "randomisation", "normality"))
	data.name <- paste(deparse(substitute(x)), "\nweights:",
	    deparse(substitute(listw)), "\n")
	res <- list(statistic=statistic, p.value=PrC, estimate=vec, 
	    alternative=alternative, method=method, data.name=data.name)
	class(res) <- "htest"
	res
}

geary.mc <- function(x, listw, nsim, zero.policy=FALSE,
	alternative="less", spChk=NULL) {
	alternative <- match.arg(alternative, c("less", "greater"))
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.numeric(x)) stop(paste(deparse(substitute(x)),
		"is not a numeric vector"))
	if(missing(nsim)) stop("nsim must be given")
	if (any(is.na(x))) stop("NA in X")
	n <- length(listw$neighbours)
	if (n != length(x)) stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
	if(nsim > gamma(n+1)) stop("nsim too large for this n")
	wc <- spweights.constants(listw, zero.policy)
	res <- numeric(length=nsim+1)
	for (i in 1:nsim) res[i] <- geary(sample(x), listw, n, wc$n1, wc$S0,
	    zero.policy)$C
	res[nsim+1] <- geary(x, listw, n, wc$n1, wc$S0, zero.policy)$C
	rankres <- rank(res)
	xrank <- rankres[length(res)]
	diff <- nsim - xrank
	diff <- ifelse(diff > 0, diff, 0)
        if (alternative == "less") 
        	pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
    	else if (alternative == "greater") 
        	pval <- punif((diff + 1)/(nsim + 1))
	if (pval < 0 || pval > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
	statistic <- res[nsim+1]
	names(statistic) <- "statistic"
	parameter <- xrank
	names(parameter) <- "observed rank"
	method <- "Monte-Carlo simulation of Geary's C"
	data.name <- paste(deparse(substitute(x)), "\nweights:",
	    deparse(substitute(listw)), "\nnumber of simulations + 1:",
	    nsim+1, "\n")
	lres <- list(statistic=statistic, parameter=parameter,
	    p.value=pval, alternative=alternative, method=method, 
	    data.name=data.name, res=res)
	class(lres) <- c("htest", "mc.sim")
	lres
}

# Copyright 2002-3 by Hisaji ONO and Roger Bivand 
#
# General G Statistics
#
#
globalG.test <- function(x, listw, zero.policy=FALSE,
	alternative="greater", spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
	if (!inherits(listw, "listw"))
	stop(paste(deparse(substitute(listw)), "is not a listw object"))
	if (listw$style != "B") stop("Only binary weights allowed")
	if (!is.numeric(x))
	stop(paste(deparse(substitute(x)), "is not a numeric vector"))
	if (any(is.na(x))) stop(paste("NA in ", deparse(substitute(x))))
	if (any(x < 0.0)) 
		stop(paste("Negative value in ", deparse(substitute(x))))
	n <- length(listw$neighbours)
	if (n != length(x))stop("Different numbers of observations")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")

	wc <- spweights.constants(listw, zero.policy=zero.policy)
	n1 <- n - 1
	n2 <- n - 2
	n3 <- n - 3
	nn <- n * n
	S0 <- wc$S0
	S1 <- wc$S1
	S2 <- wc$S2
	S02 <- S0*S0
	G <- (t(x) %*% lag.listw(listw, x, zero.policy=zero.policy)) /
		(sum(x %x% x) - (t(x) %*% x))

	E.G <- S0 / (n * n1)

	B0 <- ((nn - 3*n + 3)*S1) - (n*S2) + (3*S02)
	B1 <- -(((nn - n)*S1) - (2*n*S2) + (6*S02))
	B2 <- -((2*n*S1) - ((n+3)*S2) + (6*S02))
	B3 <- (4*n1*S1) - (2*(n+1)*S2) + (8*S02)
	B4 <- S1 - S2 + S02
	sx <- sum(x)
	sx2 <- sum(x^2)
	sx3 <- sum(x^3)
	sx4 <- sum(x^4)

	var.G <- ((B0*(sx2^2) + B1*sx4 + B2*(sx^2)*sx2 + B3*sx*sx3 +
		 B4*(sx^4)) / ((((sx^2) - sx2)^2)*n*n1*n2*n3)) - (E.G^2)

	statistic <- (G - E.G) / sqrt(var.G)
	names(statistic) <- "Getis-Ord global G statistic"
	if (alternative == "two.sided") PrG <- 2 * pnorm(-abs(statistic), 
		lower.tail=FALSE)
        else if (alternative == "greater")
            PrG <- pnorm(statistic, lower.tail=FALSE)
        else PrG <- pnorm(statistic)
	if (PrG < 0 || PrG > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
	vec <- c(G, E.G, var.G)
	names(vec) <- c("Global G statistic", "Expectation", "Variance")
	data.name <- paste(deparse(substitute(x)), "\nweights:",
	    deparse(substitute(listw)), "\n")
	res <- list(statistic=statistic, p.value=PrG, estimate=vec, 
	    alternative=alternative, data.name=data.name)
	class(res) <- "htest"
	res
}

# Copyright 2001-3 by Nicholas Lewin-Koh and Roger S. Bivand.
#


graph2nb <- function(gob, row.names=NULL,sym=FALSE) {
	if (!inherits(gob, "Graph")) stop("Not a Graph object")
	res <- vector(mode="list", length=gob$np)
    	if (!is.null(row.names)) {
		if(length(row.names) != gob$np)
            		stop("row.names wrong length")
		if (length(unique(row.names)) != length(row.names))
	    		stop("non-unique row.names given")
    	}
    	if (is.null(row.names)) row.names <- as.character(1:gob$np)
        if(sym){
          for (i in 1:gob$np) {
		res[[i]] <- sort(unique(c(gob$to[gob$from==i],
                                       gob$from[gob$to==i])))
	  	if(length(res[[i]]) == 0) res[[i]] <- 0
	  }
        }
        else{
	  for (i in 1:gob$np) {
		res[[i]] <- sort(gob$to[gob$from==i])
	  	if(length(res[[i]]) == 0) res[[i]] <- 0
	  }
        }
        attr(res, "region.id") <- row.names
 	attr(res, "call") <- attr(gob, "call")
 	attr(res, "type") <- attr(gob, "type")
	class(res) <- "nb"
	res <- sym.attr.nb(res)
	invisible(res)
}
# Copyright 2001-4 by Roger Bivand 
#

joincount <- function(dums, listw) {
	nc <- which(colSums(dums) > 1)
	n <- length(listw$neighbours)
	cardnb <- card(listw$neighbours)
	res <- as.numeric(rep(0, ncol(dums)))
	for (lev in nc) {
		res[lev] <- .Call("jcintern", listw$neighbours,
			listw$weights, as.integer(dums[,lev]),
			as.integer(cardnb), PACKAGE="spdep")
	}
	res
}

joincount.test <- function(fx, listw, zero.policy=FALSE,
	alternative="greater", spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if (!is.factor(fx)) stop(paste(deparse(substitute(x)),
		"is not a factor"))
	if (any(is.na(fx))) stop("NA in factor")
	n <- length(listw$neighbours)
	if (n != length(fx)) stop("objects of different length")
	cards <- card(listw$neighbours)
	if (!zero.policy && any(cards == 0))
		stop("regions with no neighbours found")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(fx, listw))
		stop("Check of data and weights ID integrity failed")

	wc <- spweights.constants(listw, zero.policy=zero.policy)
	S02 <- wc$S0*wc$S0

	ff <- ~ fx - 1
	dums <- model.matrix(ff, model.frame(ff))
	BB <- joincount(dums, listw)
	nBB <- length(BB)
	res <- vector(mode="list", length=nBB)
	tab <- table(fx)
	BB5 <- 0.5 * BB
	ntab <- as.vector(tab)
	Ejc <- (wc$S0*(ntab*(ntab-1))) / (2*n*wc$n1)
	Vjc <- (wc$S1*(ntab*(ntab-1))) / (n*wc$n1)
	Vjc <- Vjc + (((wc$S2 - 2*wc$S1)*ntab*(ntab-1)*(ntab-2)) /
		(n*wc$n1*wc$n2))
	Vjc <- Vjc + (((S02 + wc$S1 - wc$S2)*ntab*(ntab-1)*(ntab-2)*
		(ntab-3)) / (n*wc$n1*wc$n2*wc$n3))
	Vjc <- (0.25 * Vjc) - Ejc^2
	for (i in 1:nBB) {
		estimate <- c(BB5[i], Ejc[i], Vjc[i])
		names(estimate) <- c("Same colour statistic",
			"Expectation", "Variance")
		statistic <- (BB5[i] - Ejc[i]) / sqrt(Vjc[i])
		names(statistic) <- paste("Std. deviate for", names(tab)[i])
		p.value <- NA
		if (is.finite(statistic)) {
		    if (alternative == "two.sided") 
			p.value <- 2 * pnorm(abs(statistic), lower.tail=FALSE)
		    else if (alternative == "greater")
			p.value <- pnorm(statistic, lower.tail=FALSE)
		    else p.value <- pnorm(statistic)
		    if (p.value < 0 || p.value > 1) 
		      warning("Out-of-range p-value: reconsider test arguments")
		}
		method <- "Join count test under nonfree sampling"
		data.name <- paste(deparse(substitute(fx)), "\nweights:",
			deparse(substitute(listw)), "\n")
		res[[i]] <- list(statistic=statistic, p.value=p.value,
			estimate=estimate, method=method,
			alternative=alternative, data.name=data.name)
		class(res[[i]]) <- "htest"
	}
	class(res) <- "jclist"
	res
}

print.jclist <- function(x, ...) {
	for (i in 1:length(x)) print(x[[i]], ...)
	invisible(x)
}

joincount.mc <- function(fx, listw, nsim, zero.policy=FALSE,
	alternative="greater", spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less"))
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.factor(fx)) stop(paste(deparse(substitute(fx)),
		"is not a factor"))
	if(missing(nsim)) stop("nsim must be given")
	if (any(is.na(fx))) stop("NA in factor")
	n <- length(listw$neighbours)
	if (n != length(fx)) stop("objects of different length")
	cards <- card(listw$neighbours)
	if (!zero.policy && any(cards == 0))
		stop("regions with no neighbours found")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(fx, listw))
		stop("Check of data and weights ID integrity failed")
	if(nsim > gamma(n+1)) stop("nsim too large for this n")

	ff <- ~ fx - 1
	dums <- model.matrix(ff, model.frame(ff))
	nc <- ncol(dums)
	res <- matrix(0, nrow=nsim+1, ncol=nc)
	res[nsim+1,] <- 0.5 * joincount(dums, listw)
	tab <- table(fx)
	for (i in 1:nsim) {
		fxi <- sample(fx)
		ff <- ~ fxi - 1
		dums <- model.matrix(ff, model.frame(ff))
		res[i,] <- 0.5 * joincount(dums, listw)
	}
	rankres <- apply(res, 2, rank)
	xrank <- rankres[nrow(rankres),]
	lres <- vector(mode="list", length=nc)
	for (i in 1:nc) {
		statistic <- res[nrow(res), i]
		names(statistic) <- paste("Join-count statistic for",
			names(tab)[i])
		parameter <- xrank[i]
		names(parameter) <- "rank of observed statistic"
		diff <- nsim - xrank[i]
		diff <- ifelse(diff > 0, diff, 0)
        	if (alternative == "less") 
        		pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
    		else if (alternative == "greater") 
        		pval <- punif((diff + 1)/(nsim + 1))
		if (pval < 0 || pval > 1) 
		    warning("Out-of-range p-value: reconsider test arguments")

		method <- "Monte-Carlo simulation of join-count statistic"
		data.name <- paste(deparse(substitute(fx)), "\nweights:",
			deparse(substitute(listw)),
			"\nnumber of simulations + 1:", nsim+1, "\n")
		estimate <- c(mean(res[-(nrow(res)), i]),
			var(res[-(nrow(res)), i]))
		names(estimate) <- c("mean of simulation",
			"variance of simulation")
		lres[[i]] <- list(statistic=statistic, parameter=parameter,
			method=method, data.name=data.name, p.value=pval, 
			alternative=alternative, estimate=estimate, res=res[,i])
		class(lres[[i]]) <- c("htest", "mc.sim")
		
	}
	class(lres) <- "jclist"
	lres
}



joincount.multi <- function(fx, listw, zero.policy=FALSE,
	spChk=NULL) {
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.factor(fx)) stop(paste(deparse(substitute(fx)),
		"is not a factor"))
	if (any(is.na(fx))) stop("NA in factor")
	n <- length(listw$neighbours)
	if (n != length(fx)) stop("objects of different length")
	cards <- card(listw$neighbours)
	if (!zero.policy && any(cards == 0))
		stop("regions with no neighbours found")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(fx, listw))
		stop("Check of data and weights ID integrity failed")
	ifx <- as.integer(fx)
	k <- length(levels(fx))
	sn <- listw2sn(listw)
	y <- factor(paste(ifx[sn[,1]], ":", ifx[sn[,2]], sep=""), 
		levels=as.vector(outer(1:k, 1:k, 
			FUN=function(X,Y) paste(X,Y,sep=":"))))
	res <- matrix(tapply(sn[,3], y, sum), ncol=k)/2
		
	res[is.na(res)] <- 0
	rownames(res) <- colnames(res) <- levels(fx)

	tab <- table(fx)
	ntab <- as.vector(tab)
	wc <- spweights.constants(listw, zero.policy=zero.policy)
	S02 <- wc$S0*wc$S0

	Ejc <- (wc$S0*(ntab*(ntab-1))) / (2*n*wc$n1)

	Vjc <- (wc$S1*(ntab*(ntab-1))) / (n*wc$n1)
	Vjc <- Vjc + (((wc$S2 - 2*wc$S1)*ntab*(ntab-1)*(ntab-2)) /
		(n*wc$n1*wc$n2))
	Vjc <- Vjc + (((S02 + wc$S1 - wc$S2)*ntab*(ntab-1)*(ntab-2)*
		(ntab-3)) / (n*wc$n1*wc$n2*wc$n3))
	Vjc <- (0.25 * Vjc) - Ejc^2

	nrns <- function(x, op="*") {
		k <- length(x)
		res <- numeric(((k^2) - k)/2)
		ii <- 1
		for (i in 2:k) {
			for (j in 1:(i-1)) {
				if (is.character(op) && op == "*") {
					res[ii] <- x[i]*x[j]
				} else if (is.character(op) && op == "+") {
					res[ii] <- x[i]+x[j]
				}
				ii <- ii+1
			}
		}
		res
	}

	ldiag <- numeric(((k^2) - k)/2)
	diffcolnames <- character(((k^2) - k)/2)
	ii <- 1
	for (i in 2:k) {
		for (j in 1:(i-1)) {
			ldiag[ii] <- res[i,j] + res[j,i]
			diffcolnames[ii] <- paste(levels(fx)[i], 
				levels(fx)[j], sep=":")
			ii <- ii+1
		}
	}

	Exp <- (wc$S0*(nrns(ntab, op="*"))) / (n*wc$n1)
	Var <- (2*wc$S1*nrns(ntab, op="*"))/(n*wc$n1)
	Var <- Var + (((wc$S2 - 2*wc$S1)*nrns(ntab, op="*")*
		(nrns(ntab, op="+")-2))/(n*wc$n1*wc$n2))
	Var <- Var + ((4*(S02 + wc$S1 - wc$S2)*nrns((ntab*(ntab-1)), op="*")) /
		(n*wc$n1*wc$n2*wc$n3))
	Var <- (0.25 * Var) - Exp^2
	Jtot <- sum(ldiag)
	JtotExp <- sum(Exp)
	Jvar <- ((wc$S2/(n*wc$n1))-((4*(S02 + wc$S1 - wc$S2)*wc$n1) /
		(n*wc$n1*wc$n2*wc$n3)))*sum(nrns(ntab, op="*"))
	Jvar <- Jvar + 4*(((wc$S1 - wc$S2)/(n*wc$n1*wc$n2*wc$n3)) + 
		((2*S02*(2*n-3))/((n*wc$n1)*(n*wc$n1*wc$n2*wc$n3))))*
		sum(nrns(ntab^2, op="*"))
	if(k>2) {
		ntnsnr <- 0
		for (r in 1:(k-2)) {
			for (s in (r+1):(k-1)) {
				for (t in (s+1):(k)) {
					ntnsnr <- ntnsnr +
						ntab[r]*ntab[s]*ntab[t]
				}
			}
		}
		Jvar <- Jvar + (((2*wc$S1 - 5*wc$S2)/(n*wc$n1*wc$n2))+
		((12*(S02 + wc$S1 - wc$S2))/(n*wc$n1*wc$n2*wc$n3))+
		((8*S02)/((n*wc$n1*wc$n2)*wc$n1)))*ntnsnr
	}
	if(k>3) {
		nuntnsnr <- 0
		for (r in 1:(k-3)) {
			for (s in (r+1):(k-2)) {
				for (t in (s+1):(k-1)) {
					for (u in (t+1):(k)) {
						nuntnsnr <- nuntnsnr +
						ntab[r]*ntab[s]*ntab[t]*ntab[u]
					}
				}
			}
		}
		Jvar <- Jvar - 8*(((wc$S1 - wc$S2)/(n*wc$n1*wc$n2*wc$n3))+
		((2*S02*(2*n-3))/((n*wc$n1)*(n*wc$n1*wc$n2*wc$n3))))*nuntnsnr
	}
	Jvar <- (0.25 * Jvar)
	statistic <- (c(diag(res), ldiag, Jtot) - c(Ejc, Exp, JtotExp)) / 
		sqrt(c(Vjc, Var, Jvar))
	lres <- cbind(c(diag(res), ldiag, Jtot), c(Ejc, Exp, JtotExp), 
		c(Vjc, Var, Jvar), statistic)
	colnames(lres) <- c("Joincount", "Expected", "Variance", 
		"z-value")
	rownames(lres) <- c(paste(levels(fx), ":", levels(fx), sep=""), 
		diffcolnames, "Jtot")
	class(lres) <- c("jcmulti", "matrix")
	lres
}

print.jcmulti <- function(x, ...) {
	printCoefmat(x, ...)
}



#joincount.diffcol <- function(fx, listw, zero.policy=FALSE,
#	spChk=NULL) {
#	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
#		"is not a listw object"))
#	if(!is.factor(fx)) stop(paste(deparse(substitute(fx)),
#		"is not a factor"))
#	if (any(is.na(fx))) stop("NA in factor")
#	n <- length(listw$neighbours)
#	if (n != length(fx)) stop("objects of different length")
#	cards <- card(listw$neighbours)
#	if (!zero.policy && any(cards == 0))
#		stop("regions with no neighbours found")
#	if (is.null(spChk)) spChk <- get.spChkOption()
#	if (spChk && !chkIDs(fx, listw))
#		stop("Check of data and weights ID integrity failed")
#	ifx <- as.integer(fx)
#	
#	cardnb <- card(listw$neighbours)
#	res <- .Call("jcdiffs", listw$neighbours,
#			listw$weights, as.integer(ifx),
#			as.integer(cardnb), PACKAGE="spdep")
#	
#	res <- list(diffcount=res, tot=sum(cardnb)/2)
#	res
#}

# Copyright 2001-3 by Roger S. Bivand. 
#

knearneigh <- function(x, k=1, lonlat=FALSE)
{
    if (!is.numeric(x)) stop("Data non-numeric")
    if (!is.matrix(x)) stop("Data not in matrix form")
    if (any(is.na(x))) stop("Data include NAs")
    if (!is.double(x)) x <- as.double(x)
    np <- nrow(x)
    dimension <- ncol(x)
    if (dimension > 2) stop("Only 2D data accepted")
    if (k >= np) stop("Fewer data points than k")
    xx <- x[,1]
    if (dimension > 1) for (i in 2:dimension) xx <- c(xx, x[,i])
    nn <- integer(np*k)
    dnn <- double(np*k)
    z <- .C("knearneigh", k=as.integer(k), np=as.integer(np),
        dimension=as.integer(dimension),
        xx=as.double(xx), nn=as.integer(nn), dnn=as.double(dnn),
	as.integer(lonlat), PACKAGE="spdep")
    res <- list(nn=matrix(z$nn, np, k, byrow=TRUE), np=np, k=k,
    	dimension=dimension, x=x)
    class(res) <- "knn"
    attr(res, "call") <- match.call()
    invisible(res)
}
# Copyright 2001 by Roger Bivand
#


knn2nb <- function(knn, row.names=NULL, sym=FALSE) {
	if (class(knn) != "knn") stop("Not a knn object")
	res <- vector(mode="list", length=knn$np)
    	if (!is.null(row.names)) {
		if(length(row.names) != knn$np)
            		stop("row.names wrong length")
		if (length(unique(row.names)) != length(row.names))
	    		stop("non-unique row.names given")
    	}
    	if (is.null(row.names)) row.names <- as.character(1:knn$np)
        if(sym){
          to<-as.vector(knn$nn)
          from<-rep(1:knn$np,knn$k)
          for (i in 1:knn$np)res[[i]] <- sort(unique(c(to[from==i],
                                                       from[to==i]))) 
        } else {
          for (i in 1:knn$np) res[[i]] <- sort(knn$nn[i,])
        }
 	attr(res, "region.id") <- row.names
 	attr(res, "call") <- attr(knn, "call")
        attr(res, "sym") <- sym
	attr(res, "type") <- "knn"
 	attr(res, "knn-k") <- knn$k
	class(res) <- "nb"
	invisible(res)
}
# Copyright 1998-2004 by Roger Bivand and Andrew Bernat
#

lagsarlm <- function(formula, data = list(), listw, 
	na.action=na.fail, type="lag", method="eigen", quiet=TRUE, 
	zero.policy=FALSE, tol.solve=1.0e-7, tol.opt=.Machine$double.eps^0.5, 
	sparsedebug=FALSE) {
	mt <- terms(formula, data = data)
	mf <- lm(formula, data, na.action=na.action, 
		method="model.frame")
	na.act <- attr(mf, "na.action")
	if (!inherits(listw, "listw")) stop("No neighbourhood list")
	if (!is.null(na.act)) {
	    subset <- !(1:length(listw$neighbours) %in% na.act)
	    listw <- subset(listw, subset, zero.policy=zero.policy)
	}
	switch(type, lag = if (!quiet) cat("\nSpatial lag model\n"),
	    mixed = if (!quiet) cat("\nSpatial mixed autoregressive model\n"),
	    stop("\nUnknown model type\n"))
	if (!quiet) cat("Jacobian calculated using ")
	switch(method, 
		eigen = if (!quiet) cat("neighbourhood matrix eigenvalues\n"),
		sparse = if (!quiet) cat("sparse matrix techniques\n"),
		stop("...\nUnknown method\n"))
	y <- model.response(mf, "numeric")
#	if (any(is.na(y))) stop("NAs in dependent variable")
	x <- model.matrix(mt, mf)
#	if (any(is.na(x))) stop("NAs in independent variable")
	if (NROW(x) != length(listw$neighbours))
		stop("Input data and weights have different dimensions")
	n <- NROW(x)
	m <- NCOL(x)
	xcolnames <- colnames(x)
	K <- ifelse(xcolnames[1] == "(Intercept)", 2, 1)
	wy <- lag.listw(listw, y, zero.policy=zero.policy)
	if (any(is.na(wy))) stop("NAs in lagged dependent variable")
	if (type != "lag") {
		# check if there are enough regressors
	        if (m > 1) {
			WX <- matrix(nrow=n,ncol=(m-(K-1)))
			for (k in K:m) {
				wx <- lag.listw(listw, x[,k], 
				    zero.policy=zero.policy)
				if (any(is.na(wx))) 
				    stop("NAs in lagged independent variable")
				WX[,(k-(K-1))] <- wx
			}
		}
		if (K == 2) {
         	    # unnormalized weight matrices
                	if (!(listw$style == "W")) {
 	      			intercept <- as.double(rep(1, n))
       	        		wx <- lag.listw(listw, intercept, 
					zero.policy = zero.policy)
                    		if (m > 1) {
                        		WX <- cbind(wx, WX)
                    		} else {
			      		WX <- matrix(wx, nrow = n, ncol = 1)
                    		}
                	} 
            	}   
		m1 <- m + 1
		mm <- NCOL(x) + NCOL(WX)
            	xxcolnames <- character(mm)
		for (k in 1:m) xxcolnames[k] <- xcolnames[k]
		for (k in m1:mm) 
		    xxcolnames[k] <- paste("lag.", xcolnames[k-mm+m], sep="")
		x <- cbind(x, WX)
		colnames(x) <- xxcolnames
		m <- NCOL(x)
		rm(wx, WX)
	}
	if (method == "eigen") {
		if (!quiet) cat("Computing eigenvalues ...\n")
		eig <- eigenw(listw)
		if (!quiet) cat("\n")
#range inverted 031031, email from Salvati Nicola (and Rein Halbersma)
		if (is.complex(eig)) eig.range <- 1/range(Re(eig))
		else eig.range <- 1/range(eig)
		lm.null <- lm(y ~ x - 1)
		lm.w <- lm.fit(x, wy)
		e.null <- lm.null$residuals
		e.w <- lm.w$residuals
		e.a <- t(e.null) %*% e.null
		e.b <- t(e.w) %*% e.null
		e.c <- t(e.w) %*% e.w
		opt <- optimize(sar.lag.mixed.f, interval=eig.range,
			maximum=TRUE, tol=tol.opt, eig=eig,
			e.a=e.a, e.b=e.b, e.c=e.c, n=n, quiet=quiet)
	} else {
		opt <- dosparse(listw, y, x, wy, K, quiet, tol.opt, 
			sparsedebug)
	}
	rho <- c(opt$maximum)
	names(rho) <- "rho"
	LL <- c(opt$objective)
	lm.lag <- lm((y - rho*wy) ~ x - 1)
	r <- residuals(lm.lag)
	fit <- y - r
	names(r) <- names(fit)
	coef.rho <- coefficients(lm.lag)
	names(coef.rho) <- colnames(x)
	SSE <- deviance(lm.lag)
	s2 <- SSE/n
	if (method != "eigen") {
		LLs <- opt$LLs
		lm.null <- opt$lm.null
		rest.se <- NULL
		rho.se <- NULL
		LMtest <- NULL
		ase <- FALSE
		varb <- FALSE
	} else {
		LLs <- NULL
		tr <- function(A) sum(diag(A))
# beware of complex eigenvalues!
		O <- (eig/(1-rho*eig))^2
		omega <- sum(O)
		if (is.complex(omega)) omega <- Re(omega)
		W <- listw2mat(listw)
		A <- solve(diag(n) - rho*W, tol=tol.solve)
		AW <- A %*% W
		zero <- rbind(rep(0,length(coef.rho)))
		xtawxb <- s2*(t(x) %*% AW %*% x %*% coef.rho)
		V <- s2*(s2*tr(t(AW) %*% AW) +
			t(AW %*% x %*% coef.rho) %*%
			(AW %*% x %*% coef.rho)) + omega*s2^2
		inf1 <- rbind(n/2, s2*tr(AW), t(zero))
		inf2 <- rbind(s2*tr(AW), V, xtawxb)
		xtx <- s2*t(x) %*% x
		inf3 <- rbind(zero, t(xtawxb), xtx)
		inf <- cbind(inf1, inf2, inf3)
		varb <- (s2^2) * solve(inf, tol=tol.solve)
		rownames(varb) <- colnames(varb) <- 
			c("sigma", "rho", colnames(x))
		rest.se <- sqrt(diag(varb))[-c(1:2)]
		rho.se <- sqrt(varb[2,2])
		TW <- (W %*% W) + (t(W) %*% W)
		T22 <- sum(diag(TW))
		T21A <- sum(diag(TW %*% A))
		LMtest <- ((t(r) %*% W %*% r)/s2)^2
		LMtest <- LMtest/(T22 - ((T21A^2)*(rho.se^2)))
		ase <- TRUE
	}
	call <- match.call()
	ret <- structure(list(type=type, rho=rho, 
		coefficients=coef.rho, rest.se=rest.se, 
		LL=LL, s2=s2, SSE=SSE, parameters=(m+2), lm.model=lm.null,
		method=method, call=call, residuals=r, 
		lm.target=lm.lag, fitted.values=fit,
		se.fit=NULL, formula=formula,
		ase=ase, LLs=LLs, rho.se=rho.se, LMtest=LMtest, 
		resvar=varb, zero.policy=zero.policy), class=c("sarlm"))
	if (zero.policy) {
		zero.regs <- attr(listw$neighbours, 
			"region.id")[which(card(listw$neighbours) == 0)]
		if (length(zero.regs) > 0)
			attr(ret, "zero.regs") <- zero.regs
	}
	if (!is.null(na.act))
		ret$na.action <- na.act
	ret
}

sar.lag.mixed.f <- function(rho, eig, e.a, e.b, e.c, n, quiet)
{
	SSE <- e.a - 2*rho*e.b + rho*rho*e.c
	s2 <- SSE/n
	if (is.complex(eig)) det <- Re(prod(1 - rho*eig)) 
	else det <- prod(1 - rho*eig)
	ret <- (log(det) - ((n/2)*log(2*pi)) - (n/2)*log(s2)
		- (1/(2*s2))*SSE)
	if (!quiet) cat("Rho:\t", rho, "\tfunction value:\t", ret, "\n")
	ret
}

sar.lag.mixed.f.s <- function(rho, sn, e.a, e.b, e.c, n, quiet, sparsedebug)
{
	SSE <- e.a - 2*rho*e.b + rho*rho*e.c
	s2 <- SSE/n
	ret <- (logSpwdet(sparseweights=sn, rho=rho, debug=sparsedebug)
		- ((n/2)*log(2*pi)) - (n/2)*log(s2) - (1/(2*s2))*SSE)
	if (!quiet) cat("Rho:\t", rho, "\tfunction value:\t", ret, "\n")
	ret
}


dosparse <- function (listw, y, x, wy, K, quiet, tol.opt, sparsedebug) {
	sn <- listw2sn(listw)
	m <- ncol(x)
	n <- nrow(x)
	LLs <- vector(mode="list", length=length(K:m))
	j <- 1
	for (i in K:m) {
		thisx <- x[,-i]
		lm.null <- lm.fit(thisx, y)
		lm.w <- lm.fit(thisx, wy)
		e.null <- lm.null$residuals
		e.w <- lm.w$residuals
		e.a <- t(e.null) %*% e.null
		e.b <- t(e.w) %*% e.null
		e.c <- t(e.w) %*% e.w
		LLs[[j]] <- optimize(sar.lag.mixed.f.s, interval=c(-1,1),
		maximum=TRUE, tol=tol.opt, sn=sn,
		e.a=e.a, e.b=e.b, e.c=e.c, n=n, quiet=quiet,
		sparsedebug=sparsedebug)$objective
		attr(LLs[[j]], "nall") <- n
		attr(LLs[[j]], "nobs") <- n
		attr(LLs[[j]], "df") <- (m+2)-1
		attr(LLs[[j]], "name") <- colnames(x)[i]
		class(LLs[[j]]) <- "logLik"
		j <- j + 1
	}
	lm.null <- lm(y ~ x - 1)
	lm.w <- lm.fit(x, wy)
	e.null <- lm.null$residuals
	e.w <- lm.w$residuals
	e.a <- t(e.null) %*% e.null
	e.b <- t(e.w) %*% e.null
	e.c <- t(e.w) %*% e.w
	sn <- listw2sn(listw)
	opt <- optimize(sar.lag.mixed.f.s, interval=c(-1,1),
		maximum=TRUE, tol=tol.opt, sn=sn,
		e.a=e.a, e.b=e.b, e.c=e.c, n=n, quiet=quiet, 
		sparsedebug=sparsedebug)
	maximum <- opt$maximum
	objective <- opt$objective
	res <- list(maximum=maximum, objective=objective, LLs=LLs,
		lm.null=lm.null)
}
# Copyright 2001-3 by Roger Bivand
#

listw2sn <- function(listw) {
	if(!inherits(listw, "listw")) stop("not a listw object")
	cardw <- card(listw$neighbours)
	scard <- sum(cardw)
	z <- .Call("listw2sn", listw$neighbours, listw$weights,
		as.integer(cardw), as.integer(scard), PACKAGE="spdep")
	res <- as.data.frame(list(from=z[[1]], to=z[[2]], weights=z[[3]]))
	class(res) <- c(class(res), "spatial.neighbour")
	attr(res, "n") <- length(attr(listw, "region.id"))
	attr(res, "region.id") <- attr(listw, "region.id")
	neighbours.attrs <- names(attributes(listw$neighbours))
	attr(res, "neighbours.attrs") <- neighbours.attrs
	weights.attrs <- names(attributes(listw$weights))
	attr(res, "weights.attrs") <- weights.attrs
	if (!(is.null(attr(listw, "GeoDa"))))
		attr(res, "GeoDa") <- attr(listw, "GeoDa")
	attr(res, "listw.call") <- attr(listw, "call")
	invisible(res)
}

sn2listw <- function(sn) {
	if(!inherits(sn, "spatial.neighbour")) 
	    stop("not a spatial.neighbour object")
	n <- attr(sn, "n")
	region.id <- attr(sn, "region.id")
	nlist <- vector(mode="list", length=n)
	class(nlist) <- "nb"
	attr(nlist, "region.id") <- region.id
	vlist <- vector(mode="list", length=n)
	rle.sn <- rle(sn[,1])
	if (n != length(rle.sn$lengths))
	    warning(paste(paste(region.id[which(!(1:n %in% rle.sn$values))], 
		collapse=", "), "are not origins"))
	cs1.sn <- cumsum(rle.sn$lengths)
	cs0.sn <- c(1, cs1.sn[1:(n-1)]+1)
	ii <- 1
	for (i in 1:n) {
		if (rle.sn$value[ii] == i) {
			nlist[[i]] <- as.integer(sn[cs0.sn[ii]:cs1.sn[ii],2])
			vlist[[i]] <- as.double(sn[cs0.sn[ii]:cs1.sn[ii],3])
			ii <- ii+1
		} else {
			nlist[[i]] <- as.integer(0)
		}
	}
	res <- list(style=as.character(NA), neighbours=nlist, weights=vlist)
	class(res) <- c("listw", "nb")
	if (!(is.null(attr(sn, "GeoDa"))))
		attr(res, "GeoDa") <- attr(sn, "GeoDa")
	attr(res, "region.id") <- region.id
	attr(res, "call") <- match.call()
	invisible(res)
}

# Copyright 2001-3 by Roger Bivand 
#

lm.LMtests <- function(model, listw, zero.policy=FALSE, test="LMerr",
	spChk=NULL) {
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(class(model) != "lm") stop(paste(deparse(substitute(model)),
		"not an lm object"))
	N <- length(listw$neighbours)
	u <- resid(model)
	if (N != length(u)) stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(u, listw))
		stop("Check of data and weights ID integrity failed")
	u <- as.vector(u)

	if (is.null(attr(listw$weights, "W")) || !attr(listw$weights, "W"))
		warning("Spatial weights matrix not row standardized")
	all.tests <- c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA")
	if (test[1] == "all") test <- all.tests
	if (!all(test %in% all.tests))
		stop("Invalid test selected - must be either \"all\" or a vector of tests")		

	y <- model.response(model.frame(model))
	X <- model.matrix(terms(model), model.frame(model))
	yhat <- as.vector(fitted(model))
	p <- model$rank
	p1 <- 1:p
	XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
	sigma2 <- (t(u) %*% u) / N
	TrW <- tracew(listw)
	Wu <- lag.listw(listw, u, zero.policy)
	Wy <- lag.listw(listw, y, zero.policy)
	Wyhat <- lag.listw(listw, yhat, zero.policy)
	XtWyhat <- t(X) %*% Wyhat
	dutWu <- (t(u) %*% Wu) / sigma2
	resa <- (dutWu ^ 2) / TrW
	J <- (1/(N*sigma2)) *
		((t(Wyhat) %*% Wyhat) -
		(t(XtWyhat) %*% XtXinv %*% XtWyhat) +
		(TrW * sigma2))
	dutWy <- (t(u) %*% Wy) / sigma2
	nt <- length(test)
	tres <- vector(mode="list", length=nt)
	names(tres) <- test
	for (i in 1:nt) {
		testi <- test[i]
		zz <- switch(testi,
		LMerr = vec <- c(resa, 1),
		LMlag = vec <- c((dutWy ^ 2) / (N * J), 1),
		RLMerr = vec <- c(((dutWu - (TrW*((N*J)^-1))*dutWy)^2) /
			(TrW * (1 - TrW*((N*J)^-1))), 1),
		RLMlag = vec <- c(((dutWy - dutWu)^2)/ ((N*J) - TrW), 1),
		SARMA = vec <- c(((dutWy - dutWu)^2)/ ((N*J) - TrW) + resa, 2)
		)
		if (is.null(zz)) stop(paste(testi, ": no such test", sep=""))
		statistic <- vec[1]
		names(statistic) <- testi
		parameter <- vec[2]
		names(parameter) <- "df"
		p.value <- 1 - pchisq(statistic, parameter)
		if (p.value < 0 || p.value > 1) 
		    warning("Out-of-range p-value: reconsider test arguments")
		names(p.value) <- ""
		method <- "Lagrange multiplier diagnostics for spatial dependence"
		data.name <- paste("\n", paste(strwrap(paste("model: ",
		    gsub(" *", " ", 
	            paste(deparse(model$call), sep="", collapse="")))),
		    collapse="\n"),
    	            "\nweights: ", deparse(substitute(listw)), "\n", sep="")
		tres[[i]] <- list(statistic=statistic, parameter=parameter,
			p.value=p.value, method=method, data.name=data.name)
		class(tres[[i]]) <- "htest"
	}
	class(tres) <- "LMtestlist"
	tres
}

print.LMtestlist <- function(x, ...) {
	for (i in 1:length(x)) print(x[[i]])
	invisible(x)
}

tracew <- function (listw) {
	dlmtr <- 0
	n <- length(listw$neighbours)
	for (i in 1:n) {
		dij <- listw$neighbours[[i]]
		ndij <- length(dij)
		wdij <- listw$weights[[i]]
		for (j in 1:ndij) {
			k <- dij[j]
			if (k > i) {
			    dk <- which(listw$neighbours[[k]] == i)
			    if (length(dk) > 0 && dk > 0 &&
				dk <= length(listw$neighbours[[k]]))
				wdk <- listw$weights[[k]][dk]
				else wdk <- 0
				dlmtr <- dlmtr + (wdk * wdk) + 2 *
				(wdij[j] * wdk) + (wdij[j] * wdij[j])
			}
		}
	}
	dlmtr
}

# Copyright 2001-3 by Roger Bivand 
#

lm.morantest <- function(model, listw, zero.policy=FALSE, 
	    alternative = "greater", spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(class(model) != "lm") stop(paste(deparse(substitute(model)),
		"not an lm object"))
 	N <- length(listw$neighbours)
	u <- residuals(model)
	if (N != length(u)) 
            stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(u, listw))
		stop("Check of data and weights ID integrity failed")
	
	u <- as.vector(u)
	listw.U <- listw2U(listw)

	S0 <- sum(unlist(listw.U$weights))
	S1 <- 0.5 * sum((2*unlist(listw.U$weights))^2)
	lu <- lag.listw(listw.U, u, zero.policy=zero.policy)
	if (zero.policy) N <- length(which(card(listw$neighbours) > 0))
	I <- (N/S0) * ((t(u) %*% lu) / (t(u) %*% u))
	p <- model$rank
	p1 <- 1:p
	XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
	X <- model.matrix(terms(model), model.frame(model))
# Cliff/Ord 1981, p. 203
	Z <- lag.listw(listw.U, X, zero.policy=zero.policy)
	C1 <- t(X) %*% Z
	trA <- (sum(diag(XtXinv %*% C1)))
	EI <- -((N * trA) / ((N-p) * S0))
# minus changed from trA to EI (Luis Galvis, Dec 2, 2003)
	C2 <- t(Z) %*% Z
	C3 <- XtXinv %*% C1
	trA2 <- sum(diag(C3 %*% C3))
	trB <- sum(diag(4*(XtXinv %*% C2)))
	VI <- (((N*N)/((S0*S0)*(N-p)*(N-p+2))) *
		(S1 + 2*trA2 - trB - ((2*(trA^2))/(N-p))))
	ZI <- (I - EI) / sqrt(VI)
    	if (alternative == "two.sided") pv <- 2 * pnorm(abs(ZI), 
		lower.tail=FALSE)
    	else if (alternative == "greater")
	        pv <- pnorm(ZI, lower.tail=FALSE)
    	else pv <- pnorm(ZI)
	if (pv < 0 || pv > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
    	statistic <- ZI
    	attr(statistic, "names") <- "Moran I statistic standard deviate"
    	p.value <- pv
    	estimate <- c(I, EI, VI)
    	attr(estimate, "names") <- c("Observed Moran's I", "Expectation",
	    "Variance")
    	method <- "Global Moran's I for regression residuals"
    	data.name <- paste("\n", paste(strwrap(paste("model: ",
	    gsub(" *", " ", 
	    paste(deparse(model$call), sep="", collapse="")))), collapse="\n"),
    	    "\nweights: ", deparse(substitute(listw)), "\n", sep="")
    	res <- list(statistic = statistic, p.value = p.value,
	       estimate = estimate, method = method,
		alternative = alternative, data.name = data.name)
	class(res) <- "htest"
    	res
}
# Copyright 2001-3 by Roger Bivand 
#

localG <- function(x, listw, zero.policy=FALSE, spChk=NULL) {
	if (!inherits(listw, "listw"))
		stop(paste(deparse(substitute(listw)), "is not a listw object"))
	if (!is.numeric(x))
		stop(paste(deparse(substitute(x)), "is not a numeric vector"))
	if (any(is.na(x))) stop(paste("NA in ", deparse(substitute(x))))
	n <- length(listw$neighbours)
	if (n != length(x))stop("Different numbers of observations")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
	gstari <- FALSE
	if (!is.null(attr(listw$neighbours, "self.included")) &&
		attr(listw$neighbours, "self.included")) gstari <- TRUE
	lx <- lag.listw(listw, x, zero.policy=zero.policy)
	if (gstari) {
		xibar <- rep(mean(x), n)
		si2 <- rep(sum(scale(x, scale=FALSE)^2)/n, n)
	} else {
		xibar <- (rep(sum(x),n) - x) / (n - 1)
		si2 <- ((rep(sum(x^2),n) - x^2) / (n - 1)) - xibar^2
	}
	Wi <- sapply(listw$weights, sum)
	S1i <- sapply(listw$weights, function(x) sum(x^2))
	res <- (lx - Wi*xibar)
	if (gstari) {
		res <- res / sqrt(si2*((n*S1i - Wi^2)/(n-1)))
	} else {
		res <- res / sqrt(si2*(((n-1)*S1i - Wi^2)/(n-2)))
	}
	attr(res, "gstari") <- gstari
	attr(res, "call") <- match.call()
	class(res) <- "localG"
	invisible(res)
}
# Copyright 2001-4 by Roger Bivand 
#

localmoran <- function(x, listw, zero.policy=FALSE, na.action=na.fail, 
	alternative = "greater", p.adjust.method="none", spChk=NULL) {
	if (!inherits(listw, "listw"))
		stop(paste(deparse(substitute(listw)), "is not a listw object"))
	if (!is.null(attr(listw$neighbours, "self.included")) &&
		attr(listw$neighbours, "self.included"))
		stop("Self included among neighbours")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
	if (!is.numeric(x))
		stop(paste(deparse(substitute(x)), "is not a numeric vector"))
#	if (any(is.na(x))) stop(paste("NA in ", deparse(substitute(x))))
	NAOK <- deparse(substitute(na.action)) == "na.pass"
	x <- na.action(x)
	na.act <- attr(x, "na.action")
	if (!is.null(na.act)) {
	    subset <- !(1:length(listw$neighbours) %in% na.act)
	    listw <- subset(listw, subset, zero.policy=zero.policy)
	    excl <- class(na.act) == "exclude"
	}
	n <- length(listw$neighbours)
	if (n != length(x))stop("Different numbers of observations")
	res <- matrix(nrow=n, ncol=5)
        if (alternative == "two.sided") Prname <- "Pr(z != 0)"
        else if (alternative == "greater") Prname <- "Pr(z > 0)"
        else Prname <- "Pr(z < 0)"
	colnames(res) <- c("Ii", "E.Ii", "Var.Ii", "Z.Ii", Prname)
	xx <- mean(x, na.rm=NAOK)
	z <- x - xx
#	z <- scale(x, scale=FALSE)
	lz <- lag.listw(listw, z, zero.policy=zero.policy, NAOK=NAOK)
	s2 <- sum(z^2, na.rm=NAOK)/n
	res[,1] <- (z/s2) * lz
	Wi <- sapply(listw$weights, sum)
	res[,2] <- -Wi / (n-1)
	b2 <- (sum(z^4, na.rm=NAOK)/n)/(s2^2)
	A <- (n-b2) / (n-1)
	B <- (2*b2 - n) / ((n-1)*(n-2))
	C <- Wi^2 / ((n-1)^2)
	Wi2 <- sapply(listw$weights, function(x) sum(x^2))
	Wikh2 <- sapply(listw$weights, function(x) {
		ifelse(is.null(x), 0, 1 - crossprod(x,x))
	})
	res[,3] <- A*Wi2 + B*Wikh2 - C
	res[,4] <- (res[,1] - res[,2]) / sqrt(res[,3])
        if (alternative == "two.sided") pv <- 2 * pnorm(abs(res[,4]), 
	    lower.tail=FALSE)
        else if (alternative == "greater")
            pv <- pnorm(res[,4], lower.tail=FALSE)
        else pv <- pnorm(res[,4])
	res[,5] <- p.adjustSP(pv, listw$neighbours, method=p.adjust.method)
	if (!is.null(na.act) && excl) {
		res <- naresid(na.act, res)
	}
	attr(res, "call") <- match.call()
	if (!is.null(na.act)) attr(res, "na.action") <- na.act
	class(res) <- c("localmoran", class(res))
	invisible(res)
}


# Copyright 2001-4 by Roger Bivand 
#

moran <- function(x, listw, n, S0, zero.policy=FALSE, NAOK=FALSE) {
	n1 <- length(listw$neighbours)
	if (n1 != length(x)) stop("objects of different length")
	xx <- mean(x, na.rm=NAOK)
	z <- x - xx
	zz <- sum(z^2, na.rm=NAOK)
	K <- (length(x)*sum(z^4, na.rm=NAOK))/(zz^2)
	lz <- lag.listw(listw, z, zero.policy=zero.policy, NAOK=NAOK)
#	I <- (n / S0) * ((t(z) %*% lz) / zz)
	I <- (n / S0) * ((sum(z*lz, na.rm=NAOK)) / zz)
	res <- list(I=I, K=K)
	res
}

moran.test <- function(x, listw, randomisation=TRUE, zero.policy=FALSE,
	alternative="greater", rank = FALSE, na.action=na.fail, spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less", "two.sided"))
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if (!is.numeric(x)) stop(paste(deparse(substitute(x)),
		"is not a numeric vector"))
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
#	if (any(is.na(x))) stop("NA in X")
	xname <- deparse(substitute(x))
	wname <- deparse(substitute(listw))
	NAOK <- deparse(substitute(na.action)) == "na.pass"
	x <- na.action(x)
	na.act <- attr(x, "na.action")
	if (!is.null(na.act)) {
	    subset <- !(1:length(listw$neighbours) %in% na.act)
	    listw <- subset(listw, subset, zero.policy=zero.policy)
	}
	n <- length(listw$neighbours)
	if (n != length(x)) stop("objects of different length")
	
	wc <- spweights.constants(listw, zero.policy=zero.policy)
	S02 <- wc$S0*wc$S0
	res <- moran(x, listw, wc$n, wc$S0, zero.policy=zero.policy, 
		NAOK=NAOK)
	I <- res$I
	K <- res$K
	if (rank) K <- (3*(3*wc$n^2 -7))/(5*(wc$n^2 - 1))
	EI <- (-1) / wc$n1
	if(randomisation) {
		VI <- wc$n*(wc$S1*(wc$nn - 3*wc$n + 3) - wc$n*wc$S2 + 3*S02)
		tmp <- K*(wc$S1*(wc$nn - wc$n) - 2*wc$n*wc$S2 + 6*S02)
		VI <- (VI - tmp) / (wc$n1*wc$n2*wc$n3*S02)
		VI <- VI - EI^2
	} else {
		VI <- (wc$nn*wc$S1 - wc$n*wc$S2 + 3*S02) / (S02*(wc$nn - 1))
		VI <- VI - EI^2
	}
	ZI <- (I - EI) / sqrt(VI)
	statistic <- ZI
	names(statistic) <- "Moran I statistic standard deviate"
        if (alternative == "two.sided") 
		PrI <- 2 * pnorm(abs(ZI), lower.tail=FALSE)
        else if (alternative == "greater")
            PrI <- pnorm(ZI, lower.tail=FALSE)
        else PrI <- pnorm(ZI)
	if (PrI < 0 || PrI > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
	vec <- c(I, EI, VI)
	names(vec) <- c("Moran I statistic", "Expectation", "Variance")
	method <- paste("Moran's I test under", ifelse(randomisation,
	    "randomisation", "normality"))
	data.name <- paste(xname, ifelse(rank,
		"using rank correction",""), "\nweights:",
		wname, ifelse(is.null(na.act), "", paste("\nomitted:", 
	    paste(na.act, collapse=", "))), "\n")
	res <- list(statistic=statistic, p.value=PrI, estimate=vec, 
	    alternative=alternative, method=method, data.name=data.name)
	if (!is.null(na.act)) attr(res, "na.action") <- na.act
	class(res) <- "htest"
	res
}

moran.mc <- function(x, listw, nsim, zero.policy=FALSE,
	alternative="greater", na.action=na.fail, spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less"))
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.numeric(x)) stop(paste(deparse(substitute(x)),
		"is not a numeric vector"))
	if(missing(nsim)) stop("nsim must be given")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
#	if (any(is.na(x))) stop("NA in X")
	xname <- deparse(substitute(x))
	wname <- deparse(substitute(listw))
	if (deparse(substitute(na.action)) == "na.pass")
	    stop("na.pass not permitted")
	x <- na.action(x)
	na.act <- attr(x, "na.action")
	if (!is.null(na.act)) {
	    subset <- !(1:length(listw$neighbours) %in% na.act)
	    listw <- subset(listw, subset, zero.policy=zero.policy)
	}
	n <- length(listw$neighbours)
	if (n != length(x)) stop("objects of different length")
	if(nsim > gamma(n+1)) stop("nsim too large for this n")
	
	S0 <- Szero(listw)
	res <- numeric(length=nsim+1)
	for (i in 1:nsim) res[i] <- moran(sample(x), listw, n, S0,
	    zero.policy)$I
	res[nsim+1] <- moran(x, listw, n, S0, zero.policy)$I
	rankres <- rank(res)
	xrank <- rankres[length(res)]
	diff <- nsim - xrank
	diff <- ifelse(diff > 0, diff, 0)
	if (alternative == "less") 
        	pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
    	else if (alternative == "greater") 
        	pval <- punif((diff + 1)/(nsim + 1))
	if (pval < 0 || pval > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
	statistic <- res[nsim+1]
	names(statistic) <- "statistic"
	parameter <- xrank
	names(parameter) <- "observed rank"
	method <- "Monte-Carlo simulation of Moran's I"
	data.name <- paste(xname, "\nweights:",
	    wname, ifelse(is.null(na.act), "", paste("\nomitted:", 
	    paste(na.act, collapse=", "))), "\nnumber of simulations + 1:",
	    nsim+1, "\n")
	lres <- list(statistic=statistic, parameter=parameter,
	    p.value=pval, alternative=alternative, method=method, 
	    data.name=data.name, res=res)
	if (!is.null(na.act)) attr(lres, "na.action") <- na.act
	class(lres) <- c("htest", "mc.sim")
	lres
}


# Copyright 2001 by Roger Bivand 
#

moran.plot <- function(x, listw, zero.policy=FALSE, spChk=NULL,
 labels=NULL, xlab=NULL, ylab=NULL, quiet=FALSE, ...)
{
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	xname <- deparse(substitute(x))
	if (!is.numeric(x)) stop(paste(xname, "is not a numeric vector"))
	if (any(is.na(x))) stop("NA in X")
	n <- length(listw$neighbours)
	if (n != length(x)) stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(x, listw))
		stop("Check of data and weights ID integrity failed")
	labs <- TRUE
	if (is.logical(labels) && !labels) labs <- FALSE
	if (is.null(labels) || length(labels) != n)
		labels <- as.character(attr(listw, "region.id"))
	wx <- lag.listw(listw, x, zero.policy=zero.policy)
	if (is.null(xlab)) xlab <- xname
	if (is.null(ylab)) ylab <- paste("spatially lagged", xname)
	plot(x, wx, xlab=xlab, ylab=ylab, ...)
	if (zero.policy) {
		n0 <- wx == 0.0
		symbols(x[n0], wx[n0], inches=FALSE, 
		circles=rep(diff(range(x))/50, length(which(n0))),
		bg="grey", add=TRUE)
	}
	xwx.lm <- lm(wx ~ x)
	abline(xwx.lm)
	abline(h=mean(wx), lty=2)
	abline(v=mean(x), lty=2)
	infl.xwx <- influence.measures(xwx.lm)
	is.inf <- which(apply(infl.xwx$is.inf, 1, any))
	points(x[is.inf], wx[is.inf], pch=9, cex=1.2)
	if (labs)
	    text(x[is.inf], wx[is.inf], labels=labels[is.inf], pos=2, cex=0.7)
	rownames(infl.xwx$infmat) <- labels
	if (!quiet) summary(infl.xwx)
	invisible(infl.xwx)
}

# Copyright 2002-4 by Roger Bivand and Michael Tiefelsdorf,
# with contributions by Danlin Yu
#

localmoran.sad <- function (model, select, nb, glist = NULL, style = "W",
    zero.policy = FALSE, alternative = "greater", spChk=NULL,
    save.Vi = FALSE, tol = .Machine$double.eps^0.5,
    maxiter = 1000, tol.bounds=0.0001) {
# need to impose check on weights TODO!!
    if (class(nb) != "nb") 
        stop(paste(deparse(substitute(nb)), "not an nb object"))
#    if (class(model) != "lm") 
#        stop(paste(deparse(substitute(model)), "not an lm object"))
    clobj <- class(model)
    cond.sad <- FALSE
    n <- length(nb)
    if (clobj == "sarlm") {
    	errorsarLambda <- c(model$lambda)
	if (is.null(errorsarLambda)) stop(paste(deparse(substitute(model)), 
	    "not an error sarlm object"))
	model <- model$lm.model
	W <- nb2mat(nb, glist=glist, style=style, zero.policy=zero.policy)
	sqrtOmega <- solve(diag(n) - errorsarLambda * W)
	cond.sad <- TRUE
     } else if (clobj != "lm")
     	stop(paste(deparse(substitute(model)), "not an lm object"))
    u <- residuals(model)
    if (n != length(u)) 
        stop("objects of different length")
    if (is.null(spChk)) spChk <- get.spChkOption()
    if (spChk && !chkIDs(u, nb2listw(nb, zero.policy=zero.policy)))
	stop("Check of data and weights ID integrity failed")
    if (!(alternative %in% c("greater", "less", "two.sided")))
	stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
    u <- as.vector(u)
    select <- unique(as.integer(select))
    if (any(select < 1 || select > n))
        stop("select out of range")
    utu <- c(t(u) %*% u)
    p <- model$rank
    p1 <- 1:p
    m <- n - p - 2
    XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
    X <- model.matrix(terms(model), model.frame(model))
    B <- listw2U(nb2listw(nb, glist=glist, style="B",
	zero.policy=zero.policy))
    D <- NULL
    a <- NULL
    if (style == "W") {
        D <- 1/sapply(B$weights, sum)
    } else if (style == "S") {
        D <- 1 / sqrt(sapply(B$weights, function(x) sum(x^2)))
#        a <- sum(unlist(B$weights))
# correction by Danlin Yu, 25 March 2004
	a <- sum(sapply(B$weights, function(x) sqrt(sum(x^2))))
    } else if (style == "C") a <- sum(unlist(B$weights))
    res <- vector(mode="list", length=length(select))
    for (i in 1:length(select)) {
        Vi <- listw2star(B, select[i], style=style, n, D, a,
	    zero.policy=zero.policy)
        Viu <- lag.listw(Vi, u, zero.policy=TRUE)
	Ii <- c((t(u) %*% Viu) / utu)
	if (cond.sad) {
	    M <- diag(n) - X %*% XtXinv %*% t(X)
    	    ViI <- listw2mat(Vi) - Ii * diag(n)
	    innerTerm <- t(sqrtOmega) %*% M %*% ViI %*% M %*% sqrtOmega
	    evalue <- eigen(innerTerm, only.values=TRUE)$values
#	    is.na(evalue) <- abs(evalue) < 1.0e-7
#	    evalue <- na.omit(evalue)
#	    saveNaAction <- attr(evalue, "na.action")
	    saveNaAction <- NA
	    tau <- c(evalue)
	    e1 <- tau[1]
	    en <- tau[length(tau)]
#	    taumi <- tau
#    	    taumi <- tau - Ii
    	    low <- (1 / (2*tau[length(tau)])) + tol.bounds #+ 0.01
    	    high <- (1 / (2*tau[1])) - tol.bounds #- 0.01
    	    f <- function(omega, tau) {sum(tau/(1 - (2*omega*tau)))}
    	    root <- uniroot(f, lower=low, upper=high, tol=tol, maxiter=maxiter,
        	tau=tau)
    	    omega <- root$root
# 0 should be expectation - maybe use try()
	    if (omega < 0 ) sad.r <- try(-sqrt(sum(log(1 - 2*omega*tau))))
    	    else sad.r <- try(sqrt(sum(log(1 - 2*omega*tau))))
	    if (inherits(sad.r, "try.error")) {
	    	warning (paste("In zone:", select[i], "sad.r not a number"))
	    	sad.r <- sad.u <- sad.p <- NaN
	    } else { 
		sad.u <- omega * sqrt(2*sum(tau^2 / (1 - (2*omega*tau))^2))
    	    	sad.p <- sad.r - ((1/sad.r)*log(sad.r/sad.u))
	    }
	} else {
	    ViX <- lag.listw(Vi, X, zero.policy=TRUE)
	    MViM <- t(X) %*% ViX %*% XtXinv
	    t1 <- -sum(diag(MViM))
	    sumsq.Vi <- function(x) {
            	if (is.null(x)) NA
	    	else sum(x^2)
	    }
	    trVi2 <- sum(sapply(Vi$weights, sumsq.Vi), na.rm=TRUE)
	    t2a <- sum(diag(t(ViX) %*% ViX %*% XtXinv))
	    t2b <- sum(diag(MViM %*% MViM))
	    t2 <- trVi2 - 2*t2a + t2b
	    e1 <- 0.5 * (t1 + sqrt(2*t2 - t1^2))
	    en <- 0.5 * (t1 - sqrt(2*t2 - t1^2))
	    l <- en
	    h <- e1
	    mi <- Ii
	    aroot= m*mi*(l+h-2*mi)+mi*(3*l+3*h-4*mi)-2*l*h
            broot= (m+2)*mi*(l-mi)*(h-mi)
            c1root= l**2 * mi**2 * (m+1)**2 + h**2 * mi**2 * (m+1)**2
            c2root= 2*l*h * (2*l*h - 2*l*mi - 2*h*mi - 2*m*mi**2 -
	    	m**2 * mi**2 + mi**2)
            omega= 0.25*((aroot-sqrt(c1root+c2root))/broot)
	    if (is.nan(omega)) {
	    	warning (paste("In zone:", select[i], "omega not a number"))
	    	sad.r <- sad.u <- sad.p <- NaN
	    } else { 
            	tau <- c(c(e1), rep(0, m), c(en))
	    	taumi <- tau - Ii
            	if (omega < 0 ) sad.r <- -sqrt(sum(log(1 - 2*omega*taumi)))
            	else sad.r <- sqrt(sum(log(1 - 2*omega*taumi)))
            	sad.u <- omega * sqrt(2*sum(taumi^2 / (1 - (2*omega*taumi))^2))
            	sad.p <- sad.r - ((1/sad.r)*log(sad.r/sad.u))
	    }
	}
        if (alternative == "two.sided") p.sad <- 2 * pnorm(abs(sad.p), 
	    lower.tail=FALSE)
        else if (alternative == "greater")
            p.sad <- pnorm(sad.p, lower.tail=FALSE)
        else p.sad <- pnorm(sad.p)
        statistic <- sad.p
        attr(statistic, "names") <- "Saddlepoint approximation"
        p.value <- p.sad
        estimate <- c(Ii)
        attr(estimate, "names") <- "Observed Moran's Ii"
        internal1 <- c(omega, sad.r, sad.u)
        attr(internal1, "names") <- c("omega", "sad.r", "sad.u")
        method <- paste("Saddlepoint approximation for local Moran's I",
            "(Barndorff-Nielsen formula)")
        data.name <- paste("region:", select[i],
	    attr(nb, "region.id")[select[i]],
	    "\n", paste(strwrap(paste("model: ", gsub(" *", " ", 
	    paste(deparse(model$call), sep="", collapse="")))),
	    collapse="\n"),
            "\nneighbours:", deparse(substitute(nb)),
	    "style:", style, "\n")
        obj <- list(statistic = statistic, p.value = p.value,
            estimate = estimate, method = method,
	    alternative = alternative, data.name = data.name,
	    internal1 = internal1, df = (n-p), tau = c(c(e1), c(en)),
	    i = paste(select[i], attr(nb, "region.id")[select[i]]),
#	    if (save.Vi) {Vi = Vi}
	    Vi = if(save.Vi) Vi else NULL)
        class(obj) <- "moransad"
	res[[i]] <- obj
    }
    class(res) <- "localmoransad"
    invisible(res)
}

print.localmoransad <- function(x, ...) {
    extract <- function(x, i) {x[[i]]}
    regnames <- sapply(x, extract, 10)
    est <- sapply(x, extract, 3)
    sad <- sapply(x, extract, 1)
    pval <- sapply(x, extract, 2)
    res <- as.matrix(cbind(est, sad, pval))
    rownames(res) <- regnames
    colnames(res) <- c("Local Morans I", "Saddlepoint", "Pr. (Sad)")
    print(res, ...)
    invisible(res)
}
as.data.frame.localmoransad <- function(x, row.names=NULL, optional=FALSE) {
    n <- length(x)
    res <- matrix(0, nrow=n, ncol=14)
    regnames <- NULL
    if (!is.null(row.names)) 
	if (length(row.names) == n) regnames <- row.names
    if (is.null(regnames))for (i in 1:n) regnames <- c(regnames, x[[i]]$i)
    for (i in 1:n) {
        tau <- x[[i]]$tau
	df <- x[[i]]$df
        tau <- c(tau[1], rep(0, df-2), tau[2])
        max.I <- tau[1]
        min.I <- tau[length(tau)]
        E.I <- sum(tau)/df
        tau <- tau - E.I
        V.I <- (2*sum(tau^2)) / (df*(df+2))
        Z.I <- (x[[i]]$estimate - E.I) / sqrt(V.I)
	if (x[[i]]$alternative == "two.sided") 
	    P.I <- 2 * (1 - pnorm(Z.I))
        else if (x[[i]]$alternative == "greater")
            P.I <- pnorm(Z.I, lower.tail=FALSE)
        else P.I <- pnorm(Z.I)
        Sk.I <- ((8*sum(tau^3))/(df*(df+2)*(df+4))) / (V.I^(3/2))
        Kur.I <- ((48*sum(tau^4) + 12*(sum(tau^2))^2) /
            (df*(df+2)*(df+4)*(df+6))) / (V.I^2)
	res[i,] <- c(x[[i]]$estimate, Z.I, P.I, x[[i]]$statistic,
	    x[[i]]$p.value, E.I, V.I, Sk.I, Kur.I, min.I, max.I,
	    x[[i]]$internal1)
    }
    colnames(res) <- c("Local Morans I", "Stand. dev. (N)", "Pr. (N)",
        "Saddlepoint", "Pr. (Sad)", "Expectation", "Variance",
        "Skewness", "Kurtosis", "Minimum", "Maximum",
        "omega", "sad.r", "sad.u")
    rownames(res) <- regnames
    res <- as.data.frame(res)
    res
}


summary.localmoransad <- function(object, ...) {
    res <- as.data.frame(object)
    class(res) <- c("summary.localmoransad", class(res)) 
    res
}

print.summary.localmoransad <- function(x, ...) {
	print(as.data.frame(x), ...)
	invisible(x)
}



# Copyright 2002 by Roger Bivand and Michael Tiefelsdorf
#

lm.morantest.sad <- function (model, listw, zero.policy = FALSE, 
    alternative = "greater", spChk=NULL, tol = .Machine$double.eps^0.5,
    maxiter = 1000, tol.bounds=0.0001) 
{
    if (!inherits(listw, "listw")) 
        stop(paste(deparse(substitute(listw)), "is not a listw object"))
    if (class(model) != "lm") 
        stop(paste(deparse(substitute(model)), "not an lm object"))
    N <- length(listw$neighbours)
    u <- residuals(model)
    if (N != length(u)) 
        stop("objects of different length")
    if (is.null(spChk)) spChk <- get.spChkOption()
    if (spChk && !chkIDs(u, listw))
	stop("Check of data and weights ID integrity failed")
    if (!(alternative %in% c("greater", "less", "two.sided")))
	stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
    u <- as.vector(u)
    listw.U <- listw2U(listw)
    S0 <- sum(unlist(listw.U$weights))
    lu <- lag.listw(listw.U, u, zero.policy = zero.policy)
    Nnn <- N
    if (zero.policy) Nnn <- length(which(card(listw$neighbours) > 0))
    I <- (Nnn/S0) * ((t(u) %*% lu)/(t(u) %*% u))
    p <- model$rank
    p1 <- 1:p
    XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
    X <- model.matrix(terms(model), model.frame(model))
    M <- diag(N) - X %*% XtXinv %*% t(X)
    U <- listw2mat(listw.U)
    MVM <- M %*% U %*% M
    MVM <- 0.5 * (t(MVM) + MVM)
    evalue <- eigen(MVM, only.values=TRUE)$values
    idxpos <- (which(abs(evalue) < 1.0e-7)[1]) - 1
    tau <- evalue[1:idxpos]
    tau <- c(tau, evalue[(idxpos+1+p):N])
    taumi <- tau - I
    low <- (1 / (2*taumi[length(taumi)])) + tol.bounds
    high <- (1 / (2*taumi[1])) - tol.bounds
    f <- function(omega, taumi) {sum(taumi/(1 - (2*omega*taumi)))}
    root <- uniroot(f, lower=low, upper=high, tol=tol, maxiter=maxiter,
        taumi=taumi)
    omega <- root$root
    if (omega < 0 ) sad.r <- -sqrt(sum(log(1 - 2*omega*taumi)))
    else sad.r <- sqrt(sum(log(1 - 2*omega*taumi)))
    sad.u <- omega * sqrt(2*sum(taumi^2 / (1 - (2*omega*taumi))^2))
    sad.p <- sad.r - ((1/sad.r)*log(sad.r/sad.u))
    if (alternative == "two.sided") p.sad <- 2 * pnorm(abs(sad.p), 
	lower.tail=FALSE)
    else if (alternative == "greater")
        p.sad <- pnorm(sad.p, lower.tail=FALSE)
    else p.sad <- pnorm(sad.p)
    if (p.sad < 0 || p.sad > 1) 
	warning("Out-of-range p-value: reconsider test arguments")
    statistic <- sad.p
    attr(statistic, "names") <- "Saddlepoint approximation"
    p.value <- p.sad
    estimate <- c(I)
    attr(estimate, "names") <- "Observed Moran's I"
    internal1 <- c(omega, sad.r, sad.u)
    attr(internal1, "names") <- c("omega", "sad.r", "sad.u")
    internal2 <- unlist(root)[2:4]
    attr(internal2, "names") <- c("f.root", "iter", "estim.prec")
    method <- paste("Saddlepoint approximation for global Moran's I",
        "(Barndorff-Nielsen formula)")
    data.name <- paste("\nmodel:", paste(strwrap(gsub(" *", " ", 
	    paste(deparse(model$call), sep="", collapse=""))), collapse="\n"),
    	    "\nweights: ", deparse(substitute(listw)), "\n", sep="")
    res <- list(statistic = statistic, p.value = p.value,
        estimate = estimate, method = method,
	alternative = alternative, data.name = data.name,
	internal1 = internal1, internal2 = internal2,
	df = (N-p), tau = tau)
    class(res) <- "moransad"
    return(res)
}

print.moransad <- function(x, ...) {
    class(x) <- c("htest", "moransad")
    print(x, ...)
    invisible(x)
}

summary.moransad <- function(object, ...) {
    res <- object
    tau <- object$tau
    df <- object$df
    if (length(tau) == 2) tau <- c(tau[1], rep(0, df-2), tau[2])
    max.I <- tau[1]
    min.I <- tau[length(tau)]
    E.I <- sum(tau)/df
    tau <- tau - E.I
    V.I <- (2*sum(tau^2)) / (df*(df+2))
    Z.I <- (object$estimate - E.I) / sqrt(V.I)
    Sk.I <- ((8*sum(tau^3))/(df*(df+2)*(df+4))) / (V.I^(3/2))
    Kur.I <- ((48*sum(tau^4) + 12*(sum(tau^2))^2) /
        (df*(df+2)*(df+4)*(df+6))) / (V.I^2)
    res$xtra <- c(E.I, V.I, Z.I, Sk.I, Kur.I, min.I, max.I)
    attr(res$xtra, "names") <- c("Expectation", "Variance",
        "Std. deviate", "Skewness", "Kurtosis", "Minimum", "Maximum")
    class(res) <- c("summary.moransad", "moransad")
    res
}

print.summary.moransad <- function(x, ...) {
    class(x) <- c("htest", "summary.moransad", "moransad")
    print(x, ...)
    print(c(x$xtra, x$internal1), ...)
    if (!is.null(x$internal2)) print(x$internal2, ...)
    invisible(x)
}

# Copyright 2004 by Roger Bivand 
#

nb2blocknb <- function(nb, ID, row.names = NULL) {
	if (!inherits(nb, "nb")) stop("not an nb object")
	nbNames <- as.character(attr(nb, "region.id"))
	entNames <- as.character(ID)
	if (!identical(sort(nbNames), sort(unique(entNames))))
		stop("names do not match exactly")
	n <- length(entNames)
	if (!is.null(row.names)) {
		if (length(row.names) != n) 
			stop("row.names wrong length")
		if (length(unique(row.names)) != length(row.names)) 
		stop("non-unique row.names given")
	} else {
		row.names <- as.character(1:n)
	}
	inter <- lapply(as.list(nbNames), 
		function(x) which(match(entNames, x) == 1))

	res <- vector(mode="list", length=n)
	for (i in 1:n) {
		ii <- match(entNames[i], nbNames)
		blocks <- c(ii, nb[[ii]])
		vec <- sort(unlist(inter[blocks]))
		res[[i]] <- vec[vec != i]
	}

	attr(res, "region.id") <- row.names
	class(res) <- "nb"
	attr(res, "block") <- TRUE
	attr(res, "call") <- match.call()
	res <- sym.attr.nb(res)
	invisible(res)
}


# Copyright 2001-3 by Roger S. Bivand. 
#

nb2listw <- function(neighbours, glist=NULL, style="W", zero.policy=FALSE)
{
	if(!inherits(neighbours, "nb")) stop("Not a neighbours list")
	if (!(style %in% c("W", "B", "C", "S", "U")))
		stop(paste("Style", style, "invalid"))
	n <- length(neighbours)
	cardnb <- card(neighbours)
	if (!zero.policy)
		if (any(cardnb == 0)) stop("Empty neighbour sets found")
	vlist <- vector(mode="list", length=n)
	if (is.null(glist)) {
		glist <- vector(mode="list", length=n)
		for (i in 1:n)
			if(cardnb[i] > 0) glist[[i]] <- rep(1, length=cardnb[i])
		attr(vlist, "binary") <- TRUE
	} else {
		attr(vlist, "general") <- TRUE
		source <- deparse(substitute(glist))
		attr(vlist, as.character(source)) <- TRUE
	}
	if (length(glist) != n) stop("glist wrong length")
	if (any(cardnb != unlist(lapply(glist, length))))
		stop("neighbours and glist do not conform")
	if (any(is.na(unlist(glist))))
		stop ("NAs in general weights list")
	attr(vlist, as.character(style)) <- TRUE
	if (style == "W") {
		d <- unlist(lapply(glist, sum))
		for (i in 1:n) {
			if (cardnb[i] > 0) vlist[[i]] <- (1/d[i]) * glist[[i]]
		}
	}
	if (style == "B") {
		for (i in 1:n) {
			if (cardnb[i] > 0) vlist[[i]] <- glist[[i]]
		}
	}
	if (style == "C" || style == "U") {
		D <- sum(unlist(glist))
		if (is.na(D) || !(D > 0))
			stop(paste("Failure in sum of weights:", D))
		if (zero.policy) {
			eff.n <- n - length(which(cardnb == 0))
			if (eff.n < 1) stop("No valid observations")
		} else eff.n <- n
		for (i in 1:n) {
			if (cardnb[i] > 0) {
				if (style == "C")
					vlist[[i]] <- (eff.n/D) * glist[[i]]
				else vlist[[i]] <- (1/D) * glist[[i]]
			}
		}
	}
	if (style == "S") {
		glist2 <- lapply(glist, function(x) x^2)
		q <- sqrt(unlist(lapply(glist2, sum)))
		for (i in 1:n) {
			if (cardnb[i] > 0)
				glist[[i]] <- (1/q[i]) * glist[[i]]
		}
		Q <- sum(unlist(glist))
		if (is.na(Q) || !(Q > 0))
		    stop(paste("Failure in sum of intermediate weights:", Q))
		if (zero.policy) {
			eff.n <- n - length(which(cardnb == 0))
			if (eff.n < 1) stop("No valid observations")
		} else eff.n <- n
		for (i in 1:n) {
			if (cardnb[i] > 0)
				vlist[[i]] <- (eff.n/Q) * glist[[i]]
		}
	}
	style <- style
	if (!zero.policy)
		if (any(is.na(unlist(vlist))))
			stop ("NAs in coding scheme weights list")
	res <- list(style=style, neighbours=neighbours, weights=vlist)
	class(res) <- c("listw", "nb")
	attr(res, "region.id") <- attr(neighbours, "region.id")
	attr(res, "call") <- match.call()
	if (!is.null(attr(neighbours, "GeoDa")))
		attr(res, "GeoDa") <- attr(neighbours, "GeoDa")
	if (!is.null(attr(res, "GeoDa")$dist)) 
		attr(res, "GeoDa")$dist <- NULL
	invisible(res)
}

# Copyright 2001-3 by Roger Bivand
#


nb2mat <- function(neighbours, glist=NULL, style="W", zero.policy=FALSE)
{
	if(class(neighbours) != "nb") stop("Not a neighbours list")
	listw <- nb2listw(neighbours, glist=glist, style=style,
		zero.policy=zero.policy)
	res <- listw2mat(listw)
	attr(res, "call") <- match.call()
	invisible(res)
}

listw2mat <- function(listw) {
	n <- length(listw$neighbours)
	cardnb <- card(listw$neighbours)
	if (any(is.na(unlist(listw$weights))))
		stop ("NAs in general weights list")
	res <- matrix(0, nrow=n, ncol=n)
	for (i in 1:n)
	    if (cardnb[i] > 0)
		res[i, listw$neighbours[[i]]] <- listw$weights[[i]]
	if (!is.null(attr(listw, "region.id")))
		row.names(res) <- attr(listw, "region.id")
	invisible(res)
}

invIrM <- function(neighbours, rho, glist=NULL, style="W") {
	if(class(neighbours) != "nb") stop("Not a neighbours list")
	n <- length(neighbours)
	V <- nb2mat(neighbours, glist, style)
	e <- eigen(V, only.values = TRUE)$values
	if (is.complex(e)) feasible <- 1/(range(Re(e)))
	else feasible <- 1/(range(e))
	if (rho <= feasible[1] || rho >= feasible[2])
		stop(paste("Rho outside feasible range:", feasible))
	mat <- diag(n) - rho * V
	res <- solve(mat)
	attr(res, "call") <- match.call()
	invisible(res)
}

invIrW <- function(listw, rho) {
	if(!inherits(listw, "listw")) stop("Not a weights list")
	n <- length(listw$neighbours)
	V <- listw2mat(listw)
	e <- eigen(V, only.values = TRUE)$values
	if (is.complex(e)) feasible <- 1/(range(Re(e)))
	else feasible <- 1/(range(e))
	if (rho <= feasible[1] || rho >= feasible[2])
		stop(paste("Rho outside feasible range:", feasible))
	mat <- diag(n) - rho * V
	res <- solve(mat)
	attr(res, "call") <- match.call()
	invisible(res)
}

mat2listw <- function(x, row.names=NULL) {
	if (!is.matrix(x)) stop("x is not a matrix")
	n <- nrow(x)
	m <- ncol(x)
	if (n != m) stop("x must be a square matrix")
	if (any(x < 0)) stop("values in x cannot be negative")
	if (any(is.na(x))) stop("NA values in x not allowed")
    	if (!is.null(row.names)) {
		if(length(row.names) != n)
            		stop("row.names wrong length")
		if (length(unique(row.names)) != length(row.names))
	    		stop("non-unique row.names given")
    	}
    	if (is.null(row.names)) {
		if (!is.null(row.names(x))) {
			row.names <- row.names(x)
		} else {
			row.names <- as.character(1:n)
		}
	}
	style <- "M"
	neighbours <- vector(mode="list", length=n)
	weights <- vector(mode="list", length=n)
	for (i in 1:n) {
		nbs  <- which(x[i,] > 0.0)
		if (length(nbs) > 0) {
			neighbours[[i]] <- nbs
			weights[[i]] <- x[i, nbs]
		} else {
			neighbours[[i]] <- 0
		}
	}
	class(neighbours) <- "nb"
	attr(neighbours, "region.id") <- row.names
 	attr(neighbours, "call") <- NA
        attr(neighbours, "sym") <- is.symmetric.nb(neighbours, 
		verbose=FALSE, force=TRUE)
	res <- list(style=style, neighbours=neighbours, weights=weights)
	class(res) <- c("listw", "nb")
	attr(res, "region.id") <- attr(neighbours, "region.id")
	attr(res, "call") <- match.call()
	invisible(res)	
}
# Copyright 2001-3 by Roger Bivand
#


nbdists <- function(nb, coords, lonlat=FALSE) {
	if (!inherits(nb, "nb")) 
        	stop("Not a neighbours list")
	if (!is.numeric(coords)) stop("Data non-numeric")
	if (!is.matrix(coords)) 
            stop("Data not in matrix form")
        if (any(is.na(coords))) 
            stop("Data include NAs")
	if (!is.double(coords)) storage.mode(coords) <- "double"
	n.nb <- length(nb)
	np <- nrow(coords)
        if (np != n.nb) 
            stop("Number of coords not equal to number of regions")
        dimension <- ncol(coords)
        dlist <- .Call("nbdists", nb, as.matrix(coords), as.integer(np), 
            as.integer(dimension), as.integer(lonlat), PACKAGE="spdep")
	attr(dlist[[1]], "call") <- match.call()
	invisible(dlist[[1]])
}

# Copyright 2001-2002 by Roger Bivand
#


nblag <- function(neighbours, maxlag)
{
	if (!is.null(attributes(neighbours)$self.included) &&
		(as.logical(attributes(neighbours)$self.included)))
		stop("No lags for neighbours lists including self")
	n <- length(neighbours)
	lags <- vector(mode="list", length=maxlag)
	lags[[1]] <- neighbours
	cds <- card(neighbours)
	for (thislag in 2:maxlag)
		lags[[thislag]] <- vector(mode="list", length=n)
	for (i in 1:n) {
		already <- i
		new <- neighbours[[i]]
		for (thislag in 2:maxlag) {
			if (cds[i] > 0) {
				already <- c(already, new)
				active <- new
				new <- NULL
				for (j in active)
					new <- c(new, neighbours[[j]])
				new <- sort(unique(new))
				res <- new[-which(new %in% already)]
				if(length(res) == 0) lags[[thislag]][[i]] <- 0
				else lags[[thislag]][[i]] <- res
			}
			else lags[[thislag]][[i]] <- 0
		}
	}
	for (i in 1:maxlag) {
		class(lags[[i]]) <- "nb"
		attr(lags[[i]], "region.id") <- attr(neighbours, "region.id")
		lags[[i]] <- sym.attr.nb(lags[[i]])
	}
	attr(lags, "call") <- match.call()
	invisible(lags)
}
# Copyright 2001-2002 by Nicholas Lewin-Koh and Roger Bivand
#


union.nb<-function(nb.obj1, nb.obj2){
  if(!inherits(nb.obj1,"nb") | !inherits(nb.obj2,"nb")){
    stop("Both arguments must be of class nb")
  }
  if(any(attr(nb.obj1,"region.id")!= attr(nb.obj2,"region.id"))){
   stop("Both neighbor objects must be \n generated from the same coordinates")
  }
  n <- length(nb.obj1)
  if (n != length(nb.obj2)) stop("Both arguments must be of same length")
  card1 <- card(nb.obj1)
  card2 <- card(nb.obj2)
  new.nb<-vector(mode="list", length=n)
  for(i in 1:n) {
    if (card1[i] == 0) {
      if (card2[i] == 0) new.nb[[i]] <- 0
      else new.nb[[i]] <- nb.obj2[[i]]
    } else {
      if (card2[i] == 0) new.nb[[i]] <- nb.obj1[[i]]
      else new.nb[[i]]<-sort(union(nb.obj1[[i]], nb.obj2[[i]]))
    }
  }
  attr(new.nb,"region.id")<-attr(nb.obj1,"region.id")
  attr(new.nb,"type")<-paste("union(",attr(nb.obj1,"type"),
                             ",",attr(nb.obj2,"type"),")")
  class(new.nb)<-"nb"
  new.nb
 }

intersect.nb<-function(nb.obj1, nb.obj2){
  if(!inherits(nb.obj1,"nb") | !inherits(nb.obj2,"nb")){
    stop("Both arguments must be of class nb")
  }
  if(any(attr(nb.obj1,"region.id")!= attr(nb.obj2,"region.id"))){
   stop("Both neighbor objects must be \n generated from the same coordinates")
  }
  n <- length(nb.obj1)
  if (n != length(nb.obj2)) stop("Both arguments must be of same length")
  card1 <- card(nb.obj1)
  card2 <- card(nb.obj2)
  new.nb<-vector(mode="list", length=n)
  for(i in 1:n) {
    if (card1[i] > 0 && card2[i] > 0) {
      res <- sort(intersect(nb.obj1[[i]], nb.obj2[[i]]))
      if(length(res) == 0) new.nb[[i]] <- 0
      else new.nb[[i]] <- res
    } else new.nb[[i]] <- 0
  }
  attr(new.nb,"region.id")<-attr(nb.obj1,"region.id")
  attr(new.nb,"type")<-paste("intersect(",attr(nb.obj1,"type"),
                             ",",attr(nb.obj2,"type"),")")
  class(new.nb)<-"nb"
  new.nb
}
setdiff.nb<-function(nb.obj1, nb.obj2){
  	if(!inherits(nb.obj1,"nb") | !inherits(nb.obj2,"nb")){
    		stop("Both arguments must be of class nb")
  	}
  	if(any(attr(nb.obj1,"region.id")!= attr(nb.obj2,"region.id"))){
   		stop("Both neighbor objects must be \n generated from the same coordinates")
  	}
  	n <- length(nb.obj1)
  	if (n != length(nb.obj2)) stop("Both arguments must be of same length")
  	card1 <- card(nb.obj1)
  	card2 <- card(nb.obj2)
  	new.nb<-vector(mode="list", length=n)
  	for(i in 1:n) {
    		if (card1[i] == 0) {
      			if (card2[i] == 0) new.nb[[i]] <- 0
      			else new.nb[[i]] <- nb.obj2[[i]]
    		} else {
      			if (card2[i] == 0) new.nb[[i]] <- nb.obj1[[i]]
      			else {
            			if (card2[i] == 0)
                			new.nb[[i]] <- nb.obj1[[i]]
            			else {
                			if (card1[i] >= card2[i]) {
                  				a <- nb.obj1[[i]]
                  				b <- nb.obj2[[i]]
                			} else {
                  				b <- nb.obj1[[i]]
                  				a <- nb.obj2[[i]]
                			}
                			res <- sort(setdiff(a, b))
			                if(length(res) == 0) new.nb[[i]] <- 0
                			else new.nb[[i]] <- res
	    			}
        		}
    		}
  	}
  	attr(new.nb,"region.id")<-attr(nb.obj1,"region.id")
  	attr(new.nb,"type")<-paste("setdiff(",attr(nb.obj1,"type"),
                             ",",attr(nb.obj2,"type"),")")
  	class(new.nb)<-"nb"
  	new.nb
}

complement.nb<-function(nb.obj){
   if(!inherits(nb.obj,"nb")){
    stop("Argument must be of class nb")
   }
  n <- length(nb.obj)
  card1 <- card(nb.obj)
  new.nb<-vector(mode="list", length=n)
  cmp<-1:n
  attributes(new.nb)<-attributes(nb.obj)
  for(i in 1:n) {
    if (card1[i] == 0) new.nb[[i]] <- cmp
    else {
      res <- sort(cmp[-nb.obj[[i]]])
      if(length(res) == 0) new.nb[[i]] <- 0
      else new.nb[[i]] <- res
    }
  }
  attr(new.nb,"type")<-paste("complement(",attr(nb.obj,"type"),")")
  class(new.nb)<-"nb"
  new.nb
 }
# Copyright 2004 by Roger Bivand and Danlin Yu
#

p.adjustSP <- function(p, nb, method="none") {
	if(class(nb) != "nb") stop("Not a neighbours list")
        n <- card(nb) + 1
        pn <- cbind(p, n)
        res <- apply(pn, 1, function(x) p.adjust(x[1], method=method, n=x[2]))
        res
}

# Copyright 2001 by Roger Bivand
#


plot.nb <- function(x, coords, col="black", points=TRUE, add=FALSE, 
	arrows=FALSE, length=0.1, ...) {
	nb <- x
	sym <- is.symmetric.nb(nb, verbose = FALSE, force = FALSE)
	x <- coords[,1]
	y <- coords[,2]
	n <- length(nb)
	xlim <- range(x)
	ylim <- range(y)
	if (!add) {
		plot.new()
        	plot.window(xlim = xlim, ylim = ylim, log="", asp=1)
	}
	cardnb <- card(nb)
	if (length(col) < n) col <- rep(col[1], n)
	for (i in 1:n) {
		if (cardnb[i] > 0) {
        		inb <- nb[[i]]
        		for (j in inb) {
				if (sym) {
					lines(c(x[i], x[j]), c(y[i], y[j]),
						col=col[i], ...)
				} else {
					if (arrows) 
						arrows(x[i], y[i], x[j], y[j], 
						col=col[i], length=length, ...)
					else lines(c(x[i], x[j]), c(y[i], y[j]),
						col=col[i], ...)
				}

			}
		}
	}
	if (points) points(x, y, ...)
}
# Copyright 2001-2004 by Roger Bivand 
#
	


poly2nb <- function(pl, row.names=NULL, snap=sqrt(.Machine$double.eps),
	queen=TRUE) {
	if (!inherits(pl, "polylist")) stop("Not a polygon list")
	if (inherits(pl, "multiparts")) stop("Convert to newer polylist format")
	n <- length(pl)
	regid <- attr(pl, "region.id")
	if (is.null(regid)) {
		if(is.null(row.names)) regid <- as.character(1:n)
		else {
			if(length(row.names) != n)
				stop("row.names wrong length")
			else if (length(unique(row.names)) != length(row.names))
	    			stop("non-unique row.names given")
			else regid <- row.names
		}
	}
	poly2bbs <- function(pl) {
		n <- length(pl)
		res <- matrix(0, nrow=n, ncol=4)
		for (i in 1:n) res[i,] <- attr(pl[[i]], "bbox")
		res
	}
	bb <- poly2bbs(pl)
	


	polypoly2 <- function(poly1, poly2, snap) {
		n1 <- nrow(poly1)
		n2 <- nrow(poly2)
		res <- .Call("polypoly", as.double(poly1), 
			as.integer(n1), as.double(poly2), 
			as.integer(n2), as.double(snap), PACKAGE="spdep")
		res
	}

	ans <- vector(mode="list", length=n)
	for (i in 1:n) ans[[i]] <- integer(0)
	criterion <- ifelse(queen, 0, 1)
	for (i in 1:(n-1)) {
		for (j in (i+1):n) {
			jhit <- .Call("spInsiders", as.double(bb[i,]), 
				as.double(bb[j,]), PACKAGE="spdep")
			if (jhit > 0) {
			    khit <- 0
			    khit <- polypoly2(na.omit(pl[[i]]), 
				na.omit(pl[[j]]), snap)

			    if (khit > criterion) {
				ans[[i]] <- c(ans[[i]], j)
				ans[[j]] <- c(ans[[j]], i)
			    }
			}
		}
	}
	for (i in 1:n) ans[[i]] <- sort(ans[[i]])
	class(ans) <- "nb"
	attr(ans, "region.id") <- regid
	attr(ans, "call") <- match.call()
	if (queen) attr(ans, "type") <- "queen"
	else attr(ans, "type") <- "rook"
	invisible(ans)
}


# Copyright 2001-2003 by Roger Bivand and Luc Anselin
#

read.gal <- function(file, region.id=NULL, override.id=FALSE) 
{
	con <- file(file, open="r")
	line <- unlist(strsplit(readLines(con, 1), " "))
	x <- subset(line, nchar(line) > 0)
	if (length(x) == 1) {
		n <- as.integer(x[1])
		shpfile <- as.character(NA)
		ind <- as.character(NA)
	} else if (length(x) == 4) {
		n <- as.integer(x[2])
		shpfile <- as.character(x[3])
		ind <- as.character(x[4])
	} else stop ("Invalid header line in GAL file")
	if (n < 1) stop("Non-positive number of regions")
	if (!is.null(region.id))
		if (length(unique(region.id)) != length(region.id))
	    		stop("non-unique region.id given")
	if (is.null(region.id)) region.id <- as.character(1:n)
	if (n != length(region.id))
		stop("Mismatch in dimensions of GAL file and region.id")
    	rn <- character(n)
	res <- vector(mode="list", length=n)
	for (i in 1:n) {
		line <- unlist(strsplit(readLines(con, 1), " "))
		x <- subset(line, nchar(line) > 0)
		rn[i] <- x[1]
		line <- unlist(strsplit(readLines(con, 1), " "))
		y <- subset(line, nchar(line) > 0)
		if(length(y) != as.integer(x[2])) {
			close(con)
			stop(paste("GAL file corrupted at region", i))
		}
		res[[i]] <- y
	}
	close(con)
	if (!override.id) {
		if (!all(rn %in% region.id)) {
			stop("GAL file IDs and region.id differ")
		}
	} else region.id <- rn
	mrn <- match(rn, region.id)
	res1 <- vector(mode="list", length=n)
	for (i in 1:n) {
		x <- match(res[[i]], region.id)
		if (any(is.na(x)) | (length(x) != length(res[[i]]))) {
			stop(paste("GAL file corrupted at region", i))
		}
		if(any(x < 0) || any(x > n))
			stop("GAL file corrupted")

		res1[[mrn[i]]] <- sort(x)
	}
	class(res1) <- "nb"
    	attr(res1, "region.id") <- region.id
	attr(res1, "GeoDa") <- list(shpfile=shpfile, ind=ind)
	attr(res1, "gal") <- TRUE
	attr(res1, "call") <- TRUE
	res1 <- sym.attr.nb(res1)
	invisible(res1)
}

write.nb.gal <- function(nb, file, oldstyle=TRUE, shpfile=NULL, ind=NULL) {
	if(class(nb) != "nb") stop("not an nb object")
	n <- length(nb)
	rn <- attr(nb, "region.id")
	if (is.null(shpfile)) {
		tmp <- attr(nb, "GeoDa")$shpfile
		if (is.null(tmp)) shpfile <- "unknown"
		else shpfile <- tmp
	}
	if (is.null(ind)) {
		tmp <- attr(nb, "GeoDa")$ind
		if (is.null(tmp)) ind <- "unknown"
		else ind <- tmp
	}
	con <- file(file, open="w")
	if (oldstyle) writeLines(paste(n), con)
	else writeLines(paste("0", n, shpfile, ind, sep=" "), con)
	for (i in 1:n) {
		if (oldstyle) 
			writeLines(paste(i, length(nb[[i]]), 
				collapse=" "), con)
		else writeLines(paste(rn[i], length(nb[[i]]),
			collapse=" "), con)
		if (oldstyle) writeLines(paste(nb[[i]], collapse=" "), con)
		else writeLines(paste(rn[nb[[i]]], collapse=" "), con)
	}
	close(con)
}


# read.geoda
# helper function to read GeoDa export files
# LA 7/10/03
# specify input file = file
# default is no row names, specify row names as second parameter
# example: balt <- read.geoda("baltim.txt")
#          balt <- read.geoda("baltim.txt","STATION")

read.geoda <- function(file, row.names=NULL, skip=0)
{
	read.csv(file=file, header=TRUE, skip=skip, row.names=row.names)
	
}

# Copyright 2003-4 by Luc Anselin and Roger Bivand
#

# LA 6/28/03 read.gwt
# LA 7/12/03 revised, sorted ids
# LA 9/30/03 use match to correct orders

read.gwt2nb <- function(file, region.id=NULL) {
	con <- file(file, open="r")   #opens the file
	firstline <- unlist(strsplit(readLines(con,1)," "))
	if (length(firstline) == 4) {
		n <- as.integer(firstline[2])
		shpfile <- firstline[3]
		ind <- firstline[4]
		if (ind != deparse(substitute(region.id)))
			warning(paste("region.id not named", ind))
	} else if (length(firstline) == 1) {
		n <- as.integer(firstline[1])
		shpfile <- as.character(NA)
		ind <- as.character(NA)
		warning("Old-style GWT file")
	} else stop("Invalid header line format for GWT file")
	close(con)
	nseq <- 1:n
	if (is.null(region.id)) region.id <- nseq
	if (n != length(region.id))
		stop("Mismatch in dimensions of GWT file and region.id")
	if (length(unique(region.id)) != length(region.id))
	    stop("non-unique region.id given")
	odij <- read.table(file, skip=1)
	    # convert region.id to order
	regodij <- match(odij[,1], region.id)
	regddij <- match(odij[,2], region.id)
	odij <- cbind(regodij, regddij, odij[,3])
	qorder <- order(odij[,1],odij[,2])
	odij <- odij[qorder,]
	origvec <- unique(odij[,1])
	if (!all(nseq %in% origvec))
		warning(paste(paste(region.id[which(!(nseq %in% origvec))], 
			collapse=", "), "are not origins"))
	destvec <- unique(odij[,2])
	if (!all(nseq %in% destvec))
		warning(paste(paste(region.id[which(!(nseq %in% destvec))], 
			collapse=", "), "are not destinations"))

	res <- vector(mode="list", length=n)
	vlist <- vector(mode="list", length=n)
	rle.sn <- rle(odij[,1])
	cs1.sn <- cumsum(rle.sn$lengths)
	cs0.sn <- c(1, cs1.sn[1:(n-1)]+1)
	ii <- 1
	for (i in 1:n) {
		if (rle.sn$value[ii] == i) {
			res[[i]] <- as.integer(odij[cs0.sn[ii]:cs1.sn[ii],2])
			vlist[[i]] <- as.double(odij[cs0.sn[ii]:cs1.sn[ii],3])
			ii <- ii+1
		} else {
			res[[i]] <- as.integer(0)
		}
	}

	class(res) <- c("nb", "GWT")
	attr(res, "region.id") <- region.id
	attr(res, "neighbours.attrs") <- as.character(NA)
	attr(res, "weights.attrs") <- as.character(NA)
	attr(res, "GeoDa") <- list(dist=vlist, shpfile=shpfile, ind=ind)
	attr(res, "call") <- match.call()
	attr(res, "n") <- n
	invisible(res)
}

write.sn2gwt <- function(sn, file, shpfile=NULL, ind=NULL) {
	if(!inherits(sn, "spatial.neighbour")) 
	    stop("not a spatial.neighbour object")
	n <- attr(sn, "n")
	if (is.null(shpfile)) {
		tmp <- attr(sn, "GeoDa")$shpfile
		if (is.null(tmp)) shpfile <- "unknown"
		else shpfile <- tmp
	}
	if (is.null(ind)) {
		tmp <- attr(sn, "GeoDa")$ind
		if (is.null(tmp)) ind <- "unknown"
		else ind <- tmp
	}
	con <- file(file, open="w")
	writeLines(paste("0", n, shpfile, ind, sep=" "), con)
	write.table(as.data.frame(sn), file=con, append=TRUE,
		row.names=FALSE, col.names=FALSE)
	close(con)
}

write.sn2dat <- function(sn, file) {
	if(!inherits(sn, "spatial.neighbour")) 
	    stop("not a spatial.neighbour object")
	write.table(data.frame(sn[order(sn[,2]), ]), 
	    file=file, col.names=FALSE, row.names=FALSE)
}

read.dat2listw <- function(file) {
	wmat <- read.table(file)
	sn <- wmat[order(wmat[,1]),]
	IDS <- unique(sn[,1])
	class(sn) <- c("spatial.neighbour", "data.frame")
	attr(sn, "n") <- length(IDS)
	attr(sn, "region.id") <- as.character(IDS)
	listw <- sn2listw(sn)
	invisible(listw)
}

# Copyright 2001 by Nicholas Lewin-Koh 
#

relativeneigh <- function(coords) {
    x <- coords
    if (!is.matrix(x)) stop("Data not in matrix form")
    if (any(is.na(x))) stop("Data cannot include NAs")
    np <- nrow(x)
    if(ncol(x)!=2) stop("Only works in 2d")
    g1<-g2<-rep(0,np*3)
    nogab <- 0
    z <- .C("compute_relative", np=as.integer(np), from=as.integer(g1),
             to=as.integer(g2), nedges=as.integer(nogab),
             x=as.double(x[,1]), y=as.double(x[,2]), PACKAGE="spdep")
    z$from<-z$from[1:z$nedges]
    z$to<-z$to[1:z$nedges]
    attr(z, "call") <- match.call()
    class(z)<-c("Graph","relative")
    invisible(z)
}

plot.relative<-function(x, show.points=FALSE, add=FALSE, linecol=par(col),...){
  if(!add) plot(rel$x,rel$y,type='n')
  segments(rel$x[rel$from],rel$y[rel$from],
           rel$x[rel$to],rel$y[rel$to],col=linecol)
  if(show.points) points(rel$x,rel$y,...)
}
# Copyright 2002-4 by Roger Bivand
#

residuals.sarlm <- function(object, ...) {
	if (is.null(object$na.action))
		object$residuals
	else napredict(object$na.action, object$residuals)
}

deviance.sarlm <- function(object, ...) {
	deviance(object$lm.target)
}

coef.sarlm <- function(object, ...) {
	ret <- object$coefficients
	if(object$type == "error") ret <- c(ret, object$lambda)
	else ret <- c(ret, object$rho)
	ret
}

fitted.sarlm <- function(object, ...) {
	if (is.null(object$na.action))
		object$fitted.values
	else napredict(object$na.action, object$fitted.values)
}

predict.sarlm <- function(object, newdata=NULL, listw=NULL, 
	zero.policy=FALSE, ...) {
	if (is.null(newdata)) {
		res <- fitted.values(object)
		X <- model.matrix(terms(object$lm.model), 
				model.frame(object$lm.model))
		B <- coefficients(object$lm.target)
		y <- model.response(model.frame(object$lm.model))
		tarX <- model.matrix(terms(object$lm.target), 
				model.frame(object$lm.target))
		tary <- model.response(model.frame(object$lm.target))
		if (object$type == "error") {
			attr(res, "trend") <- as.vector(X %*% B)
			attr(res, "signal") <- as.vector( -1 * (tary - y) - 					-1 * (tarX - X) %*% B)
		} else {
			attr(res, "trend") <- fitted.values(object$lm.target)
			attr(res, "signal") <- as.vector( -1 * (tary - y))
		}
	}
	else {
		if (object$type == "error") {
			B <- coefficients(object$lm.target)
			tt <- terms(object$lm.model) 
			X <- model.matrix(delete.response(tt), data=newdata)
			trend <- X %*% B
			signal <- rep(0, length(trend))
			res <- trend + signal
			attr(res, "trend") <- trend
			attr(res, "signal") <- signal
		} else if (object$type == "mixed") {
			if (is.null(listw) || !inherits(listw, "listw")) 
				stop ("spatial weights list required")
			if (nrow(newdata) != length(listw$neighbours))
				stop("mismatch between newdata and spatial weights")
			B <- coefficients(object$lm.target)
			mt <- terms(object$formula, data = newdata)
			mf <- lm(object$formula, newdata, method="model.frame")
			X <- model.matrix(mt, mf)
			K <- ifelse(colnames(X)[1] == "(Intercept)", 2, 1)
			m <- ncol(X)
			WX <- matrix(nrow=nrow(X),ncol=(m-(K-1)))
			for (k in K:m) {
				wx <- lag.listw(listw, X[,k], 
				    zero.policy=zero.policy)
				if (any(is.na(wx))) 
				    stop("NAs in lagged independent variable")
				WX[,(k-(K-1))] <- wx
			}
			X <- cbind(X, WX)
			trend <- X %*% B
			signal <- object$rho * lag.listw(listw, 
				(invIrW(listw, object$rho) %*% trend), 
				zero.policy=zero.policy)
			res <- trend + signal
			attr(res, "trend") <- trend
			attr(res, "signal") <- signal

		} else {
			if (is.null(listw) || !inherits(listw, "listw")) 
				stop ("spatial weights list required")
			if (nrow(newdata) != length(listw$neighbours))
				stop("mismatch between newdata and spatial weights")
			B <- coefficients(object$lm.target)
			mt <- terms(object$formula, data = newdata)
			mf <- lm(object$formula, newdata, method="model.frame")
			X <- model.matrix(mt, mf)
			trend <- X %*% B
			raw.sig <- invIrW(listw, object$rho) %*% trend
			signal <- object$rho * lag.listw(listw, 
				raw.sig, zero.policy=zero.policy)
			res <- trend + signal
			attr(res, "trend") <- trend
			attr(res, "signal") <- signal
		}
	}
	class(res) <- "sarlm.pred"
	res
}

print.sarlm.pred <- function(x, ...) {
	res <- data.frame(fit=as.vector(x), trend=attr(x, "trend"),
		signal=attr(x, "signal"))
	print(res, ...)
	invisible(res)
}

# Copyright 2001 by Nicholas Lewin-Koh
#


soi.graph <- function(tri.nb,coords){
  x <- coords
  if (!is.matrix(x)) stop("Data not in matrix form")
  if (any(is.na(x))) stop("Data cannot include NAs")
  np<-length(tri.nb)
  noedges<-0
  rad<-nearneigh<-rep(0,np)
  neigh<-unlist(tri.nb)  
  noneigh<-unlist(lapply(tri.nb,length))
  g1<-g2<-rep(0,sum(noneigh))
  answ<-.C("compute_soi", np=as.integer(np), from=as.integer(g1),
     to=as.integer(g2), nedges=as.integer(noedges),
     notri.nb=as.integer(noneigh), tri.nb=as.integer(neigh),
     nn=as.integer(nearneigh), 
     circles=as.double(rad), x=as.double(x[,1]), y=as.double(x[,2]),
     PACKAGE="spdep")
  answ$from<-answ$from[1:answ$nedges]
  answ$to<-answ$to[1:answ$nedges]
  answ<-list(np=answ$np,nedges=answ$nedges,
             from=answ$from,to=answ$to,circles=answ$circ)
  attr(answ, "call") <- match.call()
  class(answ)<-c("Graph","SOI")
  invisible(answ)
}
# Copyright 2002-3 by Roger Bivand
#

sp.correlogram <- function (neighbours, var, order = 1, method = "corr", 
    style = "W", randomisation = TRUE, zero.policy = FALSE, spChk = NULL) {
    if (class(neighbours) != "nb") 
        stop("not a neighbours list")
    if (any(is.na(var))) 
        stop("no NAs permitted in variable")
    if (is.null(spChk)) 
        spChk <- get.spChkOption()
    if (spChk && !chkIDs(var, nb2listw(neighbours, zero.policy = zero.policy))) 
        stop("Check of data and weights ID integrity failed")
    nblags <- nblag(neighbours, maxlag = order)
    cardnos <- vector(mode = "list", length = order)
    for (i in 1:order) cardnos[[i]] <- table(card(nblags[[i]]))
    if (method == "corr") {
        lags.x <- matrix(0, nrow = length(var), ncol = order)
        for (i in 1:order) lags.x[, i] <- lag.listw(nb2listw(nblags[[i]], 
            style = style, zero.policy = zero.policy), var, zero.policy = zero.policy)
        res <- cor(cbind(var, lags.x))[1, -1]
        names(res) <- 1:order
    }
    else if (method == "I") {
        res <- matrix(NA, nrow = order, ncol = 3)
        for (i in 1:order) {
            listw <- nb2listw(nblags[[i]], style = style, zero.policy = zero.policy)
            res[i,] <- moran.test(var, listw, randomisation = randomisation, 
		zero.policy = zero.policy)$estimate
        }
        rownames(res) <- 1:order
    }
    else stop("method unknown")
    obj <- list(res=res, method=method, cardnos=cardnos, var=deparse(substitute(var)))
    class(obj) <- "spcor"
    obj
}


#function(neighbours, var, order=1, method="corr",
#	style="W", zero.policy=FALSE, spChk=NULL) {
#	if (class(neighbours) != "nb") stop("not a neighbours list")
#	if (any(is.na(var))) stop("no NAs permitted in variable")
#	if (is.null(spChk)) spChk <- get.spChkOption()
#	if (spChk && !chkIDs(var, nb2listw(neighbours, 
#		zero.policy=zero.policy)))
#		stop("Check of data and weights ID integrity failed")
#	nblags <- nblag(neighbours, maxlag=order)
#	cardnos <- vector(mode="list", length=order)
#	for (i in 1:order) cardnos[[i]] <- table(card(nblags[[i]]))
#	if (method == "corr") {
#		lags.x <- matrix(0, nrow=length(var), ncol=order)
#		for (i in 1:order) lags.x[,i] <- lag.listw(nb2listw(nblags[[i]],
#			style=style, zero.policy=zero.policy), var, 
#			zero.policy=zero.policy)
#		res <- cor(cbind(var, lags.x))[1,-1]
#	} else if (method == "I") {
#		res <- numeric(length=order)
#		for (i in 1:order) {
#			listw <- nb2listw(nblags[[i]], style=style, 
#				zero.policy=zero.policy) 
#			S0 <- Szero(listw)
#			res[i] <- moran(var, listw, length(var), S0, 
#				zero.policy=zero.policy)$I
#		}			
#	} else stop("method unknown")
#	class(res) <- "spcor"
#	names(res) <- 1:order
#	attr(res, "method") <- method
#	attr(res, "cardnos") <- cardnos
#	attr(res, "var") <- deparse(substitute(var))
#	res
#}

print.spcor <- function (x, ...) 
{
    if (x$method == "I") {
        meth <- "Moran's I"
        res <- as.matrix(x$res)
        rownames(res) <- rownames(x$res)
        colnames(res) <- c("estimate", "expectation", "variance")
    } else {
        meth <- "Spatial autocorrelation"
        res <- as.vector(x$res)
        names(res) <- names(x$res)
    }
    cat("Spatial correlogram for", x$var, "\nmethod:", 
        meth, "\n")
    print(res)
    invisible(res)
}


#function(x, ...) {
#	res <- as.vector(x)
#	names(res) <- names(x)
#	if (attr(x, "method") == "I") {
#		meth <- "Moran's I"
#	} else {
#		meth <- "Spatial autocorrelation"
#	}
#	cat("Spatial correlogram for", attr(x, "var"), "\nmethod:", meth, "\n")
#	print(res)
#	invisible(res)
#}

plot.spcor <- function (x, main, ylab, ylim, ...) 
{
    if (missing(main)) 
        main <- x$var
    if (x$method == "I") {
        lags <- as.integer(rownames(x$res))
        sd2 <- sqrt(x$res[,3])
        if (missing(ylim)) {
            ylim <- range(c(x$res[,1]+sd2, x$res[,1]-sd2))
	}
        if (missing(ylab)) 
            ylab <- "Moran's I"
        plot(lags, x$res[,1], type="p", pch=18, ylim = ylim, main = main, ylab = ylab, xaxt = "n")
#        segments(lags, x$res[,1], lags, x$res[,2], lwd=4, col="grey")
        arrows(lags, x$res[,1]+sd2, lags, x$res[,1]-sd2, length=0.1, angle=90)
        arrows(lags, x$res[,1]-sd2, lags, x$res[,1]+sd2, length=0.1, angle=90)
        axis(1, at = lags)
        abline(h = x$res[1,2])
    }
    else {
        res <- as.vector(x$res)
        lags <- as.integer(names(x$res))
        if (missing(ylim)) 
            ylim <- c(-1, 1)
        if (missing(ylab)) 
            ylab <- "Spatial autocorrelation"
        plot(lags, res, type = "h", ylim = ylim, main = main, ylab = ylab, 
            lwd = 4, xaxt = "n")
        axis(1, at = lags)
        abline(h = 0)
    }
}

#function(x, main, ylab, ylim, ...) {
#	res <- as.vector(x)
#	lags <- as.integer(names(x))
#	if (missing(main)) main <- attr(x, "var")
#	if (attr(x, "method") == "I") {
#		if (missing(ylim)) ylim <- range(res)
#		if (missing(ylab)) ylab <- "Moran's I"
#	} else {
#		if (missing(ylim)) ylim <- c(-1,1)
#		if (missing(ylab)) ylab <- "Spatial autocorrelation"
#	}
#	plot(lags, res, type="h", ylim=ylim, main=main, ylab=ylab, 
#		lwd=4, xaxt="n")
#	axis(1, at=lags)
#	abline(h=0)
#}

# Copyright 2002 by Roger Bivand 
#

sp.mantel.mc <- function(var, listw, nsim, type="moran", zero.policy=FALSE,
	alternative="greater", spChk=NULL) {
	alternative <- match.arg(alternative, c("greater", "less"))
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if(!is.numeric(var)) stop(paste(deparse(substitute(var)),
		"is not a numeric vector"))
	if(missing(nsim)) stop("nsim must be given")
	n <- length(listw$neighbours)
	if(nsim > gamma(n+1)) stop("nsim too large for this n")
	if (any(is.na(var))) stop("NA in var")
	if (n != length(var)) stop("objects of different length")
	if (is.null(spChk)) spChk <- get.spChkOption()
	if (spChk && !chkIDs(var, listw))
		stop("Check of data and weights ID integrity failed")
	
	listw.U <- listw2U(listw)
	mantel.moran <- function(x, listwU, zero.policy) {
		res <- x * lag.listw(listw.U, x, zero.policy=zero.policy)
		res <- sum(res)
		res
	}
	mantel.geary <- function(x, listwU, zero.policy) {
		res <- geary.intern(x, listwU, length(x), 
			zero.policy=zero.policy, type="geary")
		res <- sum(res)
		res
	} 
	mantel.sokal <- function(x, listwU, zero.policy) {
		res <- geary.intern(x, listwU, length(x), 
			zero.policy=zero.policy, type="sokal")
		res <- sum(res)
		res
	}
	if (type == "moran") f <- mantel.moran
	else if (type == "geary") f <- mantel.geary
	else if (type == "sokal") f <- mantel.sokal
	else stop("unknown type")
	xs <- scale(var)
	res <- numeric(length=nsim+1)
	for (i in 1:nsim) {
		y <- sample(xs)
		res[i] <- f(y, listw.U, zero.policy=zero.policy)
	}
	res[nsim+1] <- f(xs, listw.U, zero.policy=zero.policy)
	rankres <- rank(res)
	xrank <- rankres[length(res)]
	diff <- nsim - xrank
	diff <- ifelse(diff > 0, diff, 0)
        if (alternative == "less") 
        	pval <- punif((diff + 1)/(nsim + 1), lower.tail=FALSE)
    	else if (alternative == "greater") 
        	pval <- punif((diff + 1)/(nsim + 1))
	if (pval < 0 || pval > 1) 
		warning("Out-of-range p-value: reconsider test arguments")
	statistic <- res[nsim+1]
	names(statistic) <- "statistic"
	parameter <- xrank
	names(parameter) <- "observed rank"
	method <- paste("Mantel permutation test for", type, "measure")
	data.name <- paste(deparse(substitute(var)), "\nweights:",
	    deparse(substitute(listw)), "\nnumber of simulations + 1:",
	    nsim+1, "\n")
	est <- c(mean(res[1:nsim]), sd(res[1:nsim]))
	names(est) <- c("mean of permutations", "sd of permutations")
	lres <- list(statistic=statistic, parameter=parameter,
	    alternative=alternative, method=method, data.name=data.name, 
	    p.value=pval, res=res, estimate=est)
	class(lres) <- c("htest", "mc.sim")
	lres
}

plot.mc.sim <- function(x, ...) {
	res <- x$res
	xlim <- range(res)
	n <- length(res)
	obs <- res[n]
	res <- res[-n]
	plot(density(res), xlim=xlim, xlab=strsplit(x$data.name, "\n")[[1]][1], 
		main="Density plot of permutation outcomes", sub=x$method)
	abline(v=obs)
}
# Copyright 2003 by Roger Bivand 

set.spChkOption <- function(check) {
	if (!is.logical(check)) stop ("logical argument required")
	res <- get("spChkID", env = .spChkOption)
	assign("spChkID", check, env = .spChkOption)
	res
}

get.spChkOption <- function() {
	get("spChkID", env = .spChkOption)
}

chkIDs <- function (x, listw) 
{
    if (!is.array(x) & !is.data.frame(x)) {
        if (is.null(xn <- names(x))) 
            stop(paste(deparse(substitute(x)), "has no names"))
    }
    else {
        if (is.null(xn <- rownames(x))) 
            stop(paste(deparse(substitute(x)), "has no row names"))
    }
    if (!inherits(listw, "nb")) 
        stop(paste(deparse(substitute(listw)), "is not an listw  or nb object"))
    if (is.null(ln <- attr(listw, "region.id"))) 
        stop(paste(deparse(substitute(listw)), "has no region IDs"))
    if (length(ln) != length(xn)) 
        stop("objects of different length")
    res <- all(ln == xn)
    res
}

spNamedVec <- function(var, data) {
	if (!is.array(data) & !is.data.frame(data))
		stop(paste(deparse(substitute(data)),
			"not an array or data frame"))
	if (!is.character(var) & !is.numeric(var)) 
		stop("variable name wrong type") 
	res <- try(data[,var])
	if (inherits(res, "try-error")) 
		stop(paste(deparse(substitute(var)), "not found"))
	nms <- rownames(data)
	if (is.null(nms)) nms <- as.character(1:length(res))
	names(res) <- nms
	res
}
# Copyright 2000-2 by Roger S. Bivand. 
#

spwdet <- function(sparseweights, rho, debug=FALSE)
{
	if(!inherits(sparseweights, "spatial.neighbour"))
             stop("Not a sparse weights object")
	if(missing(rho) || !is.numeric(rho))
		stop("rho incorrectly specified")
	n <- length(attr(sparseweights, "region.id"))
	size <- length(sparseweights$weights)
	vals <- -rho*sparseweights$weights
	if (debug) verbose <- 1
	else verbose <- 0
	z <- .C("spRdet",
			n = as.integer(n),
			size = as.integer(size),
			verbose = as.integer(verbose),
			p1 = as.integer(sparseweights$from),
			p2 = as.integer(sparseweights$to),
			value = as.double(vals),
			determinant = double(1),
			pideterminant = double(1),
			exponent = integer(1),
			PACKAGE="spdep"
	)
	list(det=z$determinant, exp=z$exponent)
}

logSpwdet <- function(sparseweights, rho, debug=FALSE)
{
	if(!inherits(sparseweights, "spatial.neighbour"))
             stop("Not a sparse weights object")
	if(missing(rho) || !is.numeric(rho) ||
		rho >= 1 || rho <= -1)
		stop("rho incorrectly specified")
	res <- spwdet(sparseweights, rho=rho, debug=debug)
	fres <- log(res$det) + log(10)*res$exp
	fres
}

# Copyright 2001-4 by Roger Bivand
#

subset.nb <- function(x, subset, ...) {
    if (class(x) != "nb") stop("not a neighbours list")
    if (!is.logical(subset)) stop("subset not a logical vector")
    n <- length(x)
    if (n != length(subset))
	stop("neighours list and subset vector different lengths")
    old.ids <- 1:n
    new.ids <- match(old.ids, which(subset))
    reg.id <- subset.default(attr(x, "region.id"), subset)
    x <- sym.attr.nb(x)
    xattrs <- names(attributes(x))
    z <- subset.default(x, subset)
    nz <- length(z)
    for (i in 1:nz) {
	zi <- z[[i]]
	res <- NULL
	for (j in 1:length(zi)) {
	    a <- new.ids[zi[j]]
	    if (!is.na(a)) res <- c(res, a)
	}
	if (is.null(res)) z[[i]] <- 0
	else z[[i]] <- sort(unique(res))
    }
    attr(z, "region.id") <- reg.id
    for (i in 1:length(xattrs)) {
	if (xattrs[i] != "region.id")
	    attr(z, xattrs[i]) <- attr(x, xattrs[i])
    }
    invisible(z)
}


subset.listw <- function(x, subset, zero.policy=FALSE, ...) {
    if (!inherits(x, "listw")) stop("not a weights list")
    if (!is.logical(subset)) stop("subset not a logical vector")
    nb <- x$neighbours
    vlist <- x$weights
    if (is.null(attr(vlist, "binary"))) 
	stop("Not yet able to subset general weights lists")
    style <- x$style
    n <- length(nb)
    if (n != length(subset))
	stop("neighours list and subset vector different lengths")
    subnb <- subset.nb(x=nb, subset=subset)
    sublistw <- nb2listw(neighbours=subnb, glist=NULL, style=style,
	zero.policy=zero.policy)
    invisible(sublistw)
}


# Copyright 2001-3 by Roger Bivand
#


summary.nb <- function(object, coords=NULL, lonlat=FALSE, scale=1, ...) {
    nb <- object
    if (!inherits(nb, "nb")) stop("Not a neighbours list")
    c.nb <- card(nb)
    n.nb <- length(nb)
    regids <- attr(nb, "region.id")
    if(is.null(regids)) regids <- as.character(1:n.nb)
    print.nb(object)
    cat("Link number distribution:\n")
    print(table(c.nb, deparse.level=0))
    if(any(c.nb > 0)) {
        min.nb <- min(c.nb[c.nb > 0])
        cat(length(c.nb[c.nb == min.nb]), " least connected region",
	    ifelse(length(c.nb[c.nb == min.nb]) < 2, "", "s"), ":\n",
	    paste(regids[which(c.nb == min.nb)], collapse=" "), " with ",
	    min.nb, " link", ifelse(min.nb < 2, "", "s"), "\n", sep="")
        max.nb <- max(c.nb)
	cat(length(c.nb[c.nb == max.nb]), " most connected region",
	    ifelse(length(c.nb[c.nb == max.nb]) < 2, "", "s"), ":\n",
	    paste(regids[which(c.nb == max.nb)], collapse=" "), " with ",
	    max.nb, " link", ifelse(max.nb < 2, "", "s"), "\n", sep="")
    }
    if(!is.null(coords)) {
        if (!is.matrix(coords)) stop("Data not in matrix form")
        if (any(is.na(coords))) stop("Data include NAs")
        np <- nrow(coords)
	if(np != n.nb) stop("Number of coords not equal to number of regions")
        dimension <- ncol(coords)
	dlist <- .Call("nbdists", nb, as.matrix(coords), as.integer(np), 
	    as.integer(dimension), as.integer(lonlat), PACKAGE="spdep")[[1]]
	cat("Summary of link distances:\n")
	print(summary(unlist(dlist)))
	stem(unlist(dlist), scale=scale)
    }
}


print.nb <- function(x, ...) {
    nb <- x
    if (!inherits(nb, "nb")) stop("Not a neighbours list")
    c.nb <- card(nb)
    n.nb <- length(nb)
    regids <- attr(nb, "region.id")
    if(is.null(regids)) regids <- as.character(1:n.nb)
    cat("Neighbour list object:\n")
    cat("Number of regions:", n.nb, "\n")
    cat("Number of nonzero links:", sum(c.nb), "\n")
    cat("Percentage nonzero weights:", (100*sum(c.nb))/(n.nb^2), "\n")
    cat("Average number of links:", mean(c.nb), "\n")
    if(any(c.nb == 0)) cat(length(c.nb[c.nb == 0]), " region", 
        ifelse(length(c.nb[c.nb == 0]) < 2, "", "s"), " with no links:\n",
	paste(regids[which(c.nb == 0)], collapse=" "), "\n", sep="")
    res <- is.symmetric.nb(x, verbose=FALSE)
    if (!res) cat("Non-symmetric neighbours list\n")
    invisible(x)
}

summary.listw <- function(object, coords=NULL, lonlat=FALSE, 
	zero.policy=FALSE, scale=1, ...) {
	cat("Characteristics of weights list object:\n")
	summary(object$neighbours, coords=coords, lonlat=lonlat, 
		scale=scale, ...)
	style <- object$style
	cat(paste("\nWeights style:", style, "\n"))
	cat("Weights constants summary:\n")
	print(data.frame(rbind(unlist(spweights.constants(object,
		zero.policy=zero.policy))[c(1, 5:8)]), row.names=style))

}

print.listw <- function(x, zero.policy=FALSE, ...) {
	cat("Characteristics of weights list object:\n")
	print.nb(x$neighbours, ...)
	style <- x$style
	cat(paste("\nWeights style:", style, "\n"))
	cat("Weights constants summary:\n")
	print(data.frame(rbind(unlist(spweights.constants(x,
		zero.policy=zero.policy))[c(1, 5:8)]), row.names=style))
	invisible(x)

}

# Copyright 1998-2004 by Roger Bivand (Wald test suggested by Rein Halbersma,
# output of correlations suggested by Michael Tiefelsdorf)
#

print.sarlm <- function(x, ...)
{
	cat("\nCall:\n")
	print(x$call)
	cat("Type:", x$type, "\n")
	cat("\nCoefficients:\n")
	print(coef(x))
	cat("\nLog likelihood:", logLik(x), "\n")
	invisible(x)
}

summary.sarlm <- function(object, correlation = FALSE, ...)
{
	if (object$type == "error" || ((object$type == "lag" || 
		object$type == "mixed") && object$ase)) {
		object$coeftitle <- "(asymptotic standard errors)"
		object$Coef <- cbind(object$coefficients, object$rest.se, 
			object$coefficients/object$rest.se,
			2*(1-pnorm(abs(object$coefficients/object$rest.se))))
		colnames(object$Coef) <- c("Estimate", "Std. Error", 
			"z value", "Pr(>|z|)")
	} else {
		object$coeftitle <- "(log likelihood/likelihood ratio)"
		m <- length(object$coefficients)
		LLs <- numeric(m)
		LRs <- numeric(m)
		Pvals <- numeric(m)
		if (length(object$LLs) == m) {
			for (i in 1:m) {
				LLs[i] <- c(object$LLs[[i]])
				res <- LR.sarlm(object, object$LLs[[i]])
				LRs[i] <- res$statistic
				Pvals[i] <- res$p.value
			}
		} else {
			LLs[1] <- LRs[1] <- Pvals[1] <- NA
			for (i in 1:(m-1)) {
				LLs[i+1] <- c(object$LLs[[i]])
				res <- LR.sarlm(object, object$LLs[[i]])
				LRs[i+1] <- res$statistic
				Pvals[i+1] <- res$p.value
			}
		}
		object$Coef <- cbind(object$coefficients, LLs, LRs, Pvals)
		colnames(object$Coef) <- c("Estimate", "Log likelihood",
			"LR statistic", "Pr(>|z|)")
	}
	if (object$ase) {
		object$Wald1 <- Wald1.sarlm(object)
		if (correlation) {
			object$correlation <- diag((diag(object$resvar))
				^(-1/2)) %*% object$resvar %*% 
				diag((diag(object$resvar))^(-1/2))
			dimnames(object$correlation) <- dimnames(object$resvar)
		}
	}
	object$LR1 <- LR1.sarlm(object)
	rownames(object$Coef) <- names(object$coefficients)

	structure(object, class=c("summary.sarlm", class(object)))
}

LR1.sarlm <- function(object)
{
	if (!inherits(object, "sarlm")) stop("Not a sarlm object")
	LLx <- logLik(object)
	LLy <- logLik(object$lm.model)
	statistic <- 2*(LLx - LLy)
	attr(statistic, "names") <- "Likelihood ratio"
	parameter <- abs(attr(LLx, "df") - attr(LLy, "df"))
	if (parameter < 1) 
		stop("non-positive degrees of freedom: no test possible")
	attr(parameter, "names") <- "df"
	p.value <- 1 - pchisq(abs(statistic), parameter)
	estimate <- c(LLx, LLy)
	if (object$type == "error") alt <- "spatial error model"
	else alt <- "spatial lag model"
	attr(estimate, "names") <- c(paste("Log likelihood of",
		alt), paste("Log likelihood of OLS fit",
		deparse(substitute(y))))
	method <- "Likelihood Ratio diagnostics for spatial dependence"
	res <- list(statistic=statistic, parameter=parameter,
		p.value=p.value, estimate=estimate, method=method)
	class(res) <- "htest"
	res
}

Wald1.sarlm <- function(object) {
	if (!inherits(object, "sarlm")) stop("Not a sarlm object")
	if (!object$ase) 
		stop("Cannot compute Wald statistic: parameter a.s.e. missing")
	LLx <- logLik(object)
	LLy <- logLik(object$lm.model)
	if (object$type == "lag" || object$type == "mixed") {
		estimate <- object$rho
		statistic <- (object$rho / object$rho.se)^2
	} else {
		estimate <- object$lambda
		statistic <- (object$lambda / object$lambda.se)^2
	}
	attr(statistic, "names") <- "Wald statistic"
	parameter <- abs(attr(LLx, "df") - attr(LLy, "df"))
	if (parameter < 1) 
		stop("non-positive degrees of freedom: no test possible")
	attr(parameter, "names") <- "df"
	p.value <- 1 - pchisq(abs(statistic), parameter)
	method <- "Wald diagnostics for spatial dependence"
	res <- list(statistic=statistic, parameter=parameter,
		p.value=p.value, estimate=estimate, method=method)
	class(res) <- "htest"
	res

}


print.summary.sarlm <- function(x, digits = max(5, .Options$digits - 3),
	signif.stars = FALSE, ...)
{
	cat("\nCall:", deparse(x$call),	sep = "", fill=TRUE)
	cat("\nResiduals:\n")
	resid <- residuals(x)
	nam <- c("Min", "1Q", "Median", "3Q", "Max")
	rq <- if (length(dim(resid)) == 2) 
		structure(apply(t(resid), 1, quantile), dimnames = list(nam, 
			dimnames(resid)[[2]]))
	else structure(quantile(resid), names = nam)
	print(rq, digits = digits, ...)
	cat("\nType:", x$type, "\n")
	if (x$zero.policy) {
		zero.regs <- attr(x, "zero.regs")
		if (!is.null(zero.regs))
			cat("Regions with no neighbours included:\n",
			zero.regs, "\n")
	}
	cat("Coefficients:", x$coeftitle, "\n")
	printCoefmat(x$Coef, signif.stars=signif.stars, digits=digits,
		na.print="")
#	res <- LR.sarlm(x, x$lm.model)
	res <- x$LR1
	if (x$type == "error") {
		cat("\nLambda:", format(signif(x$lambda, digits)),
			"LR test value:", format(signif(res$statistic, digits)),
			"p-value:", format.pval(res$p.value, digits), "\n")
		if (x$ase) {
			cat("Asymptotic standard error:", 
			format(signif(x$lambda.se, digits)),
			"z-value:",format(signif((x$lambda/
				x$lambda.se), digits)),
			"p-value:", format.pval(2*(1-pnorm(abs(x$lambda/
				x$lambda.se))), digits), "\n")
			cat("Wald statistic:", format(signif(x$Wald1$statistic, 
			digits)), "p-value:", format.pval(x$Wald1$p.value, 
			digits), "\n")
		}
	} else {
		cat("\nRho:", format(signif(x$rho, digits)),
			"LR test value:", format(signif(res$statistic, digits)),
			"p-value:", format.pval(res$p.value, digits), "\n")
		if (x$ase) {
			cat("Asymptotic standard error:", 
			format(signif(x$rho.se, digits)),
			"z-value:",format(signif((x$rho/x$rho.se), digits)),
			"p-value:", format.pval(2 * (1 - pnorm(abs(x$rho/
				x$rho.se))), digits), "\n")
			cat("Wald statistic:", format(signif(x$Wald1$statistic, 
			digits)), "p-value:", format.pval(x$Wald1$p.value, 
			digits), "\n")
		}

	}
	cat("\nLog likelihood:", logLik(x), "for", x$type, "model\n")
	cat("ML residual variance (sigma squared): ", 
		format(signif(x$s2, digits)), ", (sigma: ", 
		format(signif(sqrt(x$s2), digits)), ")\n", sep="")
	cat("Number of observations:", length(x$residuals), "\n")
	cat("Number of parameters estimated:", x$parameters, "\n")
	cat("AIC: ", format(signif(AIC(x), digits)), ", (AIC for lm: ",
		format(signif(AIC(x$lm.model), digits)), ")\n", sep="")
	if (x$type != "error" && x$ase) {
		cat("LM test for residual autocorrelation\n")
		cat("test value:", format(signif(x$LMtest, digits)),
			"p-value:", format.pval((1 - pchisq(x$LMtest, 1)), 
			digits), "\n")
	}
    	correl <- x$correlation
    	if (!is.null(correl)) {
        	p <- NCOL(correl)
        	if (p > 1) {
            		cat("\nCorrelation of Coefficients:\n")
                	correl <- format(round(correl, 2), nsmall = 2, 
                  	digits = digits)
                	correl[!lower.tri(correl)] <- ""
                	print(correl[-1, -p, drop = FALSE], quote = FALSE)
            	}
    	}
    	cat("\n")
        invisible(x)
}
# Copyright 2001-4 by Roger Bivand
#


tri2nb <- function(coords, row.names = NULL) {
	require(tripack)
	n <- nrow(coords)
	if (n < 3) stop("too few coordinates")
#	left <- function(x) {
#		res <- (x[3,1]-x[2,1])*(x[1,2]-x[2,2]) >= 
#			(x[1,1]-x[2,1])*(x[3,2]-x[2,2])
#		res
#	}
#	if (left(coords[1:3,])) stop("first three coordinates collinear")
    	if (!is.null(row.names)) {
		if(length(row.names) != n)
            		stop("row.names wrong length")
		if (length(unique(row.names)) != length(row.names))
	    		stop("non-unique row.names given")
    	}
    	if (is.null(row.names)) row.names <- as.character(1:n)
	tri <- tri.mesh(x=coords[,1], y=coords[,2])
	nb <- neighbours(tri)
 	attr(nb, "region.id") <- row.names
	class(nb) <- "nb"
	attr(nb, "tri") <- TRUE
	attr(nb, "call") <- match.call()
	nb <- sym.attr.nb(nb)
	invisible(nb)
}

# Copyright 2001-3 by Roger Bivand 
#

spweights.constants <- function(listw, zero.policy=FALSE) {
	if(!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	cards <- card(listw$neighbours)
	if (!zero.policy && any(cards == 0))
		stop("regions with no neighbours found")
	n <- length(which(cards > 0))
	n1 <- n - 1
	n2 <- n - 2
	n3 <- n - 3
	nn <- n*n
	S0 <- Szero(listw)
	S1 <- 0
	S2 <- 0
	for (i in 1:length(listw$neighbours)) {
		cond <- TRUE
		if (zero.policy && cards[i] == 0) cond <- FALSE
		if (cond) {
			if (cards[i] == 0)
				stop(paste("region", i,
					"has no neighbours"))
			ij <- listw$neighbours[[i]]
			wij <- listw$weights[[i]]
			dm0 <- 0
			dm1 <- 0
			for (j in 1:length(ij)) {
				dij <- wij[j]
				ij.j <- ij[j]
				ij.lkup <- which(listw$neighbours[[ij.j]] == i)
				if (length(ij.lkup) == 1)
					dji <- listw$weights[[ij.j]][ij.lkup]
				else dji <- 0
				dm0 <- dm0 + dij
				dm1 <- dm1 + dji
				S1 <- S1 + (dij + dji)^2
			}
			S2 <- S2 + (dm0 + dm1)^2
		}
	}
	S1 <- S1 * 0.5
	invisible(list(n=n, n1=n1, n2=n2, n3=n3, nn=nn, S0=S0, S1=S1, S2=S2))
}

Szero <- function(listw) {
	sum(unlist(listw$weights))
}

lag.listw <- function(x, var, zero.policy=FALSE, NAOK=FALSE, ...) {
	listw <- x
	x <- var
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	if (!is.numeric(x)) stop(paste(deparse(substitute(listw)),
		"not numeric"))
	if (!is.logical(NAOK)) stop("NAOK must be logical")
#	if (any(is.na(x))) stop("NA in X")
	n <- length(listw$neighbours)
	cardnb <- card(listw$neighbours)
	if (is.null(dim(x))) {
		if (length(x) != n) stop("object lengths differ")
		res <- .Call("lagw", listw$neighbours, listw$weights,
			as.double(x), as.integer(cardnb),
			as.logical(zero.policy), NAOK=NAOK, PACKAGE="spdep")
	} else {
		if (nrow(x) != n) stop("object lengths differ")
		res <- matrix(0, nrow=nrow(x), ncol=ncol(x))
		for (i in 1:ncol(x)) {
			res[,i] <- .Call("lagw", listw$neighbours,
				listw$weights, as.double(x[,i]),
				as.integer(cardnb), as.logical(zero.policy),
				NAOK=NAOK, PACKAGE="spdep")

		}
	} 
	if (any(is.na(res))) warning("NAs in lagged values")
	invisible(res)
}

listw2U <- function(listw) {
	if (!inherits(listw, "listw")) stop(paste(deparse(substitute(listw)),
		"is not a listw object"))
	nb <- listw$neighbours
	wts <- listw$weights
	style <- paste(listw$style, "U", sep="")
	sym <- is.symmetric.nb(nb, FALSE, TRUE)
	n <- length(listw$neighbours)
	cardnb <- card(listw$neighbours)
	nlist <- vector(mode="list", length=n)
	attr(nlist, "region.id") <- attr(nb, "region.id")
	class(nlist) <- "nb"
	vlist <- vector(mode="list", length=n)
	attr(vlist, as.character(style)) <- TRUE
	if (sym) {
		nlist <- vector(mode="list", length=n)
		attr(nlist, "region.id") <- attr(nb, "region.id")
		class(nlist) <- "nb"
		for (i in 1:n) {
			inb <- nb[[i]]
			nlist[[i]] <- inb
			iwt <- wts[[i]]
			icd <- cardnb[i]
			if (icd > 0) {
			    for (j in 1:icd) {
				vlist[[i]][j] <- 0.5 *
				(iwt[j]+wts[[inb[j]]][which(nb[[inb[j]]] == i)])
			    }
			}
		}
	} else {
		nlist <- make.sym.nb(nb)
		for (i in 1:n) {
			inb <- nb[[i]]
			inl <- nlist[[i]]
			if (inl[1] > 0) {
			    iwt <- wts[[i]]
			    vlist[[i]] <- numeric(length=length(inl))
			    for (j in 1:length(inl)) {
				if (inl[j] %in% inb) 
				    a <- iwt[which(inb == inl[j])]
				else a <- 0
				if (i %in% nb[[inl[j]]]) 
				    b <- wts[[inl[j]]][which(nb[[inl[j]]] == i)]
				else b <- 0
				vlist[[i]][j] <- 0.5 * (a + b)
			    }
			}
		}
	}
	res <- list(style=style, neighbours=nlist, weights=vlist)
	class(res) <- "listw"
	attr(res, "region.id") <- attr(nb, "region.id")
	attr(res, "call") <- match.call()
	attr(res, "U") <- TRUE
	invisible(res)
}


listw2star <- function(listw, ireg, style, n, D, a, zero.policy=FALSE) {
    nb <- vector(mode="list", length=n)
    class(nb) <- "nb"
    wts <- vector(mode="list", length=n)
    for (i in 1:n) nb[[i]] <- 0
    inb <- listw$neighbours[[ireg]]
    iwts <- listw$weights[[ireg]]
    cond <- TRUE
    if (inb == 0 || length(inb) == 0 || is.null(iwts)) cond <- FALSE
    if (!cond && !zero.policy) stop("No-neighbour region found")
    if (style == "W") iwts <- (n*D[ireg]*iwts) / 2
    else if (style == "S") iwts <- ((n^2)*D[ireg]*iwts) / (2*a)
    else if (style == "C") iwts <- ((n^2)*iwts) / (2*a)
    if (cond) {
    	nb[[ireg]] <- inb
    	wts[[ireg]] <- iwts
    	for (j in 1:length(inb)) {
            jj <- inb[j]
            nb[[jj]] <- ireg
            wts[[jj]] <- iwts[j]
	}
    }
    res <- list(style=style, neighbours=nb, weights=wts)
    class(res) <- c("listw", "star")
    attr(res, "region.id") <- attr(listw, "region.id")
    res
}

spdep <- function(build=FALSE) {
	if (build) return(c(.spdep.Version, .spdep.Build))
	else return(.spdep.Version)
}
# Copyright 2001-4 by Roger Bivand 
#


is.symmetric.nb <- function(nb, verbose=TRUE, force=FALSE)
{
	if(!inherits(nb, "nb")) stop("Not neighbours list")
	nbsym <- attr(nb, "sym")
	if(!is.null(nbsym)) res <- nbsym
	if(force || is.null(nbsym)) {
		res <- .Call("symtest", nb=nb, card=as.integer(card(nb)),
			verbose=as.logical(verbose), PACKAGE="spdep")
	}
	if(!res && verbose) cat("Non-symmetric neighbours list\n")
	invisible(res)
}

sym.attr.nb <- function(nb) {
	if(!inherits(nb, "nb")) stop("Not neighbours list")
	nbsym <- attr(nb, "sym")
	if(is.null(nbsym))
		attr(nb, "sym") <- is.symmetric.nb(nb, verbose=FALSE,
			force=TRUE)
	invisible(nb)
}

include.self <- function(nb) {
	if (!is.null(attributes(nb)$self.included) &&
		(as.logical(attributes(nb)$self.included)))
		stop("Self already included")
	n <- length(nb)
	nc <- card(nb)
	for (i in 1:n) {
		if (nc[i] > 0) {
			nb[[i]] <- sort(c(i, nb[[i]]))
		} else {
			nb[[i]] <- i
		}
	}
		
	attr(nb, "self.included") <- TRUE
	invisible(nb)
}

# Copyright 2001 by Nicholas Lewin-Koh 

make.sym.nb <- function(nb){
	if(!inherits(nb, "nb")) stop("Not neighbours list")
	if (is.symmetric.nb(nb, FALSE, TRUE)) {
		res <- nb
	} else {
        	k <- unlist(lapply(nb,length))
        	to <- unlist(nb)
        	from <- NULL
        	res <- vector(mode="list", length=length(nb))
        	for(i in 1:length(nb)){
        		from <- c(from,rep(i,k[i]))
        	}
        	for(i in 1:length(nb)){
        		res[[i]] <- sort(unique(c(to[from==i],from[to==i])))
        		if(length(res[[i]]) == 0) res[[i]] <- 0
        	}
        	attr(res, "region.id") <- attr(nb,"region.id")
        	attr(res, "call") <- attr(nb, "call")
        	attr(res, "type") <- attr(nb, "type")
        	attr(res, "sym") <- TRUE
        	class(res) <- "nb"
	}
	invisible(res)
}

# Copyright 2001-3 by Roger Bivand 
#

.spChkOption <- new.env(FALSE, globalenv())
assign("spChkID", FALSE, env = .spChkOption)
.DESC <- packageDescription("spdep")
.spdep.Version <- paste(.DESC[["Package"]], ", version ", .DESC[["Version"]],
	 ", ", .DESC[["Date"]], sep="")
.spdep.Build <- paste("build:", .DESC[["Built"]])
#.onLoad <- function(pkg, lib) {
cat(paste(spdep()[1], ":\n a package for analysing spatial dependence,\n",
	" use help(get.spChkOption) for help on integrity checking\n",
	sep=""))
#require(maptools)
#.First.lib <- function(lib, pkg) {
#	library.dynam("spdep", pkg, lib)
#}
.noGenerics <- TRUE

