[r-cran-spatstat.utils] 01/02: New upstream version 1.7-0
Andreas Tille
tille at debian.org
Thu Sep 7 11:36:04 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-spatstat.utils.
commit 19757457854a6aafb1e124c2b10962253a7e1de7
Author: Andreas Tille <tille at debian.org>
Date: Thu Sep 7 13:34:52 2017 +0200
New upstream version 1.7-0
---
DESCRIPTION | 30 ++
MD5 | 60 +++
NAMESPACE | 192 +++++++++
R/formulae.R | 242 +++++++++++
R/locator.R | 36 ++
R/primefactors.R | 201 +++++++++
R/resolve.defaults.R | 178 ++++++++
R/tapplysum.R | 104 +++++
R/utilarg.R | 385 +++++++++++++++++
R/utilindex.R | 100 +++++
R/utilmatrix.R | 182 ++++++++
R/utiloptim.R | 34 ++
R/utilseq.R | 406 ++++++++++++++++++
R/utiltext.R | 508 ++++++++++++++++++++++
R/xycircle.R | 137 ++++++
R/xypolygon.R | 446 ++++++++++++++++++++
R/xysegment.R | 232 ++++++++++
inst/doc/packagesizes.txt | 4 +
man/articlebeforenumber.Rd | 39 ++
man/cat.factor.Rd | 44 ++
man/check.1.integer.Rd | 59 +++
man/check.named.vector.Rd | 96 +++++
man/check.nmatrix.Rd | 79 ++++
man/check.nvector.Rd | 78 ++++
man/check.range.Rd | 97 +++++
man/commasep.Rd | 38 ++
man/do.call.matched.Rd | 83 ++++
man/do.call.without.Rd | 47 +++
man/expand.polynom.Rd | 57 +++
man/ifelseAB.Rd | 76 ++++
man/macros/defns.Rd | 19 +
man/methods.xypolygon.Rd | 127 ++++++
man/optimizeWithTrace.Rd | 65 +++
man/ordinal.Rd | 43 ++
man/paren.Rd | 45 ++
man/primefactors.Rd | 87 ++++
man/resolve.defaults.Rd | 87 ++++
man/revcumsum.Rd | 48 +++
man/simplenumber.Rd | 48 +++
man/spatstat.utils-internal.Rd | 212 ++++++++++
man/spatstat.utils-package.Rd | 121 ++++++
man/spatstatLocator.Rd | 49 +++
man/splat.Rd | 62 +++
man/tapplysum.Rd | 60 +++
man/termsinformula.Rd | 102 +++++
man/verbalogic.Rd | 67 +++
src/chunkloop.h | 37 ++
src/circxseg.c | 930 +++++++++++++++++++++++++++++++++++++++++
src/distseg.c | 175 ++++++++
src/fastinterv.c | 37 ++
src/init.c | 47 +++
src/inxyp.c | 75 ++++
src/matchindices.c | 197 +++++++++
src/matchpoints.c | 52 +++
src/ply.c | 31 ++
src/ply.h | 85 ++++
src/primefax.c | 175 ++++++++
src/proto.h | 37 ++
src/revcum.c | 41 ++
tests/numerical.R | 25 ++
tests/segments.R | 13 +
61 files changed, 7469 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..633e7cc
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,30 @@
+Package: spatstat.utils
+Version: 1.7-0
+Date: 2017-07-27
+Title: Utility Functions for 'spatstat'
+Authors at R: c(person("Adrian", "Baddeley",
+ role = c("aut", "cre"),
+ email = "Adrian.Baddeley at curtin.edu.au"),
+ person("Rolf", "Turner",
+ role = "aut",
+ email="r.turner at auckland.az.nz"),
+ person("Ege", "Rubak",
+ role = "aut",
+ email = "rubak at math.aau.dk"))
+Maintainer: Adrian Baddeley <Adrian.Baddeley at curtin.edu.au>
+Depends: R (>= 3.3.0), stats, graphics, grDevices, utils, methods
+Suggests: spatstat
+Description: Contains utility functions for the 'spatstat' package
+ which may also be useful for other purposes.
+License: GPL (>= 2)
+URL: http://www.spatstat.org
+LazyData: true
+NeedsCompilation: yes
+ByteCompile: true
+BugReports: https://github.com/spatstat/spatstat.utils/issues
+Packaged: 2017-07-27 02:47:01 UTC; adrian
+Author: Adrian Baddeley [aut, cre],
+ Rolf Turner [aut],
+ Ege Rubak [aut]
+Repository: CRAN
+Date/Publication: 2017-08-09 04:20:25 UTC
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..6d3d6c3
--- /dev/null
+++ b/MD5
@@ -0,0 +1,60 @@
+dc76981f156fe6e10a46af4def640de0 *DESCRIPTION
+56c024709d12aa75f0f6057b5352b033 *NAMESPACE
+81a0013c01515a69d9c2345e2aa799c9 *R/formulae.R
+45a76e7b281f65e8d3ea2ebe0bf29c17 *R/locator.R
+fc4f08d1f8aca7ec08bca87a6293f977 *R/primefactors.R
+f8ef809caaa54f24863ca5c670e746e3 *R/resolve.defaults.R
+5442cdc8912d6e43991f00da76ea1078 *R/tapplysum.R
+dee0e2fd5e7ffd4c26e9df03c240f059 *R/utilarg.R
+fa9f66757b4edab07677672cd52e92b9 *R/utilindex.R
+e167b5360b4d4b0e569ad16bb0ee5955 *R/utilmatrix.R
+61dfab95b18e2c29dbe25ab48cf6f4e6 *R/utiloptim.R
+ef8af43f877a3d48c4ac044becda855f *R/utilseq.R
+3d57f5cead75a89fa4b55cb896e2b68a *R/utiltext.R
+3715ef2a8ce85d8caf1c886616e8882a *R/xycircle.R
+f3297f80813c7b1694bafead35972012 *R/xypolygon.R
+c178f8dd4f7476186c33de119f4ffb87 *R/xysegment.R
+3dc3ea12d3871d7f48f4170e860d3b39 *inst/doc/packagesizes.txt
+0c11e69f495fa53d182d3094d84ec3e2 *man/articlebeforenumber.Rd
+1c2612c451d7c808eb03c574f17cb36a *man/cat.factor.Rd
+743ff83b879b31b961d86fda75891cf7 *man/check.1.integer.Rd
+45e23ead6342e7dbbae73ee6955d5eac *man/check.named.vector.Rd
+150d9ed940337f5a0459e22e086ef5a4 *man/check.nmatrix.Rd
+55bd824bcd6953c9525d1bc1ffac4918 *man/check.nvector.Rd
+0c1ca20a2d58ef5fed14476580159e17 *man/check.range.Rd
+39c227b91e4fd4e2e551872f8badb6c2 *man/commasep.Rd
+ed75ba27cff27d36ba6d5aa20595acca *man/do.call.matched.Rd
+6cefb985857ea0c2f8a9b6de8eb20189 *man/do.call.without.Rd
+49dd46740ae040193c0e92c9e6d62519 *man/expand.polynom.Rd
+e410bcae4113c62d2e78d2e9b0fc9c7b *man/ifelseAB.Rd
+8ff61b15a83accee7cc5be6b9363b347 *man/macros/defns.Rd
+54d77c8281fa5b072b081baa4b75d5de *man/methods.xypolygon.Rd
+3d0f1e38d2215913c33112a50c0aa225 *man/optimizeWithTrace.Rd
+30afc78c74fd7b3dffe4d83b61e3174b *man/ordinal.Rd
+ab2589bfd6de3acb12001c1c85b6aa55 *man/paren.Rd
+cb121d5caa9f123a43d5a238c816e322 *man/primefactors.Rd
+912a6809ce279e63cef321a62a7453c7 *man/resolve.defaults.Rd
+777dc7f9df370468156fd7c9deea2d6a *man/revcumsum.Rd
+6a72c33940fa00538b191d17f37172b0 *man/simplenumber.Rd
+968b84c1c593c5f13a6e2a4681f63805 *man/spatstat.utils-internal.Rd
+7d9c345f637717be85ab5b22652fada5 *man/spatstat.utils-package.Rd
+57a5190c70cfcd77c8e154e21cac0223 *man/spatstatLocator.Rd
+ab7258a3b52edb9dc0be7166d511d812 *man/splat.Rd
+f0bcbfddef2871327ef318e27dcb9ae8 *man/tapplysum.Rd
+9f37f779d074314b5943a7f53879ed78 *man/termsinformula.Rd
+e348838985fa1a05b03e652a8bd6a4e6 *man/verbalogic.Rd
+3c6ee73e77c45154372fe8dbb893902f *src/chunkloop.h
+6b8422024408b9a11522b02502be9220 *src/circxseg.c
+2214769b031225151f12603e0f37fd0c *src/distseg.c
+e8baebe0e736ff943f75f68eb7cde4c7 *src/fastinterv.c
+17d0bdbbbe7ca241c9bd341f4e5c1e58 *src/init.c
+f8e6302140d1ee2c3e2185a32b2bddac *src/inxyp.c
+a6543c427067f6d163c65e6132c0ca43 *src/matchindices.c
+5b0cb4b7dde82600b7781671d068bd6c *src/matchpoints.c
+12c672c622be0ab2eed99c9e86de83b6 *src/ply.c
+62033117b1f5209472467cf9a80cd74c *src/ply.h
+c7cb679d7a7d2e5fe33b2a66a129ee09 *src/primefax.c
+178f74ee02d5a6fb89aaa846785e6f36 *src/proto.h
+3597e6c576823a63c3bd94e4c0b3eef9 *src/revcum.c
+b4fca8cfc0c31d47078dd130ca8323c9 *tests/numerical.R
+91e31de1f0a6229db0239e6c6f373582 *tests/segments.R
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..2f59f3c
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,192 @@
+# spatstat.utils NAMESPACE file
+
+import(stats,graphics,grDevices,utils,methods)
+
+# Do not edit the following.
+# It is generated automatically.
+
+
+
+# ..................................................
+# load dynamic library
+# (native routines are now registered in init.c)
+# ..................................................
+
+useDynLib(spatstat.utils, .registration=TRUE, .fixes="C_")
+
+# ..................................................
+# Automatically-generated list of documented objects
+# ..................................................
+export("adjustthinrange")
+export("apply23sum")
+export("Area.xypolygon")
+export("articlebeforenumber")
+export("as2vector")
+export("asNumericMatrix")
+export("assignDFcolumn")
+export("badprobability")
+export("bdrylength.xypolygon")
+export("blockdiagarray")
+export("blockdiagmatrix")
+export("can.be.formula")
+export("cat.factor")
+export("check.1.integer")
+export("check.1.real")
+export("check.1.string")
+export("check.finite")
+export("check.in.range")
+export("check.named.list")
+export("check.named.thing")
+export("check.named.vector")
+export("check.nmatrix")
+export("check.nvector")
+export("check.range")
+export("choptext")
+export("choptextline")
+export("commasep")
+export("complaining")
+export("distpl")
+export("distppl")
+export("distppll")
+export("distppllmin")
+export("divisors")
+export("do.call.matched")
+export("do.call.without")
+export("dont.complain.about")
+export("dotexpr.to.call")
+export("dropifsingle")
+export("dround")
+export("ensure2vector")
+export("ensure3Darray")
+export("equispaced")
+export("eratosthenes")
+export("evenly.spaced")
+export("exceedsMaxArraySize")
+export("exhibitStringList")
+export("expand.polynom")
+export("explain.ifnot")
+export("fakecallstring")
+export("fastFindInterval")
+export("fave.order")
+export("fillseq")
+export("findfirstfactor")
+export("firstfactor")
+export("flat.deparse")
+export("fontify")
+export("forbidNA")
+export("geomseq")
+export("getdataobjects")
+export("good.names")
+export("graphicsPars")
+export("greatest.common.divisor")
+export("gsubdot")
+export("identical.formulae")
+export("ifelse0NA")
+export("ifelse1NA")
+export("ifelseAB")
+export("ifelseAX")
+export("ifelseNegPos")
+export("ifelseXB")
+export("ifelseXY")
+export("indexCartesian")
+export("inject.expr")
+export("insertinlist")
+export("inside.range")
+export("inside.triangle")
+export("inside.xypolygon")
+export("intersect.ranges")
+export("is.blank")
+export("is.hole.xypolygon")
+export("is.parseable")
+export("is.prime")
+export("least.common.multiple")
+export("lhs.of.formula<-")
+export("lhs.of.formula")
+export("lty2char")
+export("make.parseable")
+export("mapstrings")
+export("matchIntegerDataFrames")
+export("matchNameOrPosition")
+export("matcolall")
+export("matcolany")
+export("matcolsum")
+export("matrixsample")
+export("matrowall")
+export("matrowany")
+export("matrowsum")
+export("natozero")
+export("niceround")
+export("NNdist2segments")
+export("numalign")
+export("nzpaste")
+export("offsetsinformula")
+export("optimizeWithTrace")
+export("orderstats")
+export("orderwhich")
+export("ordinal")
+export("ordinalsuffix")
+export("overlap.trapezium")
+export("overlap.xypolygon")
+export("padtowidth")
+export("paren")
+export("passthrough")
+export("paste.expr")
+export("pasteFormula")
+export("pasteN")
+export("prange")
+export("prettydiscrete")
+export("prettyinside")
+export("primefactors")
+export("primesbelow")
+export("prolongseq")
+export("ratiotweak")
+export("relatively.prime")
+export("resolve.1.default")
+export("resolve.defaults")
+export("revcumsum")
+export("reverse.xypolygon")
+export("rhs.of.formula<-")
+export("rhs.of.formula")
+export("romansort")
+export("samefunction")
+export("sensiblevarname")
+export("short.deparse")
+export("simplenumber")
+export("simplify.xypolygon")
+export("singlestring")
+export("spatstatLocator")
+export("splat")
+export("startinrange")
+export("strsplitretain")
+export("substringcount")
+export("sympoly")
+export("tapplysum")
+export("termsinformula")
+export("there.is.no.try")
+export("trap.extra.arguments")
+export("truncline")
+export("unparen")
+export("uptrimat")
+export("validposint")
+export("variablesinformula")
+export("variablesintext")
+export("verbalogic")
+export("verify.xypolygon")
+export("warn.ignored.args")
+export("xysegMcircle")
+export("xysegPcircle")
+export("xysegXcircle")
+
+# ....... Special cases ...........
+export("%orifnull%")
+# ....... End of special cases ...
+
+# .........................................
+# Automatically generated list of S3 methods
+# .........................................
+# .........................................
+# Assignment methods
+# .........................................
+# .........................................
+# End of methods
+# .........................................
diff --git a/R/formulae.R b/R/formulae.R
new file mode 100755
index 0000000..26e819e
--- /dev/null
+++ b/R/formulae.R
@@ -0,0 +1,242 @@
+#
+#
+# formulae.R
+#
+# THIS FILE IS NOW PART OF spatstat.utils
+#
+# Functions for manipulating model formulae
+#
+# $Revision: 1.27 $ $Date: 2017/07/27 02:46:37 $
+#
+# identical.formulae()
+# Test whether two formulae are identical
+#
+# termsinformula()
+# Extract the terms from a formula
+#
+# sympoly()
+# Create a symbolic polynomial formula
+#
+# polynom()
+# Analogue of poly() but without dynamic orthonormalisation
+#
+# -------------------------------------------------------------------
+#
+
+identical.formulae <- function(x, y) {
+ # workaround for bug in all.equal.formula in R 2.5.0
+ if(is.null(y) && !is.null(x))
+ return(FALSE)
+ return(identical(all.equal(x,y), TRUE))
+}
+
+termsinformula <- function(x) {
+ if(is.null(x)) return(character(0))
+ if(class(x) != "formula")
+ stop("argument is not a formula")
+ attr(terms(x), "term.labels")
+}
+
+variablesinformula <- function(x) {
+ if(is.null(x)) return(character(0))
+ if(class(x) != "formula")
+ stop("argument is not a formula")
+ all.vars(as.expression(x))
+}
+
+offsetsinformula <- function(x) {
+ if(is.null(x)) return(character(0))
+ if(class(x) != "formula")
+ stop("argument is not a formula")
+ tums <- terms(x)
+ offs <- attr(tums, "offset")
+ if(length(offs) == 0) return(character(0))
+ vars <- attr(tums, "variables")
+ termnames <- unlist(lapply(vars, deparse))[-1L]
+ termnames[offs]
+}
+
+lhs.of.formula <- function(x) {
+ if(!inherits(x, "formula"))
+ stop("x must be a formula")
+ if(length(as.list(x)) == 3) {
+ # formula has a response: return it
+ return(x[[2L]])
+ }
+ return(NULL)
+}
+
+rhs.of.formula <- function(x, tilde=TRUE) {
+ if(!inherits(x, "formula"))
+ stop("x must be a formula")
+ if(length(as.list(x)) == 3) {
+ # formula has a response: strip it
+ x <- x[-2L]
+ }
+ if(!tilde) # remove the "~"
+ x <- x[[2L]]
+ return(x)
+}
+
+#' assignment operators
+
+"lhs.of.formula<-" <- function (x, value)
+{
+ if (!inherits(x, "formula"))
+ stop("x must be a formula")
+ if (length(as.list(x)) == 2)
+ x[[3L]] <- x[[2L]]
+ x[[2L]] <- value
+ return(x)
+}
+
+"rhs.of.formula<-" <- function (x, value)
+{
+ if (!inherits(x, "formula"))
+ stop("x must be a formula")
+ x[[3L]] <- value
+ return(x)
+}
+
+sympoly <- function(x,y,n) {
+
+ if(nargs()<2) stop("Degree must be supplied.")
+ if(nargs()==2) n <- y
+ eps <- abs(n%%1)
+ if(eps > 0.000001 | n <= 0) stop("Degree must be a positive integer")
+
+ x <- deparse(substitute(x))
+ temp <- NULL
+ left <- "I("
+ rght <- ")"
+ if(nargs()==2) {
+ for(i in 1:n) {
+ xhat <- if(i==1) "" else paste("^",i,sep="")
+ temp <- c(temp,paste(left,x,xhat,rght,sep=""))
+ }
+ }
+ else {
+ y <- deparse(substitute(y))
+ for(i in 1:n) {
+ for(j in 0:i) {
+ k <- i-j
+ xhat <- if(k<=1) "" else paste("^",k,sep="")
+ yhat <- if(j<=1) "" else paste("^",j,sep="")
+ xbit <- if(k>0) x else ""
+ ybit <- if(j>0) y else ""
+ star <- if(j*k>0) "*" else ""
+ term <- paste(left,xbit,xhat,star,ybit,yhat,rght,sep="")
+ temp <- c(temp,term)
+ }
+ }
+ }
+ as.formula(paste("~",paste(temp,collapse="+")))
+ }
+
+
+
+expand.polynom <- local({
+
+ Iprefix <- function(x) { paste0("I", x) }
+ Iparen <- function(x) { Iprefix(paren(x)) }
+
+ powername <- function(x, n) {
+ ifelse(n == 0, "",
+ ifelse(n == 1,
+ x,
+ paste0(x, "^", n)))
+ }
+
+ power1name <- function(x, n, xisname) {
+ z <- powername(x, n)
+ z[n > 1] <- Iparen(z[n > 1])
+ if(!xisname)
+ z[n == 1] <- Iprefix(z[n == 1])
+ return(z)
+ }
+
+ power2name <- function(x, y, n, m, xisname, yisname) {
+ ifelse(n == 0,
+ power1name(y, m, yisname),
+ ifelse(m == 0,
+ power1name(x, n, xisname),
+ Iparen(paste(powername(x, n),
+ powername(y, m), sep="*"))))
+ }
+
+ haspolynom <- function(z) { 'polynom' %in% all.names(z) }
+
+ fiddle <- function(f) {
+ opname <- f[[1L]]
+ if(identical(opname, as.name('I'))) {
+ ## expressions enclosed in I() are protected
+ return(f)
+ }
+ if(!identical(opname, as.name('polynom'))) {
+ tbd <- unlist(lapply(f, haspolynom))
+ if(any(tbd)) {
+ ## descend recursively
+ for(i in which(tbd))
+ f[[i]] <- fiddle(f[[i]])
+ }
+ return(f)
+ }
+ ## polynom(..., d)
+ n <- length(f)
+ if(!(n %in% c(3,4)))
+ stop("Syntax of polynom() call not understood")
+ degree <- f[[n]]
+ if (!is.numeric(degree) || length(degree) != 1 ||
+ (degree%%1) != 0 || degree < 1)
+ stop("degree of polynomial should be a positive integer")
+ if(n == 3) {
+ ## polynom(x, d)
+ xlang <- f[[2L]]
+ xisname <- (length(xlang) == 1)
+ xstring <- if(xisname) paste(xlang) else paren(format(xlang))
+ xpowers <- power1name(xstring, 1:degree, xisname)
+ xpolystring <- paste(xpowers, collapse=" + ")
+ xpolylang <- as.formula(paste("~", xpolystring))[[2L]]
+ return(xpolylang)
+ } else if(n == 4) {
+ ## polynom(x, y, d)
+ xlang <- f[[2L]]
+ ylang <- f[[3L]]
+ xisname <- (length(xlang) == 1)
+ yisname <- (length(ylang) == 1)
+ xstring <- if(xisname) paste(xlang) else paren(format(xlang))
+ ystring <- if(yisname) paste(ylang) else paren(format(ylang))
+ mat <- matrix(, 1+degree, 1+degree)
+ totdeg <- col(mat) - 1
+ yd <- row(mat) - 1
+ xd <- totdeg - yd
+ xdeg <- xd[xd >= 0]
+ ydeg <- yd[xd >= 0]
+ xypowers <- power2name(xstring, ystring, xdeg, ydeg,
+ xisname, yisname)[xdeg + ydeg > 0]
+ xypolystring <- paste(xypowers, collapse=" + ")
+ xypolylang <- as.formula(paste("~", xypolystring))[[2L]]
+ return(xypolylang)
+ }
+ }
+
+ expand.polynom <- function(f) {
+ ## replaces polynom(...) by x + I(x^2) + ... inside a formula f
+ g <- fiddle(f)
+ environment(g) <- environment(f)
+ return(g)
+ }
+
+ expand.polynom
+})
+
+can.be.formula <- function(x) {
+ #' test whether x is a formula object
+ if(inherits(x, "formula")) return(TRUE)
+ #' or a character representation of a formula.
+ if(!is.character(x)) return(FALSE)
+ x <- paste(x, collapse=" ")
+ if(length(grep("~", x)) == 0) return(FALSE)
+ ok <- !inherits(try(as.formula(x), silent=TRUE), "try-error")
+ return(ok)
+}
diff --git a/R/locator.R b/R/locator.R
new file mode 100644
index 0000000..11f552b
--- /dev/null
+++ b/R/locator.R
@@ -0,0 +1,36 @@
+#'
+#' locator.R
+#'
+#' $Revision: 1.1 $ $Date: 2017/01/07 09:23:51 $
+
+spatstatLocator <- function(n, type=c("p","l","o","n"), ...) {
+ #' remedy for failure of locator(type="p") in RStudio
+ if(!identical(TRUE, dev.capabilities()$locator))
+ stop("Sorry, this graphics device does not support the locator() function")
+ # validate
+ type <- match.arg(type)
+ do.points <- type %in% c("p","o")
+ do.lines <- type %in% c("l","o")
+ argh <- list(...)
+ pointsArgs <- c("cex", "col", "pch", "fg", "bg")
+ segmentArgs <- graphicsPars("lines")
+ # go
+ res <- list(x=numeric(0), y = numeric(0))
+ i <- 1L
+ if(missing(n)) n <- Inf
+ while(i<=n){
+ tmp <- locator(n=1)
+ if(is.null(tmp)) return(res)
+ if(do.points)
+ do.call.matched(points.default, append(tmp, argh), extrargs=pointsArgs)
+ res$x <- c(res$x,tmp$x)
+ res$y <- c(res$y,tmp$y)
+ if(do.lines && i > 1L) {
+ xy <- with(res, list(x0=x[i-1L], y0=y[i-1L], x1=x[i], y1=y[i]))
+ do.call.matched(segments, append(xy, argh), extrargs=segmentArgs)
+ }
+ i <- i+1L
+ }
+ return(res)
+}
+
diff --git a/R/primefactors.R b/R/primefactors.R
new file mode 100755
index 0000000..e1b94a2
--- /dev/null
+++ b/R/primefactors.R
@@ -0,0 +1,201 @@
+#
+# primefactors.R
+#
+# $Revision: 1.8 $ $Date: 2016/12/31 08:58:36 $
+#
+
+# all primes below 2^13=8192
+
+primestable <-
+ c(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
+ 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
+ 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191,
+ 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
+ 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
+ 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
+ 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
+ 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
+ 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
+ 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757,
+ 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853,
+ 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
+ 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021,
+ 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,
+ 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
+ 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259,
+ 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
+ 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433,
+ 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
+ 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
+ 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
+ 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
+ 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831,
+ 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913,
+ 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
+ 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087,
+ 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
+ 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269,
+ 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347,
+ 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,
+ 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
+ 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
+ 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
+ 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767,
+ 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
+ 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953,
+ 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041,
+ 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163,
+ 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251,
+ 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329,
+ 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
+ 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517,
+ 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
+ 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673,
+ 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767,
+ 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853,
+ 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931,
+ 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027,
+ 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
+ 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229,
+ 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
+ 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421,
+ 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513,
+ 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603,
+ 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691,
+ 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793,
+ 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,
+ 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987,
+ 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077,
+ 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171,
+ 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
+ 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393,
+ 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,
+ 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557,
+ 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,
+ 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741,
+ 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839,
+ 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903,
+ 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043,
+ 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131,
+ 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
+ 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311,
+ 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,
+ 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521,
+ 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607,
+ 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703,
+ 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803,
+ 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899,
+ 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983,
+ 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079,
+ 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
+ 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307,
+ 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433,
+ 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523,
+ 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589,
+ 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687,
+ 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
+ 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883,
+ 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009,
+ 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101,
+ 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191)
+
+primestableMax <- 8192
+
+primesbelow <- function(nmax) {
+ if(nmax <= primestableMax) return(primestable[primestable <= nmax])
+ eratosthenes(nmax, c(primestable, primestableMax:nmax))
+}
+
+eratosthenes <- function(nmax, startset=2:nmax) {
+ # The Sieve of Eratosthenes
+ if(nmax < 2) return(numeric(0))
+ numbers <- startset
+ prime <- startset[1]
+ repeat{
+ retain <- (numbers <= prime) | (numbers %% prime != 0)
+ numbers <- numbers[retain]
+ remaining <- (numbers > prime)
+ if(!any(remaining))
+ break
+ prime <- min(numbers[remaining])
+ }
+ return(numbers)
+}
+
+primefactors <- function(n, method=c("C", "interpreted")) {
+ method <- match.arg(method)
+ switch(method,
+ interpreted = {
+ prmax <- floor(sqrt(n))
+ result <- findprimefactors(n, primesbelow(prmax))
+ },
+ C = {
+ check.1.integer(n)
+ kmax <- ceiling(log2(n))
+ z <- .C(C_primefax,
+ n=as.integer(n),
+ factors=as.integer(integer(kmax)),
+ nfactors=as.integer(integer(1L)),
+ PACKAGE = "spatstat.utils")
+ result <- z$factors[seq_len(z$nfactors)]
+ },
+ stop("Unrecognised method"))
+ return(result)
+}
+
+findprimefactors <- function(n, primes) {
+ divides.n <- (n %% primes == 0)
+ if(!any(divides.n))
+ return(n)
+ divisors <- primes[divides.n]
+ m <- n/prod(divisors)
+ if(m == 1) return(divisors)
+ mfactors <- findprimefactors(m, divisors)
+ return(sort(c(divisors, mfactors)))
+}
+
+is.prime <- function(n) { length(primefactors(n)) == 1 }
+
+relatively.prime <- function(n, m) {
+ cf <- intersect(primefactors(n), primefactors(m))
+ return(length(cf) == 0)
+}
+
+least.common.multiple <- function(n, m) {
+ nf <- primefactors(n)
+ mf <- primefactors(m)
+ p <- sort(unique(c(nf,mf)))
+ nfac <- table(factor(nf, levels=p))
+ mfac <- table(factor(mf, levels=p))
+ prod(p^pmax.int(nfac,mfac))
+}
+
+greatest.common.divisor <- function(n, m) {
+ nf <- primefactors(n)
+ mf <- primefactors(m)
+ p <- sort(unique(c(nf,mf)))
+ nfac <- table(factor(nf, levels=p))
+ mfac <- table(factor(mf, levels=p))
+ prod(p^pmin.int(nfac,mfac))
+}
+
+divisors <- local({
+
+ divisors <- function(n) {
+ p <- primefactors(n)
+ up <- sort(unique(p))
+ k <- table(factor(p, levels=up))
+ return(rest(k, up))
+ }
+
+ rest <- function(kk, uu) {
+ powers <- uu[1]^(0:(kk[1]))
+ if(length(uu) == 1)
+ return(powers)
+ rr <- rest(kk[-1], uu[-1])
+ products <- as.vector(outer(powers, rr, "*"))
+ return(sort(unique(products)))
+ }
+
+ divisors
+})
diff --git a/R/resolve.defaults.R b/R/resolve.defaults.R
new file mode 100755
index 0000000..bde0a6c
--- /dev/null
+++ b/R/resolve.defaults.R
@@ -0,0 +1,178 @@
+#
+# resolve.defaults.R
+#
+# $Revision: 1.32 $ $Date: 2017/01/02 04:48:09 $
+#
+# Resolve conflicts between several sets of defaults
+# Usage:
+# resolve.defaults(list1, list2, list3, .......)
+# where the earlier lists have priority
+#
+
+resolve.defaults <- function(..., .MatchNull=TRUE, .StripNull=FALSE) {
+ # Each argument is a list. Append them.
+ argue <- c(...)
+ # does a NULL value
+ # overwrite a non-null value occurring later in the sequence?
+ if(!.MatchNull) {
+ isnul <- sapply(argue, is.null)
+ argue <- argue[!isnul]
+ }
+ if(!is.null(nam <- names(argue))) {
+ named <- nzchar(nam)
+ arg.unnamed <- argue[!named]
+ arg.named <- argue[named]
+ if(any(discard <- duplicated(names(arg.named))))
+ arg.named <- arg.named[!discard]
+ argue <- append(arg.unnamed, arg.named)
+ }
+ # should a NULL value mean that the argument is missing?
+ if(.StripNull) {
+ isnull <- sapply(argue, is.null)
+ argue <- argue[!isnull]
+ }
+ return(argue)
+}
+
+do.call.without <- function(fun, ..., avoid) {
+ argh <- list(...)
+ nama <- names(argh)
+ if(!is.null(nama))
+ argh <- argh[!(nama %in% avoid)]
+ do.call(fun, argh)
+}
+
+do.call.matched <- function(fun, arglist, funargs,
+ extrargs=NULL,
+ matchfirst=FALSE,
+ sieve=FALSE,
+ skipargs=NULL) {
+ if(!is.function(fun) && !is.character(fun))
+ stop("Internal error: wrong argument type in do.call.matched")
+ if(is.character(fun)) {
+ fname <- fun
+ fun <- get(fname, mode="function")
+ if(!is.function(fun))
+ stop(paste("internal error: function", sQuote(fname), "not found",
+ sep=""))
+ }
+ ## determine list of argument names to be matched
+ if(missing(funargs))
+ funargs <- names(formals(fun))
+ funargs <- c(funargs, extrargs)
+ funargs <- setdiff(funargs, skipargs)
+ ## identify which arguments in the call actually match a formal argument
+ givenargs <- names(arglist)
+ matched <- givenargs %in% funargs
+ # deem the first argument to be matched?
+ if(matchfirst && !nzchar(givenargs[1]))
+ matched[1] <- TRUE
+ # apply 'fun' to matched arguments
+ out <- do.call(fun, arglist[matched])
+ # retain un-matched arguments?
+ if(sieve)
+ out <- list(result=out, otherargs=arglist[!matched])
+ return(out)
+}
+
+
+resolve.1.default <- function(.A, ...) {
+ if(is.character(.A)) {
+ ## .A is the name of the parameter to be returned
+ Aname <- .A
+ res <- resolve.defaults(...)
+ } else if(is.list(.A) && length(.A) == 1) {
+ ## .A is a list giving the name and default value of the parameter
+ Aname <- names(.A)
+ res <- resolve.defaults(..., .A)
+ } else stop("Unrecognised format for .A")
+ hit <- (names(res) == Aname)
+ if(!any(hit)) return(NULL)
+ return(res[[min(which(hit))]])
+}
+
+# extract all the arguments that match '...' rather than a named argument
+
+passthrough <- function(.Fun, ..., .Fname=NULL) {
+ if(is.null(.Fname))
+ .Fname <- deparse(substitute(.Fun))
+ # make a fake call to the named function using the arguments provided
+ cl <- eval(substitute(call(.Fname, ...)))
+ # match the call to the function
+ mc <- match.call(.Fun, cl)
+ # extract the arguments
+ mcargs <- as.list(mc)[-1]
+ # figure out which ones are actually formal arguments of the function
+ nam <- names(formals(.Fun))
+ nam <- setdiff(nam, "...")
+ known <- names(mcargs) %in% nam
+ # return the *other* arguments
+ return(mcargs[!known])
+}
+
+graphicsPars <- local({
+ ## recognised additional arguments to image.default(), axis() etc
+ PlotArgs <- c(
+ "main", "asp", "sub", "axes", "ann",
+ "cex", "font",
+ "cex.axis", "cex.lab", "cex.main", "cex.sub",
+ "col.axis", "col.lab", "col.main", "col.sub",
+ "font.axis", "font.lab", "font.main", "font.sub")
+
+ TextDefArgs <- setdiff(names(formals(text.default)), "...")
+ TextArgs <- c(TextDefArgs, "srt", "family", "xpd")
+
+ TheTable <-
+ list(plot = PlotArgs,
+ image = c(
+ "main", "asp", "sub", "axes", "ann",
+ "xlim", "ylim",
+ "box", # note 'box' is not an argument of image.default
+ "cex", "font",
+ "cex.axis", "cex.lab", "cex.main", "cex.sub",
+ "col.axis", "col.lab", "col.main", "col.sub",
+ "font.axis", "font.lab", "font.main", "font.sub",
+ "claim.title.space"),
+ axis = c(
+ "cex",
+ "cex.axis", "cex.lab",
+ "col.axis", "col.lab",
+ "font.axis", "font.lab",
+ "mgp", "xaxp", "yaxp", "tck", "tcl", "las", "fg", "xpd"),
+ owin = c(
+ "sub",
+ "xlim", "ylim",
+ "cex", "font", "col",
+ "border", "box",
+ "cex.main", "cex.sub",
+ "col.main", "col.sub",
+ "font.main", "font.sub",
+ "xaxs", "yaxs",
+ "claim.title.space"),
+ lines = c("lwd", "lty", "col", "lend", "ljoin", "lmitre"),
+ symbols = c(PlotArgs, "fg", "bg"),
+ text = TextArgs,
+ persp = c("x", "y", "z",
+ "xlim", "ylim", "zlim",
+ "xlab", "ylab", "zlab",
+ "main", "sub",
+ "theta", "phi", "r", "d", "scale",
+ "expand", "col", "border",
+ "ltheta", "lphi", "shade", "box",
+ "axes", "nticks", "ticktype")
+ )
+
+ TheTable$ppp <- unique(c(TheTable$owin,
+ TheTable$symbols,
+ "pch", "cex", "lty", "lwd",
+ "etch",
+ "annotate", "labelmap"))
+
+ graphicsPars <- function(key) {
+ n <- pmatch(key, names(TheTable))
+ if(is.na(n)) return(NULL)
+ return(TheTable[[n]])
+ }
+
+ graphicsPars
+})
diff --git a/R/tapplysum.R b/R/tapplysum.R
new file mode 100644
index 0000000..c520faf
--- /dev/null
+++ b/R/tapplysum.R
@@ -0,0 +1,104 @@
+#' tapplysum.R
+#'
+#' Faster replacement for tapply(..., FUN=sum)
+#'
+#' Adrian Baddeley and Tilman Davies
+#'
+#' $Revision: 1.11 $ $Date: 2016/12/12 09:07:06 $
+
+tapplysum <- function(x, flist, do.names=FALSE, na.rm=TRUE) {
+ stopifnot(is.numeric(x))
+ stopifnot(is.list(flist))
+ stopifnot(all(lengths(flist) == length(x)))
+ stopifnot(all(sapply(flist, is.factor)))
+ nfac <- length(flist)
+ goodx <- is.finite(x)
+ if(na.rm) goodx <- goodx | is.na(x)
+ if(!(nfac %in% 2:3) || !all(goodx)) {
+ y <- tapply(x, flist, sum)
+ y[is.na(y)] <- 0
+ return(y)
+ }
+ ifac <- flist[[1L]]
+ jfac <- flist[[2L]]
+ mi <- length(levels(ifac))
+ mj <- length(levels(jfac))
+ ii <- as.integer(ifac)
+ jj <- as.integer(jfac)
+ if(nfac == 3) {
+ kfac <- flist[[3L]]
+ mk <- length(levels(kfac))
+ kk <- as.integer(kfac)
+ }
+ #' remove NA's
+ if(nfac == 2) {
+ if(anyNA(x) || anyNA(ii) || anyNA(jj)) {
+ ok <- !(is.na(x) | is.na(ii) | is.na(jj))
+ ii <- ii[ok]
+ jj <- jj[ok]
+ x <- x[ok]
+ }
+ } else {
+ if(anyNA(x) || anyNA(ii) || anyNA(jj) || anyNA(kk)) {
+ ok <- !(is.na(x) | is.na(ii) | is.na(jj) | is.na(kk))
+ ii <- ii[ok]
+ jj <- jj[ok]
+ kk <- kk[ok]
+ x <- x[ok]
+ }
+ }
+ n <- length(ii)
+ #'
+ if(nfac == 2) {
+ result <- matrix(0, mi, mj)
+ if(n > 0) {
+ oo <- order(ii, jj)
+ zz <- .C(C_ply2sum,
+ nin = as.integer(n),
+ xin = as.double(x[oo]),
+ iin = as.integer(ii[oo]),
+ jin = as.integer(jj[oo]),
+ nout = as.integer(integer(1L)),
+ xout = as.double(numeric(n)),
+ iout = as.integer(integer(n)),
+ jout = as.integer(integer(n)),
+ PACKAGE = "spatstat.utils")
+ nout <- zz$nout
+ if(nout > 0) {
+ ijout <- cbind(zz$iout, zz$jout)[1:nout,,drop=FALSE]
+ xout <- zz$xout[1:nout]
+ result[ijout] <- xout
+ }
+ }
+ } else {
+ result <- array(0, dim=c(mi, mj, mk))
+ if(n > 0) {
+ oo <- order(ii, jj, kk)
+ zz <- .C(C_ply3sum,
+ nin = as.integer(n),
+ xin = as.double(x[oo]),
+ iin = as.integer(ii[oo]),
+ jin = as.integer(jj[oo]),
+ kin = as.integer(kk[oo]),
+ nout = as.integer(integer(1L)),
+ xout = as.double(numeric(n)),
+ iout = as.integer(integer(n)),
+ jout = as.integer(integer(n)),
+ kout = as.integer(integer(n)),
+ PACKAGE = "spatstat.utils")
+ nout <- zz$nout
+ if(nout > 0) {
+ ijkout <- cbind(zz$iout, zz$jout, zz$kout)[1:nout,,drop=FALSE]
+ xout <- zz$xout[1:nout]
+ result[ijkout] <- xout
+ }
+ }
+ }
+ if(do.names)
+ dimnames(result) <- lapply(flist, levels)
+ return(result)
+}
+
+
+
+
diff --git a/R/utilarg.R b/R/utilarg.R
new file mode 100644
index 0000000..3a25d70
--- /dev/null
+++ b/R/utilarg.R
@@ -0,0 +1,385 @@
+#'
+#' utilarg.R
+#'
+#' Utilities for checking/handling arguments
+#'
+#' $Revision: 1.2 $ $Date: 2016/12/30 03:32:21 $
+#'
+
+"%orifnull%" <- function(a, b) {
+ if(!is.null(a)) return(a)
+ # b is evaluated only now
+ return(b)
+}
+
+check.nvector <- function(v, npoints=NULL, fatal=TRUE, things="data points",
+ naok=FALSE, warn=FALSE, vname, oneok=FALSE) {
+ # vector of numeric values for each point/thing
+ if(missing(vname))
+ vname <- sQuote(deparse(substitute(v)))
+ whinge <- NULL
+ nv <- length(v)
+ if(!is.numeric(v))
+ whinge <- paste(vname, "is not numeric")
+ else if(!is.atomic(v) || !is.null(dim(v))) # vector with attributes
+ whinge <- paste(vname, "is not a vector")
+ else if(!(is.null(npoints) || (nv == npoints)) &&
+ !(oneok && nv == 1))
+ whinge <- paste("The length of", vname,
+ paren(paste0("=", nv)),
+ "should equal the number of", things,
+ paren(paste0("=", npoints)))
+ else if(!naok && anyNA(v))
+ whinge <- paste("Some values of", vname, "are NA or NaN")
+ #
+ if(!is.null(whinge)) {
+ if(fatal) stop(whinge)
+ if(warn) warning(whinge)
+ ans <- FALSE
+ attr(ans, "whinge") <- whinge
+ return(ans)
+ }
+ return(TRUE)
+}
+
+check.nmatrix <- function(m, npoints=NULL, fatal=TRUE, things="data points",
+ naok=FALSE, squarematrix=TRUE, matchto="nrow",
+ warn=FALSE) {
+ ## matrix of values for each thing or each pair of things
+ mname <- sQuote(deparse(substitute(m)))
+ whinge <- NULL
+ if(!is.matrix(m))
+ whinge <- paste(mname, "should be a matrix")
+ else if(squarematrix && (nrow(m) != ncol(m)))
+ whinge <- paste(mname, "should be a square matrix")
+ else if(!naok && anyNA(m))
+ whinge <- paste("Some values of", mname, "are NA or NaN")
+ else if(!is.null(npoints)) {
+ if(matchto=="nrow" && nrow(m) != npoints)
+ whinge <- paste("Number of rows in", mname,
+ paren(paste0("=", nrow(m))),
+ "does not match number of", things,
+ paren(paste0("=", npoints)))
+ else if(matchto=="ncol" && ncol(m) != npoints)
+ whinge <- paste("Number of columns in", mname,
+ paren(paste0("=", ncol(m))),
+ "does not match number of", things,
+ paren(paste0("=", npoints)))
+ }
+ ##
+ if(!is.null(whinge)) {
+ if(fatal) stop(whinge)
+ if(warn) warning(whinge)
+ return(FALSE)
+ }
+ return(TRUE)
+}
+
+check.named.vector <- function(x, nam, context="", namopt=character(0),
+ onError=c("fatal", "null")) {
+ xtitle <- deparse(substitute(x))
+ onError <- match.arg(onError)
+ problem <- check.named.thing(x, nam, namopt, sQuote(xtitle),
+ is.numeric(x), "vector", context,
+ fatal=(onError == "fatal"))
+ if(length(problem) > 0 && onError == "null")
+ return(NULL)
+ opt <- namopt %in% names(x)
+ return(x[c(nam, namopt[opt])])
+}
+
+check.named.list <- function(x, nam, context="", namopt=character(0),
+ onError=c("fatal", "null")) {
+ xtitle <- deparse(substitute(x))
+ onError <- match.arg(onError)
+ problem <- check.named.thing(x, nam, namopt, sQuote(xtitle),
+ is.list(x), "list", context,
+ fatal=(onError == "fatal"))
+ if(length(problem) > 0 && onError == "null")
+ return(NULL)
+ opt <- namopt %in% names(x)
+ return(x[c(nam, namopt[opt])])
+}
+
+check.named.thing <- function(x, nam, namopt=character(0), xtitle=NULL,
+ valid=TRUE, type="object", context="",
+ fatal=TRUE) {
+ if(is.null(xtitle))
+ xtitle <- sQuote(deparse(substitute(x)))
+ # check whether names(x) contains all obligatory names 'nam'
+ # and possibly some of the optional names 'namopt'
+ namesx <- names(x)
+ omitted <- !(nam %in% namesx)
+ foreign <- !(namesx %in% c(nam, namopt))
+ if(valid && !any(omitted) && !any(foreign))
+ return(character(0))
+ # some condition violated
+ if(nzchar(context))
+ xtitle <- paste(context, xtitle)
+ whinge <- paste(xtitle,
+ "must be a named", paste(type, ",", sep=""),
+ "with components", commasep(nam))
+ if(length(namopt) > 0)
+ whinge <- paste(whinge, paren(paste("and optionally", commasep(namopt))))
+ if(any(omitted)) {
+ grizzle <- paste(ngettext(sum(omitted), "parameter", "parameters"),
+ commasep(nam[omitted]),
+ "omitted")
+ whinge <- paste(whinge, grizzle, sep="; ")
+ }
+ if(any(foreign)) {
+ grizzle <- paste(ngettext(sum(foreign), "component", "components"),
+ commasep(namesx[foreign]),
+ "not recognised")
+ whinge <- paste(whinge, grizzle, sep="; ")
+ }
+ if(fatal)
+ stop(whinge, call.=FALSE)
+ return(whinge)
+}
+
+
+validposint <- function(n, caller, fatal=TRUE) {
+ nname <- deparse(substitute(n))
+ if(length(n) != 1 || n != round(n) || n <=0) {
+ if(!fatal)
+ return(FALSE)
+ prefix <- if(!missing(caller)) paste("In ", caller, ", ", sep="") else NULL
+ stop(paste0(prefix, nname, "should be a single positive integer"),
+ call.=FALSE)
+ }
+ return(TRUE)
+}
+
+
+# errors and checks
+
+forbidNA <- function(x, context="", xname, fatal=TRUE, usergiven=TRUE) {
+ if(missing(xname)) xname <- sQuote(deparse(substitute(x)))
+ if(anyNA(x)) {
+ if(usergiven) {
+ # argument came from user
+ offence <- ngettext(length(x), "be NA", "contain NA values")
+ whinge <- paste(context, xname, "must not", offence)
+ } else {
+ # argument was computed internally
+ violates <- ngettext(length(x), "is NA", "contains NA values")
+ whinge <- paste(context, xname, violates)
+ }
+ if(fatal) stop(whinge, call.=FALSE)
+ warning(whinge, call.=FALSE)
+ return(FALSE)
+ }
+ return(TRUE)
+}
+
+check.finite <- function(x, context="", xname, fatal=TRUE, usergiven=TRUE) {
+ if(missing(xname)) xname <- sQuote(deparse(substitute(x)))
+ forbidNA(x, context, xname, fatal=fatal, usergiven=usergiven)
+ if(any(!is.finite(x))) {
+ if(usergiven) {
+ # argument came from user
+ oblige <- ngettext(length(x),
+ "be a finite value", "contain finite values")
+ whinge <- paste(context, xname, "must", oblige)
+ } else {
+ # argument was computed internally
+ violates <- ngettext(length(x),
+ "is not finite", "contains non-finite values")
+ whinge <- paste(context, xname, violates)
+ }
+ if(fatal) stop(whinge, call.=FALSE)
+ warning(whinge, call.=FALSE)
+ return(FALSE)
+ }
+ return(TRUE)
+}
+
+complaining <- function(whinge, fatal=FALSE, value=NULL) {
+ if(fatal) stop(whinge, call.=FALSE)
+ warning(whinge, call.=FALSE)
+ return(value)
+}
+
+check.1.real <- function(x, context="", fatal=TRUE) {
+ xname <- deparse(substitute(x))
+ if(!is.numeric(x) || length(x) != 1) {
+ whinge <- paste(sQuote(xname), "should be a single number")
+ if(nzchar(context)) whinge <- paste(context, whinge)
+ return(complaining(whinge, fatal=fatal, value=FALSE))
+ }
+ return(TRUE)
+}
+
+check.1.integer <- function(x, context="", fatal=TRUE) {
+ xname <- deparse(substitute(x))
+ if(!is.numeric(x) || length(x) != 1 || !is.finite(x) || x %% 1 != 0) {
+ whinge <- paste(sQuote(xname), "should be a single finite integer")
+ if(nzchar(context)) whinge <- paste(context, whinge)
+ return(complaining(whinge, fatal=fatal, value=FALSE))
+ }
+ return(TRUE)
+}
+
+check.1.string <- function(x, context="", fatal=TRUE) {
+ xname <- deparse(substitute(x))
+ if(!is.character(x) || length(x) != 1) {
+ whinge <- paste(sQuote(xname), "should be a single character string")
+ if(nzchar(context)) whinge <- paste(context, whinge)
+ return(complaining(whinge, fatal=fatal, value=FALSE))
+ }
+ return(TRUE)
+}
+
+explain.ifnot <- function(expr, context="") {
+ ex <- deparse(substitute(expr))
+ ans <- expr
+ if(!(is.logical(ans) && length(ans) == 1 && ans))
+ stop(paste(context, "it must be TRUE that", sQuote(ex)), call.=FALSE)
+}
+
+warn.ignored.args <- function(..., context=NULL) {
+ if((narg <- length(list(...))) > 0) {
+ whinge <- paste(narg, "unrecognised",
+ ngettext(narg, "argument was", "arguments were"),
+ "ignored")
+ if(!is.null(context)) whinge <- paste(context, whinge)
+ warning(context)
+ }
+}
+
+trap.extra.arguments <- function(..., .Context="", .Fatal=FALSE) {
+ z <- list(...)
+ if((narg <- length(z)) == 0) return(FALSE)
+ nama <- names(z)
+ named <- nzchar(nama)
+ whinge <- paste(.Context, ":", sep="")
+ if(any(named)) {
+ # some arguments are named: ensure all are named
+ nama <- sQuote(nama)
+ if(!all(named))
+ nama[!named] <- paste("[Arg", 1:length(nama), ,"]", sep="")[!named]
+ whinge <- paste(whinge,
+ "unrecognised",
+ ngettext(narg, "argument", "arguments"),
+ commasep(nama),
+ ngettext(narg, "was", "were"), "ignored")
+ } else {
+ # all arguments unnamed
+ whinge <- paste(whinge,
+ narg, "unrecognised",
+ ngettext(narg, "argument was", "arguments were"),
+ "ignored")
+ }
+ if(.Fatal) stop(whinge, call.=FALSE) else warning(whinge, call.=FALSE)
+ return(TRUE)
+}
+
+
+## replace recognised keywords by other keywords
+mapstrings <- function(x, map=NULL) {
+ if(is.null(map)) return(x)
+ x <- as.character(x)
+ from <- names(map)
+ map <- as.character(map)
+ if(sum(nzchar(from)) != length(map))
+ stop("input names are missing in map", call.=FALSE)
+ if(any(duplicated(from)))
+ stop("input names are duplicated in map", call.=FALSE)
+ i <- match(x, from)
+ hit <- !is.na(i)
+ x[hit] <- map[i[hit]]
+ return(x)
+}
+
+there.is.no.try <- function(...) {
+ #' do, or do not
+ y <- try(..., silent=TRUE)
+ if(inherits(y, "try-error")) return(NULL)
+ return(y)
+}
+
+dont.complain.about <- function(...) {
+ #' prevents code checkers complaining about 'unused variables'
+ #' Typically needed where the variables in question
+ #' are referenced in an expression that will be evaluated elsewhere.
+ return(invisible(NULL))
+}
+
+
+matchNameOrPosition <- function(expected, avail) {
+ if(length(avail) < length(expected))
+ stop("Not enough arguments to match", call.=FALSE)
+ j <- match(expected, avail)
+ if(!anyNA(j)) return(j)
+ everything <- seq_along(avail)
+ for(k in seq_along(expected)) {
+ if(is.na(j[k]))
+ j[k] <- min(setdiff(everything, j[-k]))
+ }
+ return(j)
+}
+
+badprobability <- function(x, NAvalue=NA) {
+ ifelse(is.na(x), NAvalue, !is.finite(x) | x < 0 | x > 1)
+}
+
+# test for equivalence of two functions
+samefunction <- function(f, g) {
+ identical(deparse(f), deparse(g))
+}
+
+#' .................. calls and expressions ..................
+
+
+fakecallstring <- function(fname, parlist) {
+ cl <- do.call(call, append(list(name = fname), parlist))
+ return(format(cl))
+}
+
+dotexpr.to.call <- function(expr, dot="funX", evaluator="eval.fv") {
+ # convert an expression into a function call
+ # replacing "." by the specified variable
+ stopifnot(is.expression(expr))
+ aa <- substitute(substitute(ee, list(.=as.name(d))),
+ list(ee=expr, d=dot))
+ bb <- eval(parse(text=deparse(aa)))
+ cc <- as.call(bb)
+ cc[[1]] <- as.name(evaluator)
+ return(cc)
+}
+
+inject.expr <- function(base, expr) {
+ ## insert an expression inside a call and parse it
+ txt <- sub(".", as.character(expr), as.character(base), fixed=TRUE)
+ parse(text=txt)
+}
+
+
+## Match variable names to objects in 'data' list or environment
+getdataobjects <- function(nama, envir, datalist=NULL, fatal=FALSE) {
+ if(is.null(nama)) return(NULL)
+ stopifnot(is.character(nama))
+ n <- length(nama)
+ y <- vector(mode="list", length=n)
+ names(y) <- nama
+ if(!is.null(datalist)) {
+ hit <- nama %in% names(datalist)
+ if(any(hit))
+ y[hit] <- as.list(datalist)[nama[hit]]
+ external <- unlist(lapply(y, is.null))
+ } else external <- rep(TRUE, n)
+ y[external] <- mget(nama[external], envir=envir,
+ ifnotfound=list(NULL), inherits=TRUE)
+ if(fatal && any(bad <- unlist(lapply(y, is.null)))) {
+ nbad <- sum(bad)
+ stop(paste(ngettext(nbad, "Covariate", "Covariates"),
+ commasep(sQuote(nama[bad])),
+ ngettext(nbad, "was not found", "were not found")),
+ call.=FALSE)
+ }
+ names(y) <- nama
+ attr(y, "external") <- external
+ return(y)
+}
+
diff --git a/R/utilindex.R b/R/utilindex.R
new file mode 100644
index 0000000..b2fcf4c
--- /dev/null
+++ b/R/utilindex.R
@@ -0,0 +1,100 @@
+#' utilindex.R
+#'
+#' Copyright (c) Adrian Baddeley
+#'
+#' GNU Public Licence >= 2
+#'
+
+matchIntegerDataFrames <- function(X, Y, sort=TRUE) {
+ X <- as.data.frame(X)
+ Y <- as.data.frame(Y)
+ stopifnot(ncol(X) == ncol(Y))
+ if(!all(sapply(X, storage.mode) == "integer"))
+ X <- as.data.frame(lapply(X, as.integer))
+ if(!all(sapply(Y, storage.mode) == "integer"))
+ Y <- as.data.frame(lapply(Y, as.integer))
+ ans <- rep(NA_integer_, nrow(X))
+ switch(ncol(X),
+ {
+ #' one column
+ ans <- match(X[,1], Y[,1])
+ },
+ {
+ #' two columns
+ if(sort) {
+ #' order them
+ oX <- order(X[,1], X[,2])
+ oY <- order(Y[,1], Y[,2])
+ XX <- X[oX, , drop=FALSE]
+ YY <- Y[oY, , drop=FALSE]
+ z <- .C(C_CSmatch2int,
+ na = as.integer(nrow(XX)),
+ xa = as.integer(XX[,1]),
+ ya = as.integer(XX[,2]),
+ nb = as.integer(nrow(YY)),
+ xb = as.integer(YY[,1]),
+ yb = as.integer(YY[,2]),
+ match = as.integer(integer(nrow(XX))),
+ PACKAGE="spatstat.utils")
+ zz <- z$match
+ zz[zz == 0] <- NA
+ ans[oX] <- oY[zz]
+ } else {
+ z <- .C(C_CUmatch2int,
+ na = as.integer(nrow(X)),
+ xa = as.integer(X[,1]),
+ ya = as.integer(X[,2]),
+ nb = as.integer(nrow(Y)),
+ xb = as.integer(Y[,1]),
+ yb = as.integer(Y[,2]),
+ match = as.integer(integer(nrow(X))),
+ PACKAGE="spatstat.utils")
+ zz <- z$match
+ zz[zz == 0] <- NA
+ ans <- zz
+ }
+ },
+ {
+ #' three columns
+ if(sort) {
+ #' order them
+ oX <- order(X[,1], X[,2], X[,3])
+ oY <- order(Y[,1], Y[,2], Y[,3])
+ XX <- X[oX, , drop=FALSE]
+ YY <- Y[oY, , drop=FALSE]
+ z <- .C(C_CSmatch3int,
+ na = as.integer(nrow(XX)),
+ xa = as.integer(XX[,1]),
+ ya = as.integer(XX[,2]),
+ za = as.integer(XX[,3]),
+ nb = as.integer(nrow(YY)),
+ xb = as.integer(YY[,1]),
+ yb = as.integer(YY[,2]),
+ zb = as.integer(YY[,3]),
+ match = as.integer(integer(nrow(XX))),
+ PACKAGE="spatstat.utils")
+ zz <- z$match
+ zz[zz == 0] <- NA
+ ans[oX] <- oY[zz]
+ } else {
+ z <- .C(C_CUmatch3int,
+ na = as.integer(nrow(X)),
+ xa = as.integer(X[,1]),
+ ya = as.integer(X[,2]),
+ za = as.integer(X[,3]),
+ nb = as.integer(nrow(Y)),
+ xb = as.integer(Y[,1]),
+ yb = as.integer(Y[,2]),
+ zb = as.integer(Y[,3]),
+ match = as.integer(integer(nrow(X))),
+ PACKAGE="spatstat.utils")
+ zz <- z$match
+ zz[zz == 0] <- NA
+ ans <- zz
+ }
+ },
+ stop("Sorry, not implemented for more than 3 columns", call.=FALSE))
+ return(ans)
+}
+
+
\ No newline at end of file
diff --git a/R/utilmatrix.R b/R/utilmatrix.R
new file mode 100644
index 0000000..05a11ff
--- /dev/null
+++ b/R/utilmatrix.R
@@ -0,0 +1,182 @@
+#'
+#' utilmatrix.R
+#'
+#' Utilities for matrices and arrays
+#'
+#' $Revision: 1.1 $ $Date: 2016/12/30 03:22:46 $
+#'
+
+matrowsum <- function(x) {
+ # was: x %*% rep.int(1, ncol(x))
+ rowSums(x)
+}
+
+matcolsum <- function(x) {
+ # was: rep.int(1, nrow(x)) %*% x
+ colSums(x)
+}
+
+matrowany <- function(x) {
+ # currently faster than apply(x, 1, any) for logical arrays
+ (rowSums(x) > 0)
+}
+
+matrowall <- function(x) {
+ # currently faster than apply(x, 1, all) for logical arrays
+ (rowSums(x) == ncol(x))
+}
+
+matcolany <- function(x) {
+ # currently faster than apply(x, 2, any) for logical arrays
+ (colSums(x) > 0)
+}
+
+matcolall <- function(x) {
+ # currently faster than apply(x, 2, all) for logical arrays
+ (colSums(x) == nrow(x))
+}
+
+########
+ # hm, this is SLOWER
+
+apply23sum <- function(x) {
+ dimx <- dim(x)
+ if(length(dimx) != 3)
+ stop("x is not a 3D array")
+ result <- array(0, dimx[-1])
+
+ nz <- dimx[3]
+ for(k in 1:nz) {
+ result[,k] <- matcolsum(x[,,k])
+ }
+ result
+}
+
+######################
+#
+# matrixsample subsample or supersample a matrix
+#
+
+matrixsample <- function(mat, newdim, phase=c(0,0), scale, na.value=NA) {
+ # 'phase+1' is the position of the [1,1] corner of the new matrix
+ # expressed in the coordinates of the old matrix.
+ # 'scale' is the size of one step in the new matrix,
+ # expressed in the coordinates of the old matrix.
+ # Both 'phase' and 'scale' can take any real value.
+ olddim <- dim(mat)
+ if(missing(scale)) scale <- (olddim - 1)/(newdim - 1)
+ scale <- ensure2vector(scale)
+ newdim <- ensure2vector(newdim)
+ newmat <- matrix(na.value, newdim[1], newdim[2])
+ newrow <- 1:newdim[1]
+ newcol <- 1:newdim[2]
+ oldrow <- round(1 + phase[1] + (newrow-1) * scale[1])
+ oldcol <- round(1 + phase[2] + (newcol-1) * scale[2])
+ oldrow.ok <- (oldrow >= 1) & (oldrow <= olddim[1])
+ oldcol.ok <- (oldcol >= 1) & (oldcol <= olddim[2])
+ newmat[oldrow.ok, oldcol.ok] <- mat[oldrow[oldrow.ok],
+ oldcol[oldcol.ok]]
+ return(newmat)
+}
+
+# wrangle data.frames
+
+findfirstfactor <- function(x) {
+ if(!inherits(x, c("data.frame", "hyperframe")))
+ stop("x should be a data frame or hyperframe")
+ isfac <- unlist(lapply(as.list(x), is.factor))
+ if(!any(isfac))
+ return(NULL)
+ min(which(isfac))
+}
+
+firstfactor <- function(x) {
+ j <- findfirstfactor(x)
+ if(is.null(j)) return(NULL)
+ return(x[, j, drop=TRUE])
+}
+
+assignDFcolumn <- function(x, name, value, ...) { # for use in mapply
+ dx <- list(value)
+ names(dx) <- name
+ data.frame(append(c(as.list(x), dx), list(...)))
+}
+
+blockdiagmatrix <- function(...) {
+ x <- list(...)
+ if(!all(unlist(lapply(x, is.matrix))))
+ stop("Some of the arguments are not matrices", call.=FALSE)
+ nr <- unlist(lapply(x, nrow))
+ nc <- unlist(lapply(x, ncol))
+ result <- matrix(0, sum(nr), sum(nc))
+ rownames(result) <- unlist(lapply(x, rownames))
+ colnames(result) <- unlist(lapply(x, colnames))
+ rowend <- cumsum(nr)
+ rowstart <- c(0, rowend) + 1
+ colend <- cumsum(nc)
+ colstart <- c(0, colend) + 1
+ for(i in seq_along(x))
+ result[ (rowstart[i]):(rowend[i]) , (colstart[i]):(colend[i])] <- x[[i]]
+ return(result)
+}
+
+blockdiagarray <- function(...) {
+ x <- list(...)
+ if(!all(unlist(lapply(x, is.array))))
+ stop("Some of the arguments are not arrays", call.=FALSE)
+ dims <- lapply(x, dim)
+ dims1 <- unlist(lapply(dims, "[", i=1))
+ if(length(dim1 <- unique(dims1)) > 1)
+ stop("Arrays have different extents in first dimension")
+ dims2 <- unlist(lapply(dims, "[", i=2))
+ dims3 <- unlist(lapply(dims, "[", i=3))
+ result <- array(0, dim=c(dim1, sum(dims2), sum(dims3)))
+ dn <- lapply(x, dimnames)
+ dimnames(result)[[2]] <- unlist(lapply(dn, "[[", i=2))
+ dimnames(result)[[3]] <- unlist(lapply(dn, "[[", i=3))
+ rowend <- cumsum(dims2)
+ rowstart <- c(0, rowend) + 1
+ colend <- cumsum(dims3)
+ colstart <- c(0, colend) + 1
+ for(i in seq_along(x))
+ result[ , (rowstart[i]):(rowend[i]) , (colstart[i]):(colend[i])] <- x[[i]]
+ return(result)
+}
+
+asNumericMatrix <- function(x) {
+ ## workaround for strange artefact of as.matrix.data.frame
+ x <- as.matrix(x)
+ storage.mode(x) <- "double"
+ x
+}
+
+
+exceedsMaxArraySize <- function(...) {
+ (prod(as.numeric(c(...))) > .Machine$integer.max)
+}
+
+
+
+indexCartesian <- function(nn) {
+ # enumerate the elements of the Cartesian product of sets,
+ # where nn[i] is the size of the i-th set
+ as.matrix(do.call(expand.grid, lapply(nn, seq_len)))
+}
+
+
+
+ensure3Darray <- function(x) {
+ nd <- length(dim(x))
+ if(nd == 0) {
+ x <- array(x, dim=c(length(x), 1, 1))
+ } else if(nd == 2) {
+ x <- array(x, dim=c(dim(x), 1))
+ } else if(nd > 3) {
+ laterdims <- dim(x)[-(1:3)]
+ if(any(laterdims != 1))
+ stop("Higher-dimensional array cannot be reduced to 3 dimensions")
+ x <- array(x, dim=dim(x)[1:3])
+ }
+ return(x)
+}
+
diff --git a/R/utiloptim.R b/R/utiloptim.R
new file mode 100644
index 0000000..3d21441
--- /dev/null
+++ b/R/utiloptim.R
@@ -0,0 +1,34 @@
+#'
+#' utiloptim.R
+#'
+#' Utilities for optimization
+#'
+#' $Revision: 1.3 $ $Date: 2017/02/09 00:46:51 $
+#'
+
+optimizeWithTrace <- local({
+
+ tracer <- function(x, ..., .TheFunction, .Enviro) {
+ y <- .TheFunction(x, ...)
+ xx <- get("xx", envir=.Enviro)
+ yy <- get("yy", envir=.Enviro)
+ assign("xx", c(xx, as.numeric(x)), envir=.Enviro)
+ assign("yy", c(yy, y), envir=.Enviro)
+ return(y)
+ }
+
+ optimizeWithTrace <- function(f, interval, ...,
+ lower = min(interval), upper = max(interval)) {
+ e <- new.env()
+ assign("xx", numeric(0), envir=e)
+ assign("yy", numeric(0), envir=e)
+ result <- optimize(tracer, lower=lower, upper=upper,
+ ..., .TheFunction=f, .Enviro=e)
+ result$x <- get("xx", envir=e)
+ result$y <- get("yy", envir=e)
+ rm(e)
+ return(result)
+ }
+
+ optimizeWithTrace
+})
diff --git a/R/utilseq.R b/R/utilseq.R
new file mode 100644
index 0000000..ed80ae1
--- /dev/null
+++ b/R/utilseq.R
@@ -0,0 +1,406 @@
+#'
+#' utilseq.R
+#'
+#' Utilities for sequences, vectors, ranges of values
+#'
+#' $Revision: 1.1 $ $Date: 2017/02/12 09:08:55 $
+#'
+
+dropifsingle <- function(x) if(length(x) == 1) x[[1]] else x
+
+#' ............... ordering ......................
+
+# defines the current favorite algorithm for 'order'
+fave.order <- function(x) { sort.list(x, method="quick", na.last=NA) }
+
+# order statistic (for use in lapply calls)
+orderstats <- function(x, k, decreasing=FALSE) {
+ if(decreasing) sort(x, decreasing=TRUE, na.last=TRUE)[k] else sort(x)[k]
+}
+
+# which value is k-th smallest
+orderwhich <- function(x, k, decreasing=FALSE) {
+ if(decreasing) order(x, decreasing=TRUE, na.last=TRUE)[k] else order(x)[k]
+}
+
+
+## ................ reverse cumulative sum .....................
+
+revcumsum <- function(x) {
+ #' equivalent to rev(cumsum(rev(x)))
+ if(is.complex(x)) {
+ a <- revcumsum(Re(x))
+ b <- revcumsum(Im(x))
+ return(complex(real=a, imaginary=b))
+ }
+ n <- length(x)
+ if(identical(storage.mode(x), "integer")) {
+ z <- .C(C_irevcumsum,
+ x=as.integer(x),
+ as.integer(n),
+ PACKAGE = "spatstat.utils")
+ return(z$x)
+ } else {
+ z <- .C(C_drevcumsum,
+ x=as.double(x),
+ as.integer(n),
+ PACKAGE = "spatstat.utils")
+ return(z$x)
+ }
+}
+
+## ............. vectors of length 2 .................
+
+as2vector <- function(x) {
+ ## convert various wacky formats to numeric vector of length 2
+ ## for use as coordinates of a single point.
+ xname <- deparse(substitute(x))
+ if(is.numeric(x)) {
+ if(length(x) != 2)
+ stop(paste(xname, "should have length 2"))
+ return(x)
+ }
+ if(inherits(x, "ppp")) {
+ #' coded so that it works even if spatstat is not loaded
+ if(x$n != 1)
+ stop(paste(xname, "should consist of exactly one point"))
+ return(c(x$x, x$y))
+ }
+ if(is.list(x) && all(c("x", "y") %in% names(x))) {
+ if(length(x$x) != 1) stop(paste0(xname, "$x should have length 1"))
+ if(length(x$y) != 1) stop(paste0(xname, "$y should have length 1"))
+ return(c(x$x, x$y))
+ }
+ stop(paste("Format of", sQuote(xname), "not understood"))
+}
+
+ensure2vector <- function(x) {
+ xname <- deparse(substitute(x))
+ if(!is.numeric(x))
+ stop(paste(xname, "is not numeric"))
+ n <- length(x)
+ if(n == 0 || n > 2)
+ stop(paste(xname, "should be of length 1 or 2"))
+ if(n == 1)
+ return(rep(x,2))
+ return(x)
+}
+
+
+## ............. sequences ...................
+
+prolongseq <- function(x, newrange, step=NULL) {
+ ## Extend a sequence x so that it covers the new range.
+ stopifnot(length(newrange) == 2 && newrange[1] < newrange[2])
+ ## Check 'x' is an evenly-spaced sequence
+ if(length(x) > 1) {
+ dx <- diff(x)
+ if(any(dx <= 0))
+ stop("x must be an increasing sequence")
+ if(diff(range(dx)) > 0.01 * abs(mean(dx)))
+ stop("x must be evenly spaced")
+ }
+ ## Infer step length
+ if(!is.null(step)) {
+ check.1.real(step)
+ stopifnot(step > 0)
+ } else if(length(x) > 1) {
+ step <- mean(dx)
+ } else stop("step is needed when x is a single value")
+
+ ##
+ if(max(x) < newrange[1] || min(x) > newrange[2])
+ stop("x lies entirely outside the desired range")
+
+ ## add or trim data to left
+ if(x[1] > newrange[1]) {
+ leftbit <- seq(from=x[1], to=newrange[1], by= -step)[-1]
+ x <- c(rev(leftbit), x)
+ nleft <- length(leftbit)
+ } else {
+ nx <- length(x)
+ x <- x[x >= newrange[1]]
+ nleft <- length(x) - nx
+ }
+
+ # add or trim data to right
+ nx <- length(x)
+ if(newrange[2] > x[nx]) {
+ rightbit <- seq(from=x[nx], to=newrange[2], by= step)[-1]
+ x <- c(x, rightbit)
+ nright <- length(rightbit)
+ } else {
+ x <- x[x <= newrange[2]]
+ nright <- length(x) - nx
+ }
+ attr(x, "nleft") <- nleft
+ attr(x, "nright") <- nright
+ return(x)
+}
+
+## fill gaps in a sequence
+fillseq <- function(x, step=NULL) {
+ xname <- short.deparse(substitute(x))
+ n <- length(x)
+ if(n <= 1) return(x)
+ rx <- range(x)
+ dx <- diff(x)
+ if(any(dx < 0)) stop(paste(xname, "should be an increasing sequence"),
+ call.=FALSE)
+ ## guess step length
+ if(is.null(step)) {
+ eps <- diff(rx)/1e7
+ step <- min(dx[dx > eps])
+ }
+ ## make new sequence
+ y <- seq(rx[1], rx[2], by=step)
+ ny <- length(y)
+ ## mapping from x to y
+ i <- round((x - rx[1])/step) + 1L
+ i <- pmin(ny, pmax(1, i))
+ return(list(xnew=y, i=i))
+}
+
+geomseq <- function(from, to, length.out) {
+ if(from <= 0 || to <= 0) stop("range limits must be positive")
+ y <- exp(seq(from=log(from), to=log(to), length.out=length.out))
+ y[1] <- from #' avoid numerical error
+ y[length.out] <- to
+ return(y)
+}
+
+## ............. ranges ...................
+
+intersect.ranges <- function(r, s, fatal=TRUE) {
+ if(!is.null(r) && !is.null(s)) {
+ lo <- max(r[1],s[1])
+ hi <- min(r[2],s[2])
+ if(lo <= hi)
+ return(c(lo, hi))
+ }
+ if(fatal) stop("Intersection is empty")
+ return(NULL)
+}
+
+inside.range <- function(x, r) {
+ stopifnot(length(r) == 2 && r[1] <= r[2])
+ return(x >= r[1] & x <= r[2])
+}
+
+check.in.range <- function(x, r, fatal=TRUE) {
+ xname <- deparse(substitute(x))
+ if(inside.range(x, r))
+ return(TRUE)
+ if(fatal)
+ stop(paste(xname, "should be a number between",
+ r[1], "and", r[2]),
+ call.=FALSE)
+ return(FALSE)
+}
+
+startinrange <- function(x0, dx, r) {
+ ## find y = x0 + n * dx such that y \in r
+ if(all(inside.range(x0, r))) return(x0)
+ stopifnot(is.numeric(dx) && length(dx) == 1)
+ y <- x0 + dx * round((mean(r) - x0)/dx)
+ y[!inside.range(y, r)] <- NA
+ return(y)
+}
+
+prettyinside <- function(x, ...) {
+ r <- range(x, na.rm=TRUE)
+ if(diff(r) == 0) return(r[1])
+ p <- pretty(x, ...)
+ ok <- inside.range(p, r)
+ return(p[ok])
+}
+
+prettydiscrete <- function(x, n=10) {
+ nx <- length(x)
+ dx <- nx %/% n
+ if(dx < 1) return(x)
+ i <- 1 + (0:(n-1)) * dx
+ return(x[i])
+}
+
+
+check.range <- function(r, fatal=TRUE) {
+ rname <- deparse(substitute(r))
+ if(is.numeric(r) && identical(r, range(r, na.rm=TRUE)))
+ return(TRUE)
+ if(fatal)
+ stop(paste(rname, "should be a vector of length 2 giving (min, max)"))
+ return(FALSE)
+}
+
+evenly.spaced <- function(x, tol=1e-07) {
+ # test whether x is evenly spaced and increasing
+ dx <- diff(x)
+ if(any(dx <= .Machine$double.eps))
+ return(FALSE)
+ # The following test for equal spacing is used in hist.default
+ if(diff(range(dx)) > tol * mean(dx))
+ return(FALSE)
+ return(TRUE)
+}
+
+equispaced <- function(z, reltol=0.001) {
+ dz <- diff(z)
+ return(diff(range(dz)) < reltol * mean(dz))
+}
+
+
+adjustthinrange <- function(ur,vstep,vr) {
+ if(diff(ur) >= vstep) return(ur)
+ ur <- mean(ur) + c(-1,1) * vstep/2
+ if(ur[1] < vr[1]) ur <- vr[1] + c(0,1)*vstep
+ if(ur[2] > vr[2]) ur <- vr[2] - c(1,0)*vstep
+ return(ur)
+}
+
+fastFindInterval <- function(x, b, labels=FALSE, reltol=0.001, dig.lab = 3L) {
+ nintervals <- length(b) - 1
+ nx <- length(x)
+ if(nx == 0)
+ return(rep(0, nintervals))
+ ##
+ if(equispaced(b, reltol)) {
+ ## breaks are equally spaced
+ zz <- .C(C_fastinterv,
+ x = as.double(x),
+ n = as.integer(nx),
+ brange = as.double(range(b)),
+ nintervals = as.integer(nintervals),
+ y = as.integer(integer(nx)),
+ PACKAGE = "spatstat.utils"
+ )
+ y <- zz$y
+ } else {
+ ## use R's interval search algorithm
+ y <- findInterval(x, b, rightmost.closed=TRUE)
+ }
+ if(labels) {
+ # Digits in labels code copied from base::cut.default (with small adaptions):
+ for (dig in dig.lab:max(12L, dig.lab)) {
+ ch.br <- formatC(0 + b, digits = dig, width = 1L)
+ if (ok <- all(ch.br[-1L] != ch.br[1L:nintervals]))
+ break
+ }
+ blab <- paste0("[",
+ ch.br[1L:nintervals],
+ ",",
+ ch.br[-1],
+ c(rep(")", nintervals-1), "]"))
+ y <- as.integer(y)
+ levels(y) <- as.character(blab)
+ class(y) <- "factor"
+ }
+ return(y)
+}
+
+# ...................................................
+# efficient replacements for ifelse()
+# 'a' and 'b' are single values
+# 'x' and 'y' are vectors of the same length as 'test'
+
+# ifelse(test, a, b)
+ifelseAB <- function(test, a, b) {
+ y <- rep.int(b, length(test))
+ y[test] <- a
+ return(y)
+}
+
+# ifelse(test, a, x)
+ifelseAX <- function(test, a, x) {
+ y <- x
+ y[test] <- a
+ return(y)
+}
+
+# ifelse(test, x, b)
+ifelseXB <- function(test, x, b) {
+ y <- rep.int(b, length(test))
+ y[test] <- x[test]
+ return(y)
+}
+
+# ifelse(test, x, y)
+ifelseXY <- function(test, x, y) {
+ z <- y
+ z[test] <- x[test]
+ return(z)
+}
+
+#.... very special cases ......
+
+# ifelse(test, 1, NA)
+ifelse1NA <- function(test) {
+ y <- as.integer(test)
+ y[!test] <- NA
+ return(y)
+}
+
+# ifelse(test, 0, NA)
+ifelse0NA <- function(test) {
+ nyet <- !test
+ y <- as.integer(nyet)
+ y[nyet] <- NA
+ return(y)
+}
+
+# ifelse(test, -x, x)
+ifelseNegPos <- function(test, x) {
+ y <- x
+ y[test] <- -x[test]
+ return(y)
+}
+
+
+ratiotweak <- function(a, b, overzero=NA, zerozero=overzero) {
+ # map x/0 to 'overzero' and 0/0 to 'zerozero'
+ result <- a/b
+ bzero <- (b == 0)
+ result[ bzero ] <- overzero
+ if(!missing(zerozero)) {
+ abzero <- bzero & (a == 0)
+ result[ abzero ] <- zerozero
+ }
+ return(result)
+}
+
+natozero <- function(x) {
+ #' map NA to zero (e.g. in tapply)
+ x[is.na(x)] <- 0
+ return(x)
+}
+
+insertinlist <- function(x, i, y) {
+ ## insert a possibly longer or shorter list 'y'
+ ## into serial position 'i' in list 'x'
+ n <- length(x)
+ if(n == 0) return(y)
+ m <- seq_len(n)
+ names(m) <- names(x)
+ i <- m[[i]] # convert 'i' to integer index
+ stopifnot(length(i) == 1)
+ if(n == 1) return(y)
+ xleft <- x[seq_len(i-1)]
+ xright <- x[i + seq_len(n-i)]
+ z <- c(xleft, y, xright)
+ return(z)
+}
+
+#' ..... rounding ..............................
+
+dround <- function(x) {
+ round(x, getOption('digits'))
+}
+
+niceround <- function(x, m=c(1,2,5,10)) {
+ expo <- 10^as.integer(floor(log10(x)))
+ y <- m * expo
+ z <- y[which.min(abs(y - x))]
+ return(z)
+}
+
+
diff --git a/R/utiltext.R b/R/utiltext.R
new file mode 100644
index 0000000..1f9006c
--- /dev/null
+++ b/R/utiltext.R
@@ -0,0 +1,508 @@
+#'
+#' utiltext.R
+#'
+#' Utilities for text output, etc
+#'
+#' $Revision: 1.3 $ $Date: 2017/02/13 08:37:21 $
+#'
+
+# text magic
+
+commasep <- function(x, join=" and ", flatten=TRUE) {
+ px <- paste(x)
+ nx <- length(px)
+ if(nx <= 1) return(px)
+ commas <- c(rep(", ", length(px)-2),
+ join,
+ "")
+ out <- paste0(px, commas, collapse=if(flatten) "" else NULL)
+ return(out)
+}
+
+paren <- function(x, type="(") {
+ if(length(x) == 0) return(x)
+ if(identical(type, "")) type <- "blank"
+ switch(type,
+ "(" = {
+ out <- paste("(", x, ")", sep="")
+ },
+ "[" = {
+ out <- paste("[", x, "]", sep="")
+ },
+ "{" = {
+ out <- paste("{", x, "}", sep="")
+ },
+ blank = {
+ out <- paste(x)
+ },
+ stop(paste("Unrecognised parenthesis type:", sQuote(type)))
+ )
+ out
+}
+
+unparen <- function(x) {
+ x <- as.character(x)
+ firstchar <- substr(x, 1, 1)
+ n <- nchar(x)
+ lastchar <- substr(x, n, n)
+ enclosed <- n > 2 & (
+ (firstchar == "(" & lastchar == ")") |
+ (firstchar == "[" & lastchar == "]") |
+ (firstchar == "{" & lastchar == "}") )
+ if(any(enclosed))
+ x[enclosed] <- substr(x[enclosed], 2, n-1)
+ return(x)
+}
+
+strsplitretain <- local({
+ strsplitretain <- function(x, split=",") {
+ ## split strings after occurrence of character b, but retain b
+ y <- strsplit(x, split)
+ lapply(y, addback, b=split)
+ }
+ addback <- function(x, b=",") {
+ n <- length(x)
+ if(n <= 1) x else c(paste0(x[-n], b), x[n])
+ }
+ strsplitretain
+})
+
+truncline <- function(x, nc) {
+ if(length(x) > 1)
+ return(unlist(lapply(as.list(x), truncline, nc=nc)))
+ ## split string into words
+ y <- strsplit(x, " ", fixed=TRUE)[[1]]
+ ## find max number of whole words that take up nc characters
+ maxwords <- max(0, which(cumsum(nchar(y) + 1) <= nc+1))
+ if(maxwords == length(y))
+ return(x)
+ ## truncation will occur.
+ pad <- " [..]"
+ nc <- nc - nchar(pad)
+ maxwords <- max(0, which(cumsum(nchar(y) + 1) <= nc+1))
+ z <- paste(y[seq_len(maxwords)], collapse=" ")
+ d <- nc - nchar(z)
+ if(d < 0)
+ z <- substr(z, 1, nc)
+ z <- paste0(z, pad)
+ return(z)
+}
+
+is.blank <- function(s) {
+ y <- strsplit(s, "")
+ z <- lapply(y, "==", e2=" ")
+ ans <- sapply(z, all)
+ return(ans)
+}
+
+padtowidth <- local({
+
+ blankstring <- function(n) paste(rep(" ", n), collapse="")
+
+ padtowidth <- function(a, b, justify=c("left", "right", "centre")) {
+ justify <- match.arg(justify)
+ if(is.character(b)) b <- nchar(b) else stopifnot(is.numeric(b))
+ extra <- pmax(0, b - nchar(a))
+ rpad <- lpad <- ""
+ switch(justify,
+ left = {
+ rpad <- sapply(extra, blankstring)
+ },
+ right = {
+ lpad <- sapply(extra, blankstring)
+ },
+ centre = {
+ lpad <- sapply(floor(extra/2), blankstring)
+ rpad <- sapply(ceiling(extra/2), blankstring)
+ })
+ result <- paste0(lpad, a, rpad)
+ return(result)
+ }
+
+ padtowidth
+})
+
+## 'split cat'
+## Replacement for 'cat(paste(...))' ensuring a multi-word output string
+## doesn't extend over text margin
+
+splat <- function(..., indent=0) {
+ st <- pasteN(...) # removes NULL arguments without making a space
+ ## split at newline characters, if present
+ ss <- unlist(strsplit(st, "\n"))
+ ##
+ if(is.numeric(indent)) {
+ nindent <- indent
+ indent <- paste(rep(" ", nindent), collapse="")
+ } else if(is.character(indent)) {
+ nindent <- nchar(indent)
+ } else stop("indent should be character or numeric")
+ w <- getOption('width')
+ if(nindent >= w) {
+ warning("indentation is more than the permissible text width: ignored")
+ nindent <- 0
+ }
+ ##
+ if(nindent == 0) {
+ for(ssi in ss)
+ cat(unlist(strsplit(ssi, " ")), fill=TRUE)
+ } else {
+ wfill <- w - nindent
+ for(ssi in ss) {
+ vi <- choptextline(ssi, w=w, indent=indent)
+ for(vij in vi) {
+ cat(indent)
+ cat(vij, fill=wfill)
+ }
+ }
+ }
+ return(invisible(NULL))
+}
+
+choptext <- function(..., prefix="", indent="") {
+ s <- paste(...)
+ ## split at newline characters, if present
+ lines <- unlist(strsplit(s, "\n"))
+ ## cut into pieces that don't overreach width
+ w <- getOption('width')
+ lines <- sapply(lines, choptextline, w=w, prefix=prefix, indent=indent)
+ lines <- unname(as.vector(lines))
+ return(lines)
+}
+
+choptextline <- function(st, w=getOption('width'), prefix="", indent="") {
+ words <- unlist(strsplit(st, " "))
+ nwords <- length(words)
+ wordlengths <- nchar(words)
+ lines <- character(0)
+ prefixlength <- nchar(prefix)
+ indentlength <- nchar(indent)
+ while(nwords > 0) {
+ inset <- prefixlength + indentlength
+ wordends <- cumsum(wordlengths + c(inset, rep(1, nwords-1)))
+ n <- which.max(wordends * (wordends <= w))
+ if(n == 0) n <- 1
+ lines <- c(lines, paste(words[1:n], collapse=" "))
+ words <- words[-(1:n)]
+ wordlengths <- wordlengths[-(1:n)]
+ nwords <- nwords - n
+ prefixlength <- 0
+ }
+ return(lines)
+}
+
+exhibitStringList <- function(prefix, strings) {
+ shortblurb <- paste(prefix, paste(strings, collapse=", "), "\n")
+ if(nchar(shortblurb) < options("width")[[1]]) {
+ cat(shortblurb)
+ } else {
+ cat(paste(prefix,"\n"))
+ splat(" ", paste(strings, collapse=" "))
+ }
+ return(invisible(NULL))
+}
+
+
+## grammar, etc
+
+
+ordinal <- function(k) paste(k, ordinalsuffix(k), sep="")
+
+ordinalsuffix <- function(k) {
+ last <- abs(k) %% 10
+ lasttwo <- abs(k) %% 100
+ isteen <- (lasttwo > 10 & lasttwo < 20)
+ ending <- ifelse(isteen, "th",
+ ifelse(last == 1, "st",
+ ifelse(last == 2, "nd",
+ ifelse(last == 3, "rd",
+ "th"))))
+ return(ending)
+}
+
+articlebeforenumber <- function(k) {
+ k <- abs(k)
+ if(k == 11 || k == 18) return("an")
+ leading <- floor(k/10^floor(log10(k)))
+ if(leading == 8) return("an")
+ return("a")
+}
+
+numalign <- function(i, nmax, zero="0") {
+ stopifnot(i <= nmax)
+ nplaces <- as.integer(ceiling(log10(nmax+1)))
+ out <- paste(rep(zero, nplaces), collapse="")
+ istring <- paste(i)
+ ilen <- nchar(istring)
+ substr(out, nplaces-ilen+1, nplaces) <- istring
+ return(out)
+}
+
+singlestring <- function(s, coll="") {
+ s <- as.character(s)
+ if(length(s) > 1)
+ s <- paste(s, collapse=coll)
+ return(s)
+}
+
+verbalogic <- function(x, op="and") {
+ stopifnot(is.character(x))
+ istrue <- (x == "TRUE")
+ isfalse <- (x == "FALSE")
+ isvariable <- !istrue & !isfalse
+ y <- x[isvariable]
+ switch(op,
+ and={
+ if(any(isfalse))
+ return("FALSE")
+ if(all(istrue))
+ return("TRUE")
+ return(paste(y, collapse=" and "))
+ },
+ or={
+ if(all(isfalse))
+ return("FALSE")
+ if(any(istrue))
+ return("TRUE")
+ return(paste(y, collapse=" or "))
+ },
+ not={
+ x[isfalse] <- "TRUE"
+ x[istrue] <- "FALSE"
+ x[isvariable] <- paste("not {", y, "}")
+ return(x)
+ },
+ stop(paste("Unrecognised operation", sQuote(op))))
+}
+
+sensiblevarname <- function(guess, fallback, maxlen=12) {
+ out <- if(is.character(guess) &&
+ length(guess) == 1 &&
+ make.names(guess) == guess) guess else fallback
+ out <- substr(out, 1, maxlen)
+ return(out)
+}
+
+## deparse() can sometimes be equivalent to dumping the whole object
+short.deparse <- function(x, maxlen=60) {
+ deparse(x,
+ nlines=1,
+ width.cutoff=maxlen,
+ control="delayPromises")
+}
+
+## deparse() can produce multiple lines of text
+flat.deparse <- function(x) {
+ y <- paste(deparse(x), collapse=" ")
+ y <- gsub("\n", " ", y)
+ y <- gsub(" ", "", y)
+ return(y)
+}
+
+good.names <- function(nama, defaults, suffices) {
+ # ensure sensible, unique names
+ stopifnot(is.character(defaults))
+ if(!missing(suffices))
+ defaults <- paste(defaults, suffices, sep="")
+ result <- nama
+ if(is.null(result))
+ result <- defaults
+ else if(any(blank <- !nzchar(result)))
+ result[blank] <- defaults[blank]
+ if(anyDuplicated(result))
+ result <- make.names(result, unique=TRUE)
+ return(result)
+}
+
+nzpaste <- function(..., sep=" ", collapse=NULL) {
+ # Paste only the non-empty strings
+ v <- list(...)
+ ok <- sapply(lapply(v, nzchar), any)
+ do.call(paste, append(v[ok], list(sep=sep, collapse=collapse)))
+}
+
+pasteN <- function(...) {
+ # remove NULL arguments then paste
+ argh <- list(...)
+ argh <- argh[!sapply(argh, is.null)]
+ do.call(paste, argh)
+}
+
+substringcount <- function(x, y) {
+ ## count occurrences of 'x' in 'y'
+ yy <- paste0("a", y, "a")
+ splot <- strsplit(yy, split=x, fixed=TRUE)
+ nhits <- lengths(splot) - 1
+ return(nhits)
+}
+
+is.parseable <- local({
+ is.parseable <- function(x) sapply(x, canparse)
+
+ canparse <- function(z) {
+ !inherits(try(parse(text=z), silent=TRUE), "try-error")
+ }
+
+ is.parseable
+})
+
+make.parseable <- function(x) {
+ if(all(is.parseable(x))) x else make.names(x)
+}
+
+# paste(expression(..)) seems to be broken
+
+paste.expr <- function(x) {
+ unlist(lapply(lapply(x, deparse),
+ paste, collapse=""))
+}
+
+pasteFormula <- function(f) {
+ ## convert formula to a single string
+ sf <- paste(format(f), collapse=" ")
+ ## remove excessive blanks
+ sf <- gsub(" ", " ", sf)
+ sf <- gsub(" ", " ", sf)
+ return(sf)
+}
+
+cat.factor <- function (...) {
+ lll <- list(...)
+ chk <- sapply(lll,is.factor)
+ if(!all(chk))
+ warning("Some arguments were not factors (but were converted to factors)")
+ lll <- lapply(lll,as.data.frame,nm="v1")
+ return(do.call(rbind,lll)[,1])
+}
+
+prange <- function(r) {
+ stopifnot(length(r) == 2)
+ paren(paste(r, collapse=", "), "[")
+}
+
+# gsub(".", replacement, x) but only when "." appears as a variable
+
+gsubdot <- function(replacement, x) {
+ x <- as.character(x)
+ stopifnot(length(x) == 1)
+ # find all positions of "." in x
+ dotpos <- gregexpr("\\.", x)[[1]]
+ if(all(dotpos == -1)) return(x)
+ # find all positions of "." preceded or followed by alphanumeric
+ dotbefore <- gregexpr("\\.[0-9A-Za-z]", x)[[1]]
+ dotafter <- gregexpr("[0-9A-Za-z]\\.", x)[[1]] - 1
+ # exclude them
+ dotpos <- setdiff(dotpos, union(dotbefore, dotafter))
+ #
+ if(length(dotpos) == 0) return(x)
+ lenrep <-length(replacement)
+ while(length(dotpos) > 0) {
+ dp <- dotpos[1]
+ x <- paste0(substr(x, 0, dp-1), replacement, substr(x, dp+1, nchar(x)))
+ dotpos <- dotpos[-1] + lenrep-1
+ }
+ return(x)
+}
+
+
+
+simplenumber <- local({
+
+ iswhole <- function(x, tol=0) { abs(x %% 1) <= tol }
+
+ iszero <- function(x, tol=0) { abs(x) <= tol }
+
+ simplenumber <- function(x, unit = "", multiply="*",
+ tol=.Machine$double.eps) {
+ ## Try to express x as a simple multiple or fraction
+ if(length(x) > 1)
+ return(sapply(as.list(x), simplenumber,
+ unit=unit, multiply=multiply, tol=tol))
+ s <- if(x < 0) "-" else ""
+ x <- abs(x)
+ if(unit == "") {
+ if(iswhole(x, tol)) return(paste0(s, round(x)))
+ for(i in 1:12) {
+ if(iswhole(i/x, tol)) return(paste0(s, i, "/", round(i/x)))
+ if(iswhole(i*x, tol)) return(paste0(s, round(i*x), "/", i))
+ }
+ } else {
+ if(iszero(x, tol)) return("0")
+ if(iszero(x-1, tol)) return(paste0(s,unit))
+ if(iswhole(x, tol)) return(paste0(s, round(x), multiply, unit))
+ if(iswhole(1/x, tol)) return(paste0(s, unit, "/", round(1/x)))
+ for(i in 2:12) {
+ if(iswhole(i/x, tol))
+ return(paste0(s, i, multiply, unit, "/", round(i/x)))
+ if(iswhole(i*x, tol))
+ return(paste0(s, round(i*x), multiply, unit, "/", i))
+ }
+ }
+ return(NULL)
+ }
+
+ simplenumber
+})
+
+
+fontify <- function(x, font="italic") {
+ if(!nzchar(font) || font == "plain")
+ return(x)
+ if(is.character(x))
+ return(paste0(font, "(", x, ")"))
+ if(is.expression(x)) {
+ if((n <- length(x)) > 0) {
+ for(i in 1:n)
+ x[[i]] <- fontify(x[[i]], font)
+ }
+ return(x)
+ }
+ if(is.language(x) || is.numeric(x))
+ return(substitute(f(X), list(f=as.name(font), X=x)))
+ if(all(sapply(x, is.language)))
+ return(lapply(x, fontify))
+ return(NULL)
+}
+
+romansort <- local({
+
+ # sort character strings in order of Roman alphabet
+
+ romansort <- function(x) {
+ if(!is.character(x)) return(sort(x))
+ x <- as.vector(x)
+ ## convert each 'word' to a vector of single characters
+ cc <- strsplit(x, "")
+ ## find position of each character in Roman alphabet
+ mm <- lapply(cc, match, table=c(letters, LETTERS))
+ mmax <- max(unlist(mm), na.rm=TRUE)
+ ## encode
+ nn <- sapply(mm, powercode, base=mmax)
+ ## find ordering
+ oo <- order(nn, na.last=TRUE)
+ return(x[oo])
+ }
+
+ powercode <- function(x, base) sum(x * base^rev((seq_len(length(x))-1)))
+
+ romansort
+})
+
+variablesintext <- function(x) all.vars(as.expression(parse(text=x)))
+
+## convert numeric matrix to character, and blank out lower sub-diagonal.
+uptrimat <- function(x) {
+ stopifnot(is.matrix(x))
+ x[] <- as.character(x)
+ x[row(x) > col(x)] <- ""
+ return(noquote(x))
+}
+
+## convert lty codes to text
+lty2char <- function(i) {
+ if(is.numeric(i)) c("blank", "solid", "dashed", "dotted",
+ "dotdash", "longdash", "twodash")[(i %% 7) + 1] else i
+}
+
+
+
diff --git a/R/xycircle.R b/R/xycircle.R
new file mode 100644
index 0000000..0669fd8
--- /dev/null
+++ b/R/xycircle.R
@@ -0,0 +1,137 @@
+#'
+#' xycircle.R
+#'
+#' Low-level utilities for circle geometry
+#'
+#' $Revision: 1.9 $ $Date: 2017/02/20 02:03:33 $
+#'
+
+xysegXcircle <- function(xcentres, ycentres, radii, x0, y0, x1, y1,
+ check=TRUE) {
+ #' 'Cross' version
+ #' Find intersections between circles and segments
+ #' for all combinations of centres, radii and segments.
+ #'
+ #' xcentres, ycentres: numeric vectors of coordinates of centres
+ #' radii: numeric vector of radii
+ #' x0, y0, x1, y1: numeric vectors of segment endpoint coordinates
+ #'
+ if(check)
+ stopifnot(is.numeric(xcentres),
+ is.numeric(ycentres),
+ length(xcentres) == length(ycentres),
+ is.numeric(radii),
+ is.numeric(x0),
+ is.numeric(y0),
+ is.numeric(x1),
+ is.numeric(y1),
+ length(x0) == length(y0),
+ length(x1) == length(y1),
+ length(x0) == length(x1))
+ storage.mode(xcentres) <- storage.mode(ycentres) <- "double"
+ storage.mode(x0) <- storage.mode(y0) <- "double"
+ storage.mode(x1) <- storage.mode(y1) <- "double"
+ storage.mode(radii) <- "double"
+ z <- .Call(C_circXseg,
+ XC = xcentres,
+ YC = ycentres,
+ R = radii,
+ X0 = x0,
+ Y0 = y0,
+ X1 = x1,
+ Y1 = y1,
+ PACKAGE="spatstat.utils")
+ result <- as.data.frame(z)
+ #' indices i, j, k specify provenance of each intersection point
+ #' i = centre, j = segment, k = radius
+ names(result) <- c("x", "y", "i", "j", "k", "sinalpha")
+ return(result)
+}
+
+xysegMcircle <- function(xcentres, ycentres, radmat, x0, y0, x1, y1,
+ check=TRUE) {
+ #' 'Matrix' version
+ #' Find intersections between circles and segments
+ #' where radii are given in a matrix with rows corresponding to centres.
+ #'
+ #' xcentres, ycentres: numeric vectors of coordinates of centres
+ #' radmat: matrix of radii (rows correspond to centres)
+ #' x0, y0, x1, y1: numeric vectors of segment endpoint coordinates
+ #'
+ if(check)
+ stopifnot(is.numeric(xcentres),
+ is.numeric(ycentres),
+ length(xcentres) == length(ycentres),
+ is.numeric(radmat),
+ is.matrix(radmat),
+ nrow(radmat) == length(xcentres),
+ is.numeric(x0),
+ is.numeric(y0),
+ is.numeric(x1),
+ is.numeric(y1),
+ length(x0) == length(y0),
+ length(x1) == length(y1),
+ length(x0) == length(x1))
+ storage.mode(xcentres) <- storage.mode(ycentres) <- "double"
+ storage.mode(x0) <- storage.mode(y0) <- "double"
+ storage.mode(x1) <- storage.mode(y1) <- "double"
+ storage.mode(radmat) <- "double"
+ z <- .Call(C_circMseg,
+ XC = xcentres,
+ YC = ycentres,
+ R = radmat,
+ X0 = x0,
+ Y0 = y0,
+ X1 = x1,
+ Y1 = y1,
+ PACKAGE="spatstat.utils")
+ result <- as.data.frame(z)
+ #' indices i, j, k specify provenance of each intersection point
+ #' i = centre, j = segment, k = radius
+ names(result) <- c("x", "y", "i", "j", "k", "sinalpha")
+ return(result)
+}
+
+xysegPcircle <- function(xc, yc, rc, x0, y0, x1, y1,
+ check=TRUE) {
+ #' 'Parallel' version
+ #' Find intersections between circles and segments
+ #' for circles with centres (xc, yc) and radii (rc) corresponding.
+ #'
+ #' xc, y: numeric vectors of coordinates of centres
+ #' rc: numeric vector of radii (corresponding to xc, yc)
+ #' x0, y0, x1, y1: numeric vectors of segment endpoint coordinates
+ #'
+ if(check)
+ stopifnot(is.numeric(xc),
+ is.numeric(yc),
+ length(xc) == length(yc),
+ is.numeric(rc),
+ length(rc) == length(xc),
+ is.numeric(x0),
+ is.numeric(y0),
+ is.numeric(x1),
+ is.numeric(y1),
+ length(x0) == length(y0),
+ length(x1) == length(y1),
+ length(x0) == length(x1))
+ storage.mode(xc) <- storage.mode(yc) <- "double"
+ storage.mode(x0) <- storage.mode(y0) <- "double"
+ storage.mode(x1) <- storage.mode(y1) <- "double"
+ storage.mode(rc) <- "double"
+ z <- .Call(C_circXseg,
+ XC = xc,
+ YC = yc,
+ RC = rc,
+ X0 = x0,
+ Y0 = y0,
+ X1 = x1,
+ Y1 = y1,
+ PACKAGE="spatstat.utils")
+ result <- as.data.frame(z)
+ #' indices i, j specify provenance of each intersection point
+ #' i = circle, j = segment
+ names(result) <- c("x", "y", "i", "j", "sinalpha")
+ return(result)
+}
+
diff --git a/R/xypolygon.R b/R/xypolygon.R
new file mode 100755
index 0000000..ff611ef
--- /dev/null
+++ b/R/xypolygon.R
@@ -0,0 +1,446 @@
+#
+# xypolygon.R
+#
+# $Revision: 1.67 $ $Date: 2017/01/08 00:38:10 $
+#
+# low-level functions defined for polygons in list(x,y) format
+#
+verify.xypolygon <- function(p, fatal=TRUE) {
+ whinge <- NULL
+ if(!is.list(p) || !all(c("x","y") %in% names(p)))
+ whinge <- "polygon must be a list with components x and y"
+ else if(is.null(p$x) || is.null(p$y) || !is.numeric(p$x) || !is.numeric(p$y))
+ whinge <- "components x and y must be numeric vectors"
+ else if(anyNA(p$x) || anyNA(p$y))
+ whinge <- "x and y coordinates must not contain NA values"
+ else if(length(p$x) != length(p$y))
+ whinge <- "lengths of x and y vectors unequal"
+ else if(length(p$x) < 3)
+ whinge <- "need at least 3 vertices for a polygon"
+ ok <- is.null(whinge)
+ if(!ok && fatal)
+ stop(whinge)
+ return(ok)
+}
+
+inside.xypolygon <- function(pts, polly, test01=TRUE, method="C") {
+ # pts: list(x,y) points to be tested
+ # polly: list(x,y) vertices of a single polygon (n joins to 1)
+ # test01: logical - if TRUE, test whether all values in output are 0 or 1
+ pts <- xy.coords(pts, NULL)
+ verify.xypolygon(polly)
+
+ x <- pts$x
+ y <- pts$y
+ xp <- polly$x
+ yp <- polly$y
+
+ full.npts <- npts <- length(x)
+ nedges <- length(xp) # sic
+
+ # Check for points (x,y) that coincide with vertices (xp, yp)
+ # Handle them separately
+ z <- .C(C_Cmatchxy,
+ na=as.integer(npts),
+ xa=as.double(x),
+ ya=as.double(y),
+ nb=as.integer(nedges),
+ xb=as.double(xp),
+ yb=as.double(yp),
+ match=as.integer(integer(npts)),
+ PACKAGE = "spatstat.utils")
+ is.vertex <- (z$match != 0)
+ retain <- !is.vertex
+ # Remove vertices from subsequent consideration; replace them later
+ if(vertices.present <- !all(retain)) {
+ x <- x[retain]
+ y <- y[retain]
+ npts <- sum(retain)
+ }
+
+ #------------- MAIN ALGORITHM -------------------------------
+ score <- numeric(npts)
+ on.boundary <- rep.int(FALSE, npts)
+
+ if(anyretain<- any(retain)) {
+ switch(method,
+ C={
+ #------------------ call C routine ------------------
+ temp <- .C(C_inxyp,
+ x=as.double(x),
+ y=as.double(y),
+ xp=as.double(xp),
+ yp=as.double(yp),
+ npts=as.integer(npts),
+ nedges=as.integer(nedges),
+ score=as.integer(score),
+ onbndry=as.integer(on.boundary),
+ PACKAGE = "spatstat.utils")
+ score <- temp$score/2
+ on.boundary <- as.logical(temp$onbndry)
+ },
+# Fortran code removed!
+# Fortran={
+# #------------------ call Fortran routine ------------------
+# temp <- DOTFortran("inxyp",
+# x=as.double(x),
+# y=as.double(y),
+# xp=as.double(xp),
+# yp=as.double(yp),
+# npts=as.integer(npts),
+# nedges=as.integer(nedges),
+# score=as.double(score),
+# onbndry=as.logical(on.boundary))
+# score <- temp$score
+# on.boundary <- temp$onbndry
+# },
+ interpreted={
+ #----------------- original interpreted code --------------
+ for(i in 1:nedges) {
+ x0 <- xp[i]
+ y0 <- yp[i]
+ x1 <- if(i == nedges) xp[1] else xp[i+1]
+ y1 <- if(i == nedges) yp[1] else yp[i+1]
+ dx <- x1 - x0
+ dy <- y1 - y0
+ if(dx < 0) {
+ # upper edge
+ xcriterion <- (x - x0) * (x - x1)
+ consider <- (xcriterion <= 0)
+ if(any(consider)) {
+ ycriterion <-
+ y[consider] * dx - x[consider] * dy + x0 * dy - y0 * dx
+ # closed inequality
+ contrib <- (ycriterion >= 0) *
+ ifelseAB(xcriterion[consider] == 0, 1/2, 1)
+ # positive edge sign
+ score[consider] <- score[consider] + contrib
+ # detect whether any point lies on this segment
+ on.boundary[consider] <-
+ on.boundary[consider] | (ycriterion == 0)
+ }
+ } else if(dx > 0) {
+ # lower edge
+ xcriterion <- (x - x0) * (x - x1)
+ consider <- (xcriterion <= 0)
+ if(any(consider)) {
+ ycriterion <-
+ y[consider] * dx - x[consider] * dy + x0 * dy - y0 * dx
+ # open inequality
+ contrib <- (ycriterion < 0) *
+ ifelseAB(xcriterion[consider] == 0, 1/2, 1)
+ # negative edge sign
+ score[consider] <- score[consider] - contrib
+ # detect whether any point lies on this segment
+ on.boundary[consider] <-
+ on.boundary[consider] | (ycriterion == 0)
+ }
+ } else {
+ # vertical edge
+ consider <- (x == x0)
+ if(any(consider)) {
+ # zero score
+ # detect whether any point lies on this segment
+ yconsider <- y[consider]
+ ycriterion <- (yconsider - y0) * (yconsider - y1)
+ on.boundary[consider] <-
+ on.boundary[consider] | (ycriterion <= 0)
+ }
+ }
+ }
+ },
+ stop(paste("Unrecognised choice for", sQuote("method")))
+ )
+ }
+
+ #------------------- END SWITCH ------------------------------
+
+ # replace any polygon vertices that were temporarily removed
+ if(vertices.present) {
+ full.score <- numeric(full.npts)
+ full.on.boundary <- rep.int(FALSE, full.npts)
+ if(anyretain) {
+ full.score[retain] <- score
+ full.on.boundary[retain] <- on.boundary
+ }
+ full.score[is.vertex] <- 1
+ full.on.boundary[is.vertex] <- TRUE
+ score <- full.score
+ on.boundary <- full.on.boundary
+ npts <- full.npts
+
+ }
+
+ #-------------------------------------------------
+
+ # any point recognised as lying on the boundary gets score 1.
+ score[on.boundary] <- 1
+
+ if(test01) {
+ # check sanity
+ if(!all((score == 0) | (score == 1)))
+ warning("internal error: some scores are not equal to 0 or 1")
+ }
+
+ attr(score, "on.boundary") <- on.boundary
+
+ return(score)
+}
+
+
+is.hole.xypolygon <- function(polly) {
+ h <- polly$hole
+ if(is.null(h))
+ h <- (Area.xypolygon(polly) < 0)
+ return(h)
+}
+
+Area.xypolygon <- function(polly) {
+ #
+ # polly: list(x,y) vertices of a single polygon (n joins to 1)
+ #
+
+ # area could be pre-calculated
+ if(!is.null(pa <- polly$area) && is.numeric(pa) && length(pa)==1)
+ return(pa)
+
+ # else calculate
+ verify.xypolygon(polly)
+ xp <- polly$x
+ yp <- polly$y
+
+ nedges <- length(xp) # sic
+
+ # place x axis below polygon
+ yp <- yp - min(yp)
+
+ # join vertex n to vertex 1
+ nxt <- c(2:nedges, 1)
+
+ # x step, WITH sign
+ dx <- xp[nxt] - xp
+
+ # average height
+ ym <- (yp + yp[nxt])/2
+
+ -sum(dx * ym)
+}
+
+bdrylength.xypolygon <- function(polly) {
+ verify.xypolygon(polly)
+ xp <- polly$x
+ yp <- polly$y
+ nedges <- length(xp)
+ nxt <- c(2:nedges, 1)
+ dx <- xp[nxt] - xp
+ dy <- yp[nxt] - yp
+ sum(sqrt(dx^2 + dy^2))
+}
+
+reverse.xypolygon <- function(p, adjust=FALSE) {
+ # reverse the order of vertices
+ # (=> change sign of polygon)
+ verify.xypolygon(p)
+ p$x <- rev(p$x)
+ p$y <- rev(p$y)
+ if(adjust) {
+ if(!is.null(p$hole)) p$hole <- !p$hole
+ if(!is.null(p$area)) p$area <- -p$area
+ }
+ return(p)
+}
+
+overlap.xypolygon <- function(P, Q) {
+ # compute area of overlap of two simple closed polygons
+ verify.xypolygon(P)
+ verify.xypolygon(Q)
+
+ xp <- P$x
+ yp <- P$y
+ np <- length(xp)
+ nextp <- c(2:np, 1)
+
+ xq <- Q$x
+ yq <- Q$y
+ nq <- length(xq)
+ nextq <- c(2:nq, 1)
+
+ # adjust y coordinates so all are nonnegative
+ ylow <- min(c(yp,yq))
+ yp <- yp - ylow
+ yq <- yq - ylow
+
+ area <- 0
+ for(i in 1:np) {
+ ii <- c(i, nextp[i])
+ xpii <- xp[ii]
+ ypii <- yp[ii]
+ for(j in 1:nq) {
+ jj <- c(j, nextq[j])
+ area <- area +
+ overlap.trapezium(xpii, ypii, xq[jj], yq[jj])
+ }
+ }
+ return(area)
+}
+
+overlap.trapezium <- function(xa, ya, xb, yb, verb=FALSE) {
+ # compute area of overlap of two trapezia
+ # which have same baseline y = 0
+ #
+ # first trapezium has vertices
+ # (xa[1], 0), (xa[1], ya[1]), (xa[2], ya[2]), (xa[2], 0).
+ # Similarly for second trapezium
+
+ # Test for vertical edges
+ dxa <- diff(xa)
+ dxb <- diff(xb)
+ if(dxa == 0 || dxb == 0)
+ return(0)
+
+ # Order x coordinates, x0 < x1
+ if(dxa > 0) {
+ signa <- 1
+ lefta <- 1
+ righta <- 2
+ if(verb) cat("A is positive\n")
+ } else {
+ signa <- -1
+ lefta <- 2
+ righta <- 1
+ if(verb) cat("A is negative\n")
+ }
+ if(dxb > 0) {
+ signb <- 1
+ leftb <- 1
+ rightb <- 2
+ if(verb) cat("B is positive\n")
+ } else {
+ signb <- -1
+ leftb <- 2
+ rightb <- 1
+ if(verb) cat("B is negative\n")
+ }
+ signfactor <- signa * signb # actually (-signa) * (-signb)
+ if(verb) cat(paste("sign factor =", signfactor, "\n"))
+
+ # Intersect x ranges
+ x0 <- max(xa[lefta], xb[leftb])
+ x1 <- min(xa[righta], xb[rightb])
+ if(x0 >= x1)
+ return(0)
+ if(verb) {
+ cat(paste("Intersection of x ranges: [", x0, ",", x1, "]\n"))
+ abline(v=x0, lty=3)
+ abline(v=x1, lty=3)
+ }
+
+ # Compute associated y coordinates
+ slopea <- diff(ya)/diff(xa)
+ y0a <- ya[lefta] + slopea * (x0-xa[lefta])
+ y1a <- ya[lefta] + slopea * (x1-xa[lefta])
+ slopeb <- diff(yb)/diff(xb)
+ y0b <- yb[leftb] + slopeb * (x0-xb[leftb])
+ y1b <- yb[leftb] + slopeb * (x1-xb[leftb])
+
+ # Determine whether upper edges intersect
+ # if not, intersection is a single trapezium
+ # if so, intersection is a union of two trapezia
+
+ yd0 <- y0b - y0a
+ yd1 <- y1b - y1a
+ if(yd0 * yd1 >= 0) {
+ # edges do not intersect
+ areaT <- (x1 - x0) * (min(y1a,y1b) + min(y0a,y0b))/2
+ if(verb) cat(paste("Edges do not intersect\n"))
+ } else {
+ # edges do intersect
+ # find intersection
+ xint <- x0 + (x1-x0) * abs(yd0/(yd1 - yd0))
+ yint <- y0a + slopea * (xint - x0)
+ if(verb) {
+ cat(paste("Edges intersect at (", xint, ",", yint, ")\n"))
+ points(xint, yint, cex=2, pch="O")
+ }
+ # evaluate left trapezium
+ left <- (xint - x0) * (min(y0a, y0b) + yint)/2
+ # evaluate right trapezium
+ right <- (x1 - xint) * (min(y1a, y1b) + yint)/2
+ areaT <- left + right
+ if(verb)
+ cat(paste("Left area = ", left, ", right=", right, "\n"))
+ }
+
+ # return area of intersection multiplied by signs
+ return(signfactor * areaT)
+}
+
+
+simplify.xypolygon <- function(p, dmin) {
+ verify.xypolygon(p)
+ x <- p$x
+ y <- p$y
+ n <- length(x)
+ if(n <= 3) return(p)
+ dmin2 <- dmin^2
+ # edge lengths: len2[i] is distance from i to i+1
+ len2 <- (x - c(x[-1], x[1]))^2 + (y - c(y[-1],y[1]))^2
+ #
+ while(n > 3 && any(len2 < dmin2)) {
+ # delete the shortest edge
+ kill <- which.min(len2)
+ # edge from 'kill' to 'kill+1' will be removed
+ # Replacement vertex is midpoint of segment
+ left <- if(kill == 1) n else (kill - 1)
+ killplus1 <- if(kill == n) 1 else (kill + 1)
+ right <- if(killplus1 == n) 1 else (killplus1 + 1)
+ xmid <- (x[kill]+x[killplus1])/2
+ ymid <- (y[kill]+y[killplus1])/2
+ d2leftmid <- (xmid-x[left])^2+(ymid-y[left])^2
+ d2midright <- (xmid-x[right])^2+(ymid-y[right])^2
+ # adjust vectors: first replace segment endpoints without deleting
+ x[kill] <- xmid
+ y[kill] <- ymid
+ x[killplus1] <- xmid
+ y[killplus1] <- ymid
+ len2[left] <- d2leftmid
+ len2[kill] <- 0
+ len2[killplus1] <- d2midright
+ # now delete
+ x <- x[-kill]
+ y <- y[-kill]
+ len2 <- len2[-kill]
+ n <- n-1
+ }
+ #
+ p$x <- x
+ p$y <- y
+ p$area <- Area.xypolygon(p[c("x","y")])
+ return(p)
+}
+
+inside.triangle <- function(x, y, xx, yy) {
+ # test whether points x[i], y[i] lie in triangle xx[1:3], yy[1:3]
+ # using barycentric coordinates
+ # vector 0 is edge from A to C
+ v0x <- xx[3] - xx[1]
+ v0y <- yy[3] - yy[1]
+ # vector 1 is edge from A to B
+ v1x <- xx[2] - xx[1]
+ v1y <- yy[2] - yy[1]
+ # vector 2 is from vertex A to point P
+ v2x <- x - xx[1]
+ v2y <- y - yy[1]
+ # inner products
+ dot00 <- v0x^2 + v0y^2
+ dot01 <- v0x * v1x + v0y * v1y
+ dot02 <- v0x * v2x + v0y * v2y
+ dot11 <- v1x^2 + v1y^2
+ dot12 <- v1x * v2x + v1y * v2y
+ # unnormalised barycentric coordinates
+ Denom <- dot00 * dot11 - dot01 * dot01
+ u <- dot11 * dot02 - dot01 * dot12
+ v <- dot00 * dot12 - dot01 * dot02
+ # test
+ return((u > 0) & (v > 0) & (u + v < Denom))
+}
diff --git a/R/xysegment.R b/R/xysegment.R
new file mode 100755
index 0000000..7f03bfe
--- /dev/null
+++ b/R/xysegment.R
@@ -0,0 +1,232 @@
+#
+# xysegment.S
+#
+# $Revision: 1.18 $ $Date: 2017/02/20 06:27:10 $
+#
+# Low level utilities for analytic geometry for line segments
+#
+# author: Adrian Baddeley 2001
+# from an original by Rob Foxall 1997
+#
+# distpl(p, l)
+# distance from a single point p = (xp, yp)
+# to a single line segment l = (x1, y1, x2, y2)
+#
+# distppl(p, l)
+# distances from each of a list of points p[i,]
+# to a single line segment l = (x1, y1, x2, y2)
+# [uses only vector parallel ops]
+#
+# distppll(p, l)
+# distances from each of a list of points p[i,]
+# to each of a list of line segments l[i,]
+# [interpreted code uses large matrices and 'outer()']
+# [Fortran implementation included!]
+
+distpl <- function(p, l) {
+ xp <- p[1]
+ yp <- p[2]
+ dx <- l[3]-l[1]
+ dy <- l[4]-l[2]
+ leng <- sqrt(dx^2 + dy^2)
+ # vector from 1st endpoint to p
+ xpl <- xp - l[1]
+ ypl <- yp - l[2]
+ # distance from p to 1st & 2nd endpoints
+ d1 <- sqrt(xpl^2 + ypl^2)
+ d2 <- sqrt((xp-l[3])^2 + (yp-l[4])^2)
+ dmin <- min(d1,d2)
+ # test for zero length
+ if(leng < .Machine$double.eps)
+ return(dmin)
+ # rotation sine & cosine
+ co <- dx/leng
+ si <- dy/leng
+ # back-rotated coords of p
+ xpr <- co * xpl + si * ypl
+ ypr <- - si * xpl + co * ypl
+ # test
+ if(xpr >= 0 && xpr <= leng)
+ dmin <- min(dmin, abs(ypr))
+ return(dmin)
+}
+
+distppl <- function(p, l) {
+ xp <- p[,1]
+ yp <- p[,2]
+ dx <- l[3]-l[1]
+ dy <- l[4]-l[2]
+ leng <- sqrt(dx^2 + dy^2)
+ # vector from 1st endpoint to p
+ xpl <- xp - l[1]
+ ypl <- yp - l[2]
+ # distance from p to 1st & 2nd endpoints
+ d1 <- sqrt(xpl^2 + ypl^2)
+ d2 <- sqrt((xp-l[3])^2 + (yp-l[4])^2)
+ dmin <- pmin.int(d1,d2)
+ # test for zero length
+ if(leng < .Machine$double.eps)
+ return(dmin)
+ # rotation sine & cosine
+ co <- dx/leng
+ si <- dy/leng
+ # back-rotated coords of p
+ xpr <- co * xpl + si * ypl
+ ypr <- - si * xpl + co * ypl
+ # ypr is perpendicular distance to infinite line
+ # Applies only when xp, yp in the middle
+ middle <- (xpr >= 0 & xpr <= leng)
+ if(any(middle))
+ dmin[middle] <- pmin.int(dmin[middle], abs(ypr[middle]))
+
+ return(dmin)
+}
+
+distppll <- function(p, l, mintype=0,
+ method=c("C", "Fortran", "interpreted"), listit=FALSE) {
+ np <- nrow(p)
+ nl <- nrow(l)
+ xp <- p[,1]
+ yp <- p[,2]
+ if(is.na(match(mintype,0:2)))
+ stop(paste("Argument", sQuote("mintype"), "must be 0, 1 or 2.\n"))
+ method <- match.arg(method)
+ if(method == "Fortran") {
+ warning("method='Fortran' is no longer supported; method='C' was used")
+ method <- "C"
+ }
+ switch(method,
+ interpreted={
+ dx <- l[,3]-l[,1]
+ dy <- l[,4]-l[,2]
+ # segment lengths
+ leng <- sqrt(dx^2 + dy^2)
+ # rotation sines & cosines
+ co <- dx/leng
+ si <- dy/leng
+ co <- matrix(co, nrow=np, ncol=nl, byrow=TRUE)
+ si <- matrix(si, nrow=np, ncol=nl, byrow=TRUE)
+ # matrix of squared distances from p[i] to 1st endpoint of segment j
+ xp.x1 <- outer(xp, l[,1], "-")
+ yp.y1 <- outer(yp, l[,2], "-")
+ d1 <- xp.x1^2 + yp.y1^2
+ # ditto for 2nd endpoint
+ xp.x2 <- outer(xp, l[,3], "-")
+ yp.y2 <- outer(yp, l[,4], "-")
+ d2 <- xp.x2^2 + yp.y2^2
+ # for each (i,j) rotate p[i] around 1st endpoint of segment j
+ # so that line segment coincides with x axis
+ xpr <- xp.x1 * co + yp.y1 * si
+ ypr <- - xp.x1 * si + yp.y1 * co
+ d3 <- ypr^2
+ # test
+ lenf <- matrix(leng, nrow=np, ncol=nl, byrow=TRUE)
+ zero <- (lenf < .Machine$double.eps)
+ outside <- (zero | xpr < 0 | xpr > lenf)
+ if(any(outside))
+ d3[outside] <- Inf
+
+ dsq <- matrix(pmin.int(d1, d2, d3),nrow=np, ncol=nl)
+ d <- sqrt(dsq)
+ if(mintype >= 1)
+ min.d <- apply(d, 1, min)
+ if(mintype == 2)
+ min.which <- apply(d, 1, which.min)
+ },
+# Fortran code removed!
+# Fortran={
+# eps <- .Machine$double.eps
+# if(mintype > 0) {
+# big <- sqrt(2)*diff(range(c(p,l)))
+# xmin <- rep.int(big,np)
+# } else {
+# xmin <- 1
+# }
+# n2 <- if(mintype > 1) np else 1
+# temp <- DOTFortran("dppll",
+# x=as.double(xp),
+# y=as.double(yp),
+# l1=as.double(l[,1]),
+# l2=as.double(l[,2]),
+# l3=as.double(l[,3]),
+# l4=as.double(l[,4]),
+# np=as.integer(np),
+# nl=as.integer(nl),
+# eps=as.double(eps),
+# mint=as.integer(mintype),
+# rslt=double(np*nl),
+# xmin=as.double(xmin),
+# jmin=integer(n2))
+# d <- matrix(temp$rslt, nrow=np, ncol=nl)
+# if(mintype >= 1)
+# min.d <- temp$xmin
+# if(mintype == 2)
+# min.which <- temp$jmin
+# },
+ C = {
+ eps <- .Machine$double.eps
+ temp <- .C(C_prdist2segs,
+ x=as.double(xp),
+ y=as.double(yp),
+ npoints =as.integer(np),
+ x0=as.double(l[,1]),
+ y0=as.double(l[,2]),
+ x1=as.double(l[,3]),
+ y1=as.double(l[,4]),
+ nsegments=as.integer(nl),
+ epsilon=as.double(eps),
+ dist2=as.double(numeric(np * nl)),
+ PACKAGE = "spatstat.utils")
+ d <- matrix(sqrt(temp$dist2), nrow=np, ncol=nl)
+ if(mintype == 2) {
+ min.which <- apply(d, 1, which.min)
+ min.d <- d[cbind(1:np, min.which)]
+ } else if (mintype == 1) {
+ min.d <- apply(d, 1, min)
+ }
+ })
+ ###### end switch #####
+ if(mintype==0)
+ return(if(listit) list(d=d) else d)
+ else if(mintype==1)
+ return(list(d=d, min.d=min.d))
+ else if(mintype==2)
+ return(list(d=d, min.d=min.d, min.which=min.which))
+}
+
+# (distance to) nearest segment
+
+distppllmin <- function(p, l, big=NULL) {
+ np <- nrow(p)
+ nl <- nrow(l)
+ if(is.null(big)) {
+ xdif <- diff(range(c(p[,1],l[, c(1,3)])))
+ ydif <- diff(range(c(p[,2],l[, c(2,4)])))
+ big <- 2 * (xdif^2 + ydif^2)
+ }
+ z <- NNdist2segments(p[,1], p[,2],
+ l[,1], l[,2], l[,3], l[,4],
+ big)
+ return(list(min.d=sqrt(z$dist2), min.which=z$index))
+}
+
+NNdist2segments <- function(xp, yp, x0, y0, x1, y1, bigvalue) {
+ np <- length(xp)
+ ns <- length(x0)
+ dist2 <- rep(bigvalue, np)
+ z <- .C(C_nndist2segs,
+ xp=as.double(xp),
+ yp=as.double(yp),
+ npoints=as.integer(np),
+ x0=as.double(x0),
+ y0=as.double(y0),
+ x1=as.double(x1),
+ y1=as.double(y1),
+ nsegments=as.integer(ns),
+ epsilon=as.double(.Machine$double.eps),
+ dist2=as.double(dist2),
+ index=as.integer(integer(np)),
+ PACKAGE = "spatstat.utils")
+ return(list(dist2=z$dist2,
+ index=z$index + 1L))
+}
diff --git a/inst/doc/packagesizes.txt b/inst/doc/packagesizes.txt
new file mode 100644
index 0000000..0fe47e7
--- /dev/null
+++ b/inst/doc/packagesizes.txt
@@ -0,0 +1,4 @@
+date version nhelpfiles nobjects ndatasets Rlines srclines
+"2017-02-23" "1.3-3" 24 159 0 3038 1714
+"2017-03-21" "1.4-1" 24 159 0 3048 1714
+"2017-07-27" "1.7-0" 27 163 0 3190 1919
diff --git a/man/articlebeforenumber.Rd b/man/articlebeforenumber.Rd
new file mode 100644
index 0000000..ad746e8
--- /dev/null
+++ b/man/articlebeforenumber.Rd
@@ -0,0 +1,39 @@
+\name{articlebeforenumber}
+\alias{articlebeforenumber}
+\title{
+ Indefinite Article Preceding A Number
+}
+\description{
+ Determines the indefinite article (\emph{an} or \code{a}) which should
+ precede a given number, if the number is read out in English.
+}
+\usage{
+articlebeforenumber(k)
+}
+\arguments{
+ \item{k}{A single integer.}
+}
+\details{
+ This function applies the rule that,
+ if the English word for the number \code{k} begins with a vowel, then
+ it should be preceded by \code{an}, and otherwise by \code{a}.
+}
+\value{
+ One of the character strings \code{"an"} or \code{"a"}.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{ordinal}}
+}
+\examples{
+ f <- function(k) cat(paste(articlebeforenumber(k),
+ paste0(k, "-fold"),
+ "increase\n"))
+ f(8)
+ f(18)
+ f(28)
+}
+\keyword{manip}
+\keyword{utilities}
diff --git a/man/cat.factor.Rd b/man/cat.factor.Rd
new file mode 100644
index 0000000..bdae941
--- /dev/null
+++ b/man/cat.factor.Rd
@@ -0,0 +1,44 @@
+\name{cat.factor}
+\alias{cat.factor}
+\title{
+ Combine Several Factors
+}
+\description{
+ Combine (concatenate) several factor objects, to produce a factor.
+}
+\usage{
+cat.factor(\dots)
+}
+\arguments{
+ \item{\dots}{
+ Any number of arguments. Each argument should be a factor,
+ or will be converted to a factor.
+ }
+}
+\details{
+ The arguments \code{\dots} are concatenated as they would be
+ using \code{\link[base]{c}()} or \code{\link[base]{cat}()},
+ except that factor levels are retained
+ and merged correctly. See the Examples.
+}
+\value{
+ A factor, whose length is the sum of the lengths of all
+ arguments. The levels of the resulting factor are the union of the
+ levels of the arguments.
+}
+\author{
+ \rolf.
+}
+\seealso{
+ \code{\link[base]{c}}.
+}
+\examples{
+ f <- factor(letters[1:3])
+ g <- factor(letters[3:5])
+ f
+ g
+ cat(f,g)
+ c(f,g)
+ cat.factor(f, g)
+}
+\keyword{manip}
diff --git a/man/check.1.integer.Rd b/man/check.1.integer.Rd
new file mode 100644
index 0000000..9fc7668
--- /dev/null
+++ b/man/check.1.integer.Rd
@@ -0,0 +1,59 @@
+\name{check.1.integer}
+\alias{check.1.integer}
+\alias{check.1.real}
+\alias{check.1.string}
+\title{
+ Check Argument Type and Length
+}
+\description{
+ These utility functions check whether a given argument is a single value
+ of the required type.
+}
+\usage{
+check.1.real(x, context = "", fatal = TRUE)
+check.1.integer(x, context = "", fatal = TRUE)
+check.1.string(x, context = "", fatal = TRUE)
+}
+\arguments{
+ \item{x}{
+ The argument to be checked.
+ }
+ \item{context}{
+ Optional string specifying the context in which the argument
+ is checked.
+ }
+ \item{fatal}{
+ Logical value indicating what to do if \code{x} is not of the
+ required type.
+ }
+}
+\details{
+ These functions check whether the argument \code{x} is a single
+ atomic value of type \code{numeric}, \code{integer} or
+ \code{character}.
+
+ If \code{x} does have the required length and type, the result
+ of the function is the logical value \code{TRUE}.
+
+ If \code{x} does not have the required length and type,
+ then if \code{fatal=TRUE} (the default) an error occurs,
+ while if \code{fatal=FALSE} a warning is issued and the
+ function returns the value \code{FALSE}.
+}
+\value{
+ A logical value (or an error may occur).
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{check.named.vector}}
+}
+\examples{
+ x <- pi
+ check.1.real(x)
+ check.1.integer(x, fatal=FALSE)
+ check.1.string(x, fatal=FALSE)
+}
+\keyword{classes}
+\keyword{error}
diff --git a/man/check.named.vector.Rd b/man/check.named.vector.Rd
new file mode 100644
index 0000000..15c5d06
--- /dev/null
+++ b/man/check.named.vector.Rd
@@ -0,0 +1,96 @@
+\name{check.named.vector}
+\alias{check.named.vector}
+\alias{check.named.list}
+\alias{check.named.thing}
+\title{
+ Check Whether Object Has Required Components
+}
+\description{
+ These functions check whether the object \code{x} has
+ components with the required names, and does not have
+ any unexpected components.
+}
+\usage{
+check.named.vector(x, nam, context = "", namopt = character(0),
+ onError = c("fatal", "null"))
+
+check.named.list(x, nam, context = "", namopt = character(0),
+ onError = c("fatal", "null"))
+
+check.named.thing(x, nam, namopt = character(0),
+ xtitle = NULL, valid = TRUE, type = "object",
+ context = "", fatal = TRUE)
+}
+\arguments{
+ \item{x}{
+ The object to be checked.
+ }
+ \item{nam}{
+ Vector of character strings giving the names of all the
+ components which must be present.
+ }
+ \item{namopt}{
+ Vector of character strings giving the names of components
+ which may optionally be present.
+ }
+ \item{context}{
+ Character string describing the context in which \code{x} is being checked.
+ }
+ \item{xtitle}{
+ Optional character string to be used when referring to \code{x}.
+ }
+ \item{valid}{
+ Logical value indicating whether \code{x} belongs to the required
+ class of objects.
+ }
+ \item{type}{
+ Character string describing the required class of objects.
+ }
+ \item{onError}{
+ Character string indicating what to do if \code{x} fails the checks.
+ }
+ \item{fatal}{
+ Logical value indicating what to do if \code{x} fails the checks.
+ If \code{fatal=TRUE} (the default), an error occurs.
+ }
+}
+\details{
+ \code{check.named.thing} checks whether \code{x} has all the
+ required components, in the sense that \code{names(x)}
+ includes all the names in \code{nam},
+ and that every entry in \code{names(x)} belongs to either \code{nam} or
+ \code{namopt}. If all these checks are true, the result is a
+ zero-length character vector. Otherwise, if \code{fatal=TRUE} (the
+ default), an error occurs; otherwise the result is a character
+ vector describing the checks which failed.
+
+ \code{check.named.vector} checks whether \code{x} is a numeric vector
+ and \code{check.named.list} checks whether \code{x} is a list.
+ They then call \code{check.named.thing} to check whether all the
+ required components are present. If all these checks are true, the result is
+ a reordered version of \code{x} in which all the compulsory entries
+ appear first. Otherwise, if \code{onError="fatal"} (the default)
+ an error occurs; otherwise the result is \code{NULL}.
+}
+\value{
+ \code{check.named.vector} returns a numeric vector or \code{NULL}.
+
+ \code{check.named.list} returns a list or \code{NULL}.
+
+ \code{check.named.thing} returns a character vector.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{check.1.integer}}
+}
+\examples{
+ z <- list(a=1, b=2, e=42)
+ check.named.list(z, c("a", "b"), namopt=c("c", "d", "e"))
+ check.named.thing(z, c("a", "b"), namopt=c("c", "d", "e"))
+ zz <- unlist(z)
+ check.named.vector(zz, c("a", "b"), namopt=c("c", "d", "e"))
+ check.named.thing(z, c("b", "c"), namopt=c("d", "e"), fatal=FALSE)
+}
+\keyword{error}
diff --git a/man/check.nmatrix.Rd b/man/check.nmatrix.Rd
new file mode 100644
index 0000000..f483068
--- /dev/null
+++ b/man/check.nmatrix.Rd
@@ -0,0 +1,79 @@
+\name{check.nmatrix}
+\alias{check.nmatrix}
+\title{
+ Check for Numeric Matrix with Correct Dimensions
+}
+\description{
+ This is a programmer's utility function
+ to check whether the argument is a numeric
+ vector of the correct length.
+}
+\usage{
+check.nmatrix(m, npoints = NULL, fatal = TRUE, things = "data points",
+ naok = FALSE, squarematrix = TRUE, matchto = "nrow", warn = FALSE)
+}
+\arguments{
+ \item{m}{
+ The argument to be checked.
+ }
+ \item{npoints}{
+ The required number of rows and/or columns for the matrix \code{m}.
+ }
+ \item{fatal}{
+ Logical value indicating whether to stop with an error message
+ if \code{m} does not satisfy all requirements.
+ }
+ \item{things}{
+ Character string describing what the rows/columns of \code{m} should
+ correspond to.
+ }
+ \item{naok}{
+ Logical value indicating whether \code{NA} values are permitted.
+ }
+ \item{squarematrix}{
+ Logical value indicating whether \code{m} must be a square matrix.
+ }
+ \item{matchto}{
+ Character string (either \code{"nrow"} or \code{"ncol"})
+ indicating whether it is the rows or the columns
+ of \code{m} which must correspond to \code{npoints}.
+ }
+ \item{warn}{
+ Logical value indicating whether to issue a warning
+ if \code{v} does not satisfy all requirements.
+ }
+}
+\details{
+ This programmer's utility function checks whether \code{m} is a numeric matrix
+ of the correct dimensions, and checks for \code{NA} values.
+ If \code{matchto="nrow"} (the default) then
+ the number of rows of \code{m} must be equal to \code{npoints}.
+ If \code{matchto="ncol"} then the number of columns of \code{m}
+ must be equal to \code{npoints}. If \code{squarematrix=TRUE} (the
+ default) then the numbers of rows and columns must be equal.
+ If \code{naok = FALSE} (the default) then the entries of \code{m}
+ must not include \code{NA}.
+
+ If these requirements are all satisfied, the result is the logical
+ value \code{TRUE}.
+
+ If not, then if \code{fatal=TRUE} (the default), an error occurs;
+ if \code{fatal=FALSE}, the result is the logical value \code{FALSE}
+ with an attribute describing the requirement that was not satisfied.
+}
+\value{
+ A logical value indicating whether all the requirements were
+ satisfied.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{check.nvector}}
+}
+\examples{
+ z <- matrix(1:16, 4, 4)
+ check.nmatrix(z, 4)
+}
+\keyword{error}
+\keyword{utilities}
diff --git a/man/check.nvector.Rd b/man/check.nvector.Rd
new file mode 100644
index 0000000..2c02879
--- /dev/null
+++ b/man/check.nvector.Rd
@@ -0,0 +1,78 @@
+\name{check.nvector}
+\alias{check.nvector}
+\title{
+ Check For Numeric Vector With Correct Length
+}
+\description{
+ This is a programmer's utility function
+ to check whether the argument is a numeric
+ vector of the correct length.
+}
+\usage{
+ check.nvector(v, npoints = NULL, fatal = TRUE, things = "data points",
+ naok = FALSE, warn = FALSE, vname, oneok = FALSE)
+}
+\arguments{
+ \item{v}{
+ The argument to be checked.
+ }
+ \item{npoints}{
+ The required length of \code{v}.
+ }
+ \item{fatal}{
+ Logical value indicating whether to stop with an error message
+ if \code{v} does not satisfy all requirements.
+ }
+ \item{things}{
+ Character string describing what the entries of \code{v} should
+ correspond to.
+ }
+ \item{naok}{
+ Logical value indicating whether \code{NA} values are permitted.
+ }
+ \item{warn}{
+ Logical value indicating whether to issue a warning
+ if \code{v} does not satisfy all requirements.
+ }
+ \item{vname}{
+ Character string giving the name of \code{v} to be used in messages.
+ }
+ \item{oneok}{
+ Logical value indicating whether \code{v} is permitted to have
+ length 1.
+ }
+}
+\details{
+ This function checks whether \code{v} is a numeric vector with
+ length equal to \code{npoints} (or length equal to 1 if
+ \code{oneok=TRUE}), not containing any \code{NA} values (unless
+ \code{naok=TRUE}).
+
+ If these requirements are all satisfied, the result is the logical
+ value \code{TRUE}.
+
+ If not, then if \code{fatal=TRUE} (the default), an error occurs;
+ if \code{fatal=FALSE}, the result is the logical value \code{FALSE}
+ with an attribute describing the requirement that was not satisfied.
+}
+\value{
+ A logical value indicating whether all the requirements were
+ satisfied. If \code{FALSE}, then this value has an attribute
+ \code{"whinge"}, a character string describing the requirements that
+ were not satisfied.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{check.nmatrix}},
+ \code{\link{check.1.real}}, \code{\link{check.named.vector}}.
+}
+\examples{
+ z <- 1:10
+ check.nvector(z, 5, fatal=FALSE)
+ y <- 42
+ check.nvector(y, 5, fatal=FALSE, oneok=TRUE)
+}
+\keyword{error}
+\keyword{utilities}
diff --git a/man/check.range.Rd b/man/check.range.Rd
new file mode 100644
index 0000000..94a20b8
--- /dev/null
+++ b/man/check.range.Rd
@@ -0,0 +1,97 @@
+\name{check.range}
+\alias{check.range}
+\alias{check.in.range}
+\alias{inside.range}
+\alias{intersect.ranges}
+\alias{prange}
+\title{
+ Utilities for Ranges of Values
+}
+\description{
+ These simple functions handle an interval or range of numerical
+ values. \code{check.range(r)} checks whether \code{r} specifies a
+ range of values, that is, whether \code{r} is a vector of length 2
+ with \code{r[1] <= r[2]}. \code{intersect.ranges(r, s)} finds the intersection
+ of two ranges \code{r} and \code{s}. \code{inside.range(x, r)} returns
+ a logical vector containing \code{TRUE} if the corresponding entry of
+ \code{x} falls inside the range \code{r}, and \code{FALSE} if it does
+ not. \code{check.in.range(x, r)} checks whether a single number
+ \code{x} falls inside the specified range \code{r}.
+ Finally \code{prange(r)} produces a character string that represents
+ the range \code{r}.
+}
+\usage{
+check.range(r, fatal = TRUE)
+
+check.in.range(x, r, fatal = TRUE)
+
+inside.range(x, r)
+
+intersect.ranges(r, s, fatal = TRUE)
+
+prange(r)
+}
+\arguments{
+ \item{r}{
+ A numeric vector of length 2 specifying the endpoints of a range of values.
+ }
+ \item{x}{
+ Numeric vector of data.
+ }
+ \item{s}{
+ A numeric vector of length 2 specifying the endpoints of a range of values.
+ }
+ \item{fatal}{
+ Logical value indicating whether to stop with an error message
+ if the data do not pass the check.
+ }
+}
+\details{
+ \code{check.range} checks whether \code{r} specifies a
+ range of values, that is, whether \code{r} is a vector of length 2
+ with \code{r[1] <= r[2]}. If so, the result is \code{TRUE}. If not,
+ then if \code{fatal=TRUE}, an error occurs, while if
+ \code{fatal=FALSE} the result is \code{FALSE}.
+
+ \code{intersect.ranges(r, s)} finds the intersection
+ of two ranges \code{r} and \code{s}. If the intersection is non-empty,
+ the result is a numeric vector of length 2. If the intersection is empty,
+ then if \code{fatal=TRUE}, an error occurs, while if
+ \code{fatal=FALSE} the result is \code{NULL}.
+
+ \code{inside.range(x, r)} returns
+ a logical vector containing \code{TRUE} if the corresponding entry of
+ \code{x} falls inside the range \code{r}, and \code{FALSE} if it does
+ not.
+
+ \code{check.in.range(x, r)} checks whether a single number
+ \code{x} falls inside the specified range \code{r}.
+ If so, the result is \code{TRUE}. If not, then if \code{fatal=TRUE},
+ an error occurs, while if
+ \code{fatal=FALSE} the result is \code{FALSE}.
+
+ Finally \code{prange(r)} produces a character string that represents
+ the range \code{r}.
+}
+\value{
+ The result of \code{check.range}, \code{check.in.range}
+ and \code{inside.range}, is
+ a logical value or logical vector. The result of
+ \code{intersect.ranges} is a numerical vector of length 2, or \code{NULL}.
+ The result of \code{prange} is a character string.
+}
+\author{
+ \adrian
+}
+\examples{
+ rr <- c(0, 2)
+ ss <- c(1, 3)
+ x <- seq(0.5, 3.5, by=1)
+ check.range(rr)
+ check.range(42, fatal=FALSE)
+ inside.range(x, rr)
+ intersect.ranges(rr, ss)
+ prange(rr)
+}
+\keyword{programming}
+\keyword{utilities}
diff --git a/man/commasep.Rd b/man/commasep.Rd
new file mode 100644
index 0000000..51e21ff
--- /dev/null
+++ b/man/commasep.Rd
@@ -0,0 +1,38 @@
+\name{commasep}
+\alias{commasep}
+\title{
+ List of Items Separated By Commas
+}
+\description{
+ Convert the elements of a vector into character strings
+ and paste them together, separated by commas.
+}
+\usage{
+commasep(x, join = " and ", flatten = TRUE)
+}
+\arguments{
+ \item{x}{
+ Vector of items in the list.
+ }
+ \item{join}{
+ The string to be used to separate the last two items in the list.
+ }
+ \item{flatten}{
+ Logical value indicating whether to return a single character string
+ (\code{flatten=TRUE}, the default) or a list (\code{flatten=FALSE}).
+ }
+}
+\value{
+ A character string (if \code{flatten=TRUE}, the default)
+ or a list of character strings.
+}
+\author{
+ \adrian.
+}
+\examples{
+ commasep(letters[1:4])
+ y <- commasep(sQuote(letters[1:4]))
+ cat(y, fill=TRUE)
+}
+\keyword{utilities}
+\keyword{manip}
diff --git a/man/do.call.matched.Rd b/man/do.call.matched.Rd
new file mode 100644
index 0000000..26cce40
--- /dev/null
+++ b/man/do.call.matched.Rd
@@ -0,0 +1,83 @@
+\name{do.call.matched}
+\alias{do.call.matched}
+\title{
+ Call a Function, Passing Only Recognised Arguments
+}
+\description{
+ Call a specified function, using only those arguments which are
+ known to be acceptable to the function.
+}
+\usage{
+ do.call.matched(fun, arglist, funargs, extrargs = NULL,
+ matchfirst = FALSE, sieve = FALSE, skipargs = NULL)
+}
+\arguments{
+ \item{fun}{
+ A function, or a character string giving the name of a function,
+ to be called.
+ }
+ \item{arglist}{
+ A named list of arguments.
+ }
+ \item{funargs}{
+ Character vector giving the names of arguments that are recognised
+ by \code{fun}. Defaults to the names of the formal arguments of
+ \code{fun}.
+ }
+ \item{extrargs}{
+ Optional. Character vector giving the names of additional arguments
+ that can be handled by \code{fun}.
+ }
+ \item{skipargs}{
+ Optional. Character vector giving the names of arguments which should
+ \bold{not} be passed to \code{fun}.
+ }
+ \item{matchfirst}{
+ Logical value indicating whether the first entry
+ of \code{arglist} is permitted to have an empty name
+ and should be matched to the first argument of \code{fun}.
+ }
+ \item{sieve}{
+ Logical value indicating whether to return the
+ un-used arguments as well as the result of the function call.
+ See Details.
+ }
+}
+\details{
+ This function is a wrapper for \code{\link[base]{do.call}}
+ which avoids passing arguments that are unrecognised by \code{fun}.
+
+ In the simplest case \code{do.call.matched(fun, arglist)}
+ is like \code{do.call(fun, arglist)}, except that
+ entries of \code{arglist} which do not match any formal
+ argument of \code{fun} are removed.
+ Extra argument names can be permitted using \code{extrargs},
+ and argument names can be forbidden using \code{skipargs}.
+}
+\value{
+ If \code{sieve=FALSE} (the default), the result is the
+ return value from \code{fun}.
+
+ If \code{sieve=TRUE}, the result is a list with entries
+ \code{result} (the return value from \code{fun}) and
+ \code{otherargs} (a list of the arguments that were not passed
+ to \code{fun}).
+}
+\author{
+ \adrian
+}
+\seealso{
+ \code{\link{resolve.defaults}},
+ \code{\link{do.call.without}}.
+
+ \code{\link[base]{do.call}}
+}
+\examples{
+ f <- function(x=0,y=0, ...) { paste(x, y, ..., sep=", ") }
+ f()
+ do.call.matched(f, list(y=2))
+ do.call.matched(f, list(y=2, z=5), extrargs="z")
+ do.call.matched(f, list(y=2, z=5), extrargs="z", skipargs="y")
+}
+\keyword{programming}
+\keyword{utilities}
diff --git a/man/do.call.without.Rd b/man/do.call.without.Rd
new file mode 100644
index 0000000..2b24cb4
--- /dev/null
+++ b/man/do.call.without.Rd
@@ -0,0 +1,47 @@
+\name{do.call.without}
+\alias{do.call.without}
+\title{
+ Call a Function, Omitting Certain Arguments
+}
+\description{
+ Call a specified function, omitting some arguments which are
+ inappropriate to the function.
+}
+\usage{
+ do.call.without(fun, \dots, avoid)
+}
+\arguments{
+ \item{fun}{
+ The function to be called. A function name, a
+ character string giving the name of the function,
+ or an expression that yields a function.
+ }
+ \item{\dots}{
+ Any number of arguments.
+ }
+ \item{avoid}{
+ Vector of character strings, giving the names of arguments that should
+ \emph{not} be passed to \code{fun}.
+ }
+}
+\details{
+ This is a simple mechanism for preventing some arguments from being
+ passed in a function call. The arguments \code{\dots} are collected in
+ a list. A argument is omitted if its name exactly matches
+ one of the strings in \code{avoid}.
+}
+\value{
+ The return value of \code{fun}.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{do.call.matched}} for a more complicated and flexible call.
+}
+\examples{
+ do.call.without(paste, 1, 2, z=3, w=4, avoid="z")
+}
+\keyword{programming}
+\keyword{utilities}
+
diff --git a/man/expand.polynom.Rd b/man/expand.polynom.Rd
new file mode 100644
index 0000000..9dadff0
--- /dev/null
+++ b/man/expand.polynom.Rd
@@ -0,0 +1,57 @@
+\name{expand.polynom}
+\alias{expand.polynom}
+\alias{sympoly}
+\title{
+ Expand Symbolic Polynomials in a Formula
+}
+\description{
+ Create a formula representing a polynomial,
+ or expand polynomials in an existing formula.
+}
+\usage{
+ expand.polynom(f)
+ sympoly(x, y, n)
+}
+\arguments{
+ \item{f}{
+ A formula.
+ }
+ \item{x,y}{
+ Variable names.
+ }
+ \item{n}{
+ Integer specifying the degree of the polynomial.
+ (If \code{n} is missing, \code{y} will be interpreted as the degree.)
+ }
+}
+\details{
+ These functions expand a polynomial into its homogeneous terms
+ and return a model formula.
+
+ \code{sympoly(x, n)} creates a formula whose right-hand side represents the
+ polynomial of degree \code{n} in the variable \code{x}. Each
+ homogeneous term \code{x^k} is a separate term in the formula.
+
+ \code{sympoly(x, y, n)} creates a formula representing the
+ polynomial of degree \code{n} in the two variables \code{x} and
+ \code{y}.
+
+ If \code{f} is a formula containing a term of the form
+ \code{polynom(\dots)} then \code{expand.polynom(f)} replaces this term
+ by its expansion as a sum of homogeneous terms.
+}
+\value{
+ A \code{formula}.
+}
+\author{
+ \spatstatAuthors.
+}
+\seealso{
+ \code{\link[spatstat]{polynom}}
+}
+\examples{
+ sympoly(A, 4)
+ sympoly(A, B, 3)
+ expand.polynom(U ~ A + polynom(B, 2))
+}
+\keyword{utilities}
diff --git a/man/ifelseAB.Rd b/man/ifelseAB.Rd
new file mode 100644
index 0000000..c5f3668
--- /dev/null
+++ b/man/ifelseAB.Rd
@@ -0,0 +1,76 @@
+\name{ifelseAB}
+\alias{ifelse0NA}
+\alias{ifelse1NA}
+\alias{ifelseAB}
+\alias{ifelseAX}
+\alias{ifelseXB}
+\alias{ifelseXY}
+\alias{ifelseNegPos}
+\title{
+ Conditional Selection
+}
+\description{
+ These low-level functions provide faster alternatives to
+ some uses of \code{ifelse}.
+}
+\usage{
+ifelseAB(test, a, b)
+ifelseAX(test, a, x)
+ifelseXB(test, x, b)
+ifelseXY(test, x, y)
+ifelseNegPos(test, x)
+ifelse0NA(test)
+ifelse1NA(test)
+}
+\arguments{
+ \item{test}{A logical vector.}
+ \item{a}{A single atomic value.}
+ \item{b}{A single atomic value.}
+ \item{x}{A vector of values, of the same length as \code{test}.}
+ \item{y}{A vector of values, of the same length as \code{test}.}
+}
+\details{
+ These low-level functions provide faster alternatives to
+ some uses of \code{\link[base]{ifelse}}. They were developed by
+ trial-and-error comparison of computation times of different expressions.
+
+ \code{ifelse0NA(test)} is equivalent to \code{ifelse(test, 0, NA)}.
+
+ \code{ifelse1NA(test)} is equivalent to \code{ifelse(test, 1, NA)}.
+
+ \code{ifelseAB(test, a, b)} is equivalent to \code{ifelse(test, a, b)}
+ where \code{a} and \code{b} must be single values.
+
+ \code{ifelseAX(test, a, x)} is equivalent to \code{ifelse(test, a, x)}
+ where \code{a} must be a single value, and \code{x} a vector of the
+ same length as \code{test}.
+
+ \code{ifelseXB(test, x, b)} is equivalent to \code{ifelse(test, x, b)}
+ where \code{b} must be a single value, and \code{x} a vector of the
+ same length as \code{test}.
+
+ \code{ifelseXY(test, x, y)} is equivalent to \code{ifelse(test, x, y)}
+ where \code{x} and \code{y} must be vectors of the
+ same length as \code{test}.
+
+ \code{ifelseNegPos(test, x)} is equivalent to \code{ifelse(test, x, -x)}
+ where \code{x} must be a vector of the same length as \code{test}.
+}
+\value{
+ A vector of the same length as \code{test} containing values of the
+ same type as \code{a,b,x,y}.
+}
+\author{
+ \spatstatAuthors.
+}
+\seealso{
+ \code{\link[base]{ifelse}}
+}
+\examples{
+ x <- runif(4e5)
+ u <- (x < 0.5)
+ system.time(ifelse(u, 2, x))
+ system.time(ifelseAX(u, 2, x))
+}
+\keyword{manip}
+\keyword{utilities}
diff --git a/man/macros/defns.Rd b/man/macros/defns.Rd
new file mode 100644
index 0000000..466babb
--- /dev/null
+++ b/man/macros/defns.Rd
@@ -0,0 +1,19 @@
+%% macro definitions for spatstat man pages
+\newcommand{\adrian}{Adrian Baddeley \email{Adrian.Baddeley at curtin.edu.au}}
+\newcommand{\rolf}{Rolf Turner \email{r.turner at auckland.ac.nz}}
+\newcommand{\ege}{Ege Rubak \email{rubak at math.aau.dk}}
+\newcommand{\spatstatAuthors}{\adrian, \rolf and \ege}
+% Names with accents
+\newcommand{\Bogsted}{\ifelse{latex}{\out{B\o gsted}}{Bogsted}}
+\newcommand{\Cramer}{\ifelse{latex}{\out{Cram\'er}}{Cramer}}
+\newcommand{\Hogmander}{\ifelse{latex}{\out{H{\"o}gmander}}{Hogmander}}
+\newcommand{\Jyvaskyla}{\ifelse{latex}{\out{Jyv\"askyl\"a}}{Jyvaskyla}}
+\newcommand{\Matern}{\ifelse{latex}{\out{Mat\'ern}}{Matern}}
+\newcommand{\Moller}{\ifelse{latex}{\out{M\o ller}}{Moller}}
+\newcommand{\Oehlschlaegel}{\ifelse{latex}{\out{Oehlschl\"{a}gel}}{Oehlschlaegel}}
+\newcommand{\Prokesova}{\ifelse{latex}{\out{Proke\u{s}ov{\'{a}}}}{Prokesova}}
+\newcommand{\Sarkka}{\ifelse{latex}{\out{S\"{a}rkk\"{a}}}{Sarkka}}
+%% List of all Gibbs interactions
+\newcommand{\GibbsInteractionsList}{\code{\link{AreaInter}}, \code{\link{BadGey}}, \code{\link{Concom}}, \code{\link{DiggleGatesStibbard}}, \code{\link{DiggleGratton}}, \code{\link{Fiksel}}, \code{\link{Geyer}}, \code{\link{Hardcore}}, \code{\link{Hybrid}}, \code{\link{LennardJones}}, \code{\link{MultiStrauss}}, \code{\link{MultiStraussHard}}, \code{\link{OrdThresh}}, \code{\link{Ord}}, \code{\link{Pairwise}}, \code{\link{PairPiece}}, \code{\link{Penttinen}}, \code{\link{Poisson}}, \code [...]
+%% List of interactions recognised by RMH code
+\newcommand{\rmhInteractionsList}{\code{\link{AreaInter}}, \code{\link{BadGey}}, \code{\link{DiggleGatesStibbard}}, \code{\link{DiggleGratton}}, \code{\link{Fiksel}}, \code{\link{Geyer}}, \code{\link{Hardcore}}, \code{\link{Hybrid}}, \code{\link{LennardJones}}, \code{\link{MultiStrauss}}, \code{\link{MultiStraussHard}}, \code{\link{PairPiece}}, \code{\link{Penttinen}}, \code{\link{Poisson}}, \code{\link{Softcore}}, \code{\link{Strauss}}, \code{\link{StraussHard}} and \code{\link{Triplets}}}
diff --git a/man/methods.xypolygon.Rd b/man/methods.xypolygon.Rd
new file mode 100644
index 0000000..0ea9b6a
--- /dev/null
+++ b/man/methods.xypolygon.Rd
@@ -0,0 +1,127 @@
+\name{methods.xypolygon}
+\alias{methods.xypolygon} %DoNotExport
+\alias{verify.xypolygon}
+\alias{is.hole.xypolygon}
+\alias{Area.xypolygon}
+\alias{bdrylength.xypolygon}
+\alias{reverse.xypolygon}
+\alias{overlap.xypolygon}
+\alias{simplify.xypolygon}
+\alias{inside.xypolygon}
+\title{
+ Calculations for Polygons in the Plane
+}
+\description{
+ Compute the area or boundary length of a polygon,
+ determine whether a point falls inside a polygon,
+ compute the area of overlap between two polygons,
+ and related tasks.
+}
+\usage{
+verify.xypolygon(p, fatal = TRUE)
+is.hole.xypolygon(polly)
+Area.xypolygon(polly)
+bdrylength.xypolygon(polly)
+reverse.xypolygon(p, adjust=FALSE)
+overlap.xypolygon(P, Q)
+simplify.xypolygon(p, dmin)
+inside.xypolygon(pts, polly, test01, method)
+}
+\arguments{
+ \item{p,polly,P,Q}{
+ Data representing a polygon. See Details.
+ }
+ \item{dmin}{
+ Single numeric value giving the minimum permissible
+ length of an edge in the simplified polygon.
+ }
+ \item{fatal}{
+ Logical value indicating whether failure is a fatal error.
+ }
+ \item{pts}{
+ Coordinates of points to be tested.
+ A named list with entries \code{x,y} which are numeric vectors
+ of coordinates.
+ }
+ \item{adjust}{
+ Logical value indicating whether internal data should be adjusted.
+ See Details.
+ }
+ \item{test01,method}{
+ For developer use only.
+ }
+}
+\details{
+ In the \pkg{spatstat} family of packages, a polygon in the
+ Euclidean plane is represented as a named list
+ with the following entries:
+ \describe{
+ \item{\code{x},\code{y}}{
+ Numeric vectors giving the coordinates of the vertices.
+ The vertices should be traversed in anti-clockwise order
+ (unless the polygon is a hole, when they should be traversed
+ in clockwise order) and the last vertex should \bold{not} repeat
+ the first vertex.
+ }
+ \item{hole}{
+ Optional. A logical value indicating whether the
+ polygon is a hole.
+ }
+ \item{area}{
+ Optional. Single numeric value giving the area of the polygon
+ (negative if it is a hole).
+ }
+ }
+ The function \code{verify.xypolygon} checks whether its argument
+ satisfies this format. If so, it returns \code{TRUE}; if not,
+ it returns \code{FALSE} or (if \code{fatal=TRUE}) generates a fatal error.
+
+ The other functions listed here perform basic calculations for
+ polygons using elementary Cartesian analytic geometry in \R.
+
+ \code{is.hole.xypolygon} determines whether a polygon is a hole or not.
+
+ \code{Area.xypolygon} computes the area of the polygon using the
+ discrete Green's formula.
+
+ \code{bdrylength.xypolygon} calculates the total length of edges
+ of the polygon.
+
+ \code{reverse.xypolygon} reverses the order of the
+ coordinate vectors \code{x} and \code{y}. If \code{adjust=TRUE},
+ the other entries \code{hole} and \code{area} will be adjusted as well.
+
+ \code{overlap.xypolygon} computes the area of overlap between two
+ polygons using the discrete Green's formula. It is slow compared
+ to the code in the \pkg{polyclip} package.
+
+ \code{simplify.xypolygon} removes vertices of the polygon until
+ every edge is longer than \code{dmin}.
+
+ \code{inside.xypolygon(pts, polly)} determines whether each point
+ in \code{pts} lies inside the polygon \code{polly} and returns a
+ logical vector.
+}
+\value{
+ \code{verify.xypolygon} and
+ \code{is.hole.xypolygon} return a single logical value.
+
+ \code{inside.xypolygon} returns a logical vector.
+
+ \code{Area.xypolygon}, \code{bdrylength.xypolygon}
+ and \code{overlap.xypolygon}
+ return a single numeric value.
+
+ \code{reverse.xypolygon} and \code{simplify.xypolygon}
+ return another polygon object.
+}
+\author{
+ \adrian.
+}
+\examples{
+ p <- list(x=c(0,1,4,2), y=c(0,0,2,3))
+ is.hole.xypolygon(p)
+ Area.xypolygon(p)
+ bdrylength.xypolygon(p)
+}
+\keyword{math}
diff --git a/man/optimizeWithTrace.Rd b/man/optimizeWithTrace.Rd
new file mode 100644
index 0000000..1fa382d
--- /dev/null
+++ b/man/optimizeWithTrace.Rd
@@ -0,0 +1,65 @@
+\name{optimizeWithTrace}
+\alias{optimizeWithTrace}
+\title{
+ One Dimensional Optimization with Tracing
+}
+\description{
+ Find the minimum or maximum of a function over an interval
+ of real numbers, keeping track of the function arguments and
+ function values that were evaluated.
+}
+\usage{
+optimizeWithTrace(f, interval, \dots,
+ lower = min(interval), upper = max(interval))
+}
+\arguments{
+ \item{f}{
+ The function to be minimized or maximized.
+ }
+ \item{interval}{
+ Numeric vector of length 2 containing the end-points of the interval
+ to be searched.
+ }
+ \item{lower, upper}{
+ The lower and upper endpoints of the interval to be searched.
+ }
+ \item{\dots}{
+ Other arguments passed to \code{\link[stats]{optimize}},
+ including arguments to the function \code{f}.
+ }
+}
+\details{
+ This is a simple wrapper for the optimization routine
+ \code{\link[stats]{optimize}}. The function \code{f} will be
+ optimized by computing its value at several locations in the interval,
+ as described in the help for \code{\link[stats]{optimize}}.
+ This wrapper function stores the locations and
+ resulting function values, and returns them along with the
+ result of the optimization.
+}
+\value{
+ A list with components
+ \itemize{
+ \item \code{minimum} (or \code{maximum}), the location in the
+ search interval which yielded the optimum value;
+ \item \code{objective}, the value of the function at this
+ location;
+ \item \code{x}, the sequence of locations in the interval that
+ were considered (in the order considered);
+ \item \code{y}, the function values corresponding to \code{x}.
+ }
+}
+\author{
+ \spatstatAuthors.
+}
+\seealso{
+ \code{\link{optimize}}
+}
+\examples{
+ f <- function (x, a) (x - a)^2
+ result <- optimizeWithTrace(f, c(0, 1), tol = 0.0001, a = 1/3)
+ result
+ curve(f(x, 1/3))
+ with(result, points(x, y, pch=16))
+}
+\keyword{optimize}
diff --git a/man/ordinal.Rd b/man/ordinal.Rd
new file mode 100644
index 0000000..d30bcbf
--- /dev/null
+++ b/man/ordinal.Rd
@@ -0,0 +1,43 @@
+\name{ordinal}
+\alias{ordinal}
+\alias{ordinalsuffix}
+\title{
+ Ordinal Numbers
+}
+\description{
+ Returns the appropriate abbreviation in English for an ordinal number
+ (for example \code{ordinal(5)} is \code{"5th"}).
+}
+\usage{
+ordinal(k)
+ordinalsuffix(k)
+}
+\arguments{
+ \item{k}{An integer or vector of integers.}
+}
+\details{
+ \code{ordinal(k)} returns a character string representing the
+ \code{k}th ordinal number. \code{ordinalsuffix(k)} determines the
+ appropriate suffix.
+
+ The suffix can be either \code{"st"} (abbreviating
+ \emph{first}), \code{"nd"} (abbreviating \emph{second}),
+ \code{"rd"} (abbreviating \emph{third}) or
+ \code{"th"} (for all other ordinal numbers \code{fourth},
+ \code{fifth}, etc).
+}
+\value{
+ A character string or character vector of the same length as \code{k}.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{articlebeforenumber}}
+}
+\examples{
+ ordinal(1:7)
+ cat(paste("Happy", ordinal(21), "Birthday"), fill=TRUE)
+}
+\keyword{manip}
+\keyword{utilities}
diff --git a/man/paren.Rd b/man/paren.Rd
new file mode 100644
index 0000000..a407fd1
--- /dev/null
+++ b/man/paren.Rd
@@ -0,0 +1,45 @@
+\name{paren}
+\alias{paren}
+\alias{unparen}
+\title{
+ Add or Remove Parentheses
+}
+\description{
+ Add or remove enclosing parentheses around a string.
+}
+\usage{
+paren(x, type = "(")
+unparen(x)
+}
+\arguments{
+ \item{x}{
+ A character string, or vector of character strings.
+ }
+ \item{type}{
+ Type of parentheses: either \code{"("}, \code{"["} or \code{"{"}.
+ }
+}
+\details{
+ \code{paren(x)} adds enclosing parentheses to the beginning and end of
+ the string \code{x}.
+
+ \code{unparen(x)} removes enclosing parentheses if they are present.
+}
+\value{
+ A character string, or vector of character strings of the same length
+ as \code{x}.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link{commasep}}
+}
+\examples{
+ paren("Hello world")
+ paren(42, "[")
+ paren(letters[1:10])
+ unparen(c("(yes)", "[no]", "{42}"))
+}
+\keyword{manip}
+\keyword{utilities}
diff --git a/man/primefactors.Rd b/man/primefactors.Rd
new file mode 100644
index 0000000..ce13291
--- /dev/null
+++ b/man/primefactors.Rd
@@ -0,0 +1,87 @@
+\name{primefactors}
+\alias{primefactors}
+\alias{primesbelow}
+\alias{divisors}
+\alias{is.prime}
+\alias{relatively.prime}
+\alias{greatest.common.divisor}
+\alias{least.common.multiple}
+\title{
+ Primes, Prime Factorization, Common Divisor
+}
+\description{
+ These functions find prime numbers,
+ factorise a composite number into its prime factors,
+ determine whether a number is prime, and find the least common multiple
+ or greatest common divisor of two numbers.
+}
+\usage{
+primefactors(n, method=c("C", "interpreted"))
+divisors(n)
+is.prime(n)
+relatively.prime(n, m)
+least.common.multiple(n,m)
+greatest.common.divisor(n,m)
+primesbelow(nmax)
+}
+\arguments{
+ \item{n,m}{
+ Integers to be factorized.
+ }
+ \item{nmax}{
+ Integer. Upper limit on prime numbers to be found.
+ }
+ \item{method}{
+ Character string indicating the choice of algorithm.
+ (Developer use only.)
+ }
+}
+\details{
+ \code{is.prime(n)} returns \code{TRUE} if \code{n} is a prime number,
+ and \code{FALSE} otherwise.
+
+ \code{primefactors(n)} factorises the integer \code{n}
+ into its prime number factors, and returns an integer vector
+ containing these factors. Some factors may be repeated.
+
+ \code{divisors(n)} finds all the integers which divide
+ the integer \code{n}, and returns them as a sorted vector of integers
+ (beginning with \code{1} and ending with \code{n}).
+
+ \code{relatively.prime(n, m)} returns \code{TRUE} if the integers
+ \code{n} and \code{m} are relatively prime, that is, if they have no
+ common factors.
+
+ \code{least.common.multiple} and \code{greatest.common.divisor}
+ return the least common multiple or greatest common divisor of two
+ integers \code{n} and \code{m}.
+
+ \code{primesbelow(nmax)} returns an integer vector containing all the
+ prime numbers less than or equal to \code{nmax}.
+}
+\value{
+ \code{is.prime} and \code{relatively.prime} return a logical value.
+
+ \code{least.common.multiple} and \code{greatest.common.divisor}
+ return a single integer.
+
+ \code{primefactors} and \code{primesbelow}
+ return an integer vector.
+}
+\author{
+ \adrian.
+}
+\examples{
+ is.prime(17)
+
+ relatively.prime(2, 3)
+
+ primefactors(24) ## Note repeated factors
+
+ divisors(24)
+
+ greatest.common.divisor(60, 100)
+
+ primesbelow(20)
+}
+\keyword{math}
diff --git a/man/resolve.defaults.Rd b/man/resolve.defaults.Rd
new file mode 100644
index 0000000..2b996a7
--- /dev/null
+++ b/man/resolve.defaults.Rd
@@ -0,0 +1,87 @@
+\name{resolve.defaults}
+\alias{resolve.defaults}
+\alias{resolve.1.default}
+\title{
+ Determine Values of Variables Using Several Default Rules
+}
+\description{
+ Determine the values of variables by applying
+ several different default rules in a given order.
+}
+\usage{
+ resolve.defaults(\dots, .MatchNull = TRUE, .StripNull = FALSE)
+
+ resolve.1.default(.A, \dots)
+}
+\arguments{
+ \item{\dots}{
+ Several lists of \code{name=value} pairs.
+ }
+ \item{.MatchNull}{
+ Logical value. If \code{TRUE} (the default), an entry of the form
+ \code{name=NULL} will be treated as assigning the value
+ \code{NULL} to the variable \code{name}. If \code{FALSE},
+ such entries will be ignored.
+ }
+ \item{.StripNull}{
+ Logical value indicating whether
+ entries of the form \code{name=NULL} should be removed
+ from the result.
+ }
+ \item{.A}{
+ Either a character string giving the name of the variable
+ to be extracted, or a list consisting of one \code{name=value} pair
+ giving the variable name and its fallback default value.
+ }
+}
+\details{
+ These functions determine the values of variables by applying a series
+ of default rules, in the order specified.
+
+ Each of the arguments \code{\dots} should be a list of
+ \code{name=value} pairs giving a value
+ for a variable \code{name}. Each list could represent a
+ set of arguments given by the user, or a
+ rule assigning default values to some variables.
+ Lists that appear earlier in the sequence of arguments \code{\dots}
+ take precedence.
+
+ The arguments \code{\dots} will be concatenated into a single list.
+ The earliest occurrence of each \code{name} is then used to
+ determine the final value of the variable \code{name}.
+
+ The function \code{resolve.defaults} returns a list of
+ \code{name=value} pairs for all variables encountered.
+ It is commonly used to decide the values of
+ arguments to be passed to another function
+ using \code{\link[base]{do.call}}.
+
+ The function \code{resolve.1.default} returns the value
+ of the specified variable as determined by \code{resolve.defaults}.
+ It is commonly used inside a function to determine the value
+ of an argument.
+}
+\value{
+ The result of \code{resolve.defaults} is
+ a list of \code{name=value} pairs.
+
+ The result of \code{resolve.1.default} can be any kind of value.
+}
+\author{
+ \adrian
+}
+\seealso{
+ \code{\link[base]{do.call}}
+}
+\examples{
+ user <- list(day="Friday")
+ ruleA <- list(month="Jan", gravity=NULL)
+ ruleB <- list(day="Tuesday", month="May", gravity=42)
+ resolve.defaults(user, ruleA, ruleB)
+ resolve.defaults(user, ruleA, ruleB, .StripNull=TRUE)
+ resolve.defaults(user, ruleA, ruleB, .MatchNull=FALSE)
+
+ resolve.1.default("month", user, ruleA, ruleB)
+}
+\keyword{programming}
+\keyword{utilities}
diff --git a/man/revcumsum.Rd b/man/revcumsum.Rd
new file mode 100644
index 0000000..496fdfd
--- /dev/null
+++ b/man/revcumsum.Rd
@@ -0,0 +1,48 @@
+\name{revcumsum}
+\alias{revcumsum}
+\title{
+ Reverse Cumulative Sum
+}
+\description{
+ Returns a vector of cumulative sums of the input values,
+ running in reverse order. That is, the \code{i}th entry in the output
+ is the sum of entries \code{i} to \code{n} in the input,
+ where \code{n} is the length of the input.
+}
+\usage{
+revcumsum(x)
+}
+\arguments{
+ \item{x}{
+ A numeric or complex vector.
+ }
+}
+\details{
+ This low-level utility function is a faster alternative to
+ \code{\link[base]{rev}(\link[base]{cumsum}(\link[base]{rev}(x)))}
+ under certain conditions.
+ It computes the reverse cumulative sum of the entries of \code{x}.
+ If \code{y <- revcumsum(x)}, then \code{y[i] = sum(x[i:n])} where
+ \code{n = length(x)}.
+
+ This function should not be used if \code{x} could contain \code{NA}
+ values: this would lead to an error.
+}
+\value{
+ A vector of the same length and type as \code{x}.
+}
+\author{
+ \adrian.
+}
+\seealso{
+ \code{\link[base]{cumsum}}.
+}
+\examples{
+ revcumsum(1:5)
+ rev(cumsum(rev(1:5)))
+ x <- runif(1e6)
+ system.time(rev(cumsum(rev(x))))
+ system.time(revcumsum(x))
+}
+\keyword{arith}
+\keyword{utilities}
diff --git a/man/simplenumber.Rd b/man/simplenumber.Rd
new file mode 100644
index 0000000..6a0e060
--- /dev/null
+++ b/man/simplenumber.Rd
@@ -0,0 +1,48 @@
+\name{simplenumber}
+\alias{simplenumber}
+\title{
+ Simple Rational Number
+}
+\description{
+ Given a numeric value, try to express it
+ as a simple rational number.
+}
+\usage{
+simplenumber(x, unit = "", multiply = "*", tol = .Machine$double.eps)
+}
+\arguments{
+ \item{x}{
+ A single numeric value.
+ }
+ \item{unit}{
+ Optional. Character string giving the name of the unit
+ in which \code{x} is expressed. Typically an irrational number
+ such as \code{pi}. See Examples.
+ }
+ \item{multiply}{
+ Optional. Character string representing multiplication.
+ }
+ \item{tol}{
+ Numerical tolerance.
+ }
+}
+\details{
+ The code tries to express \code{x} as an integer \code{x=n}, or as
+ the reciprocal of an integer \code{x=1/n}, or
+ as a simple rational number \code{x = m/n}, where \code{m,n} are
+ small integers.
+}
+\value{
+ A character string representing the simple number,
+ or \code{NULL} if not successful.
+}
+\author{
+ \adrian.
+}
+\examples{
+ simplenumber(0.3)
+ simplenumber(0.333333333333333333333333)
+ x <- pi * 2/3
+ simplenumber(x/pi, "pi")
+}
+\keyword{symbolmath}
diff --git a/man/spatstat.utils-internal.Rd b/man/spatstat.utils-internal.Rd
new file mode 100755
index 0000000..71b4749
--- /dev/null
+++ b/man/spatstat.utils-internal.Rd
@@ -0,0 +1,212 @@
+\name{spatstat.utils-internal}
+\title{Internal Functions of spatstat.utils Package}
+\alias{adjustthinrange}
+\alias{apply23sum}
+\alias{as2vector}
+\alias{asNumericMatrix}
+\alias{assignDFcolumn}
+\alias{badprobability}
+\alias{blockdiagarray}
+\alias{blockdiagmatrix}
+\alias{check.finite}
+\alias{choptext}
+\alias{choptextline}
+\alias{complaining}
+\alias{distpl}
+\alias{distppl}
+\alias{distppll}
+\alias{distppllmin}
+\alias{dont.complain.about}
+\alias{dotexpr.to.call}
+\alias{dropifsingle}
+\alias{dround}
+\alias{ensure2vector}
+\alias{ensure3Darray}
+\alias{equispaced}
+\alias{eratosthenes}
+\alias{evenly.spaced}
+\alias{exceedsMaxArraySize}
+\alias{exhibitStringList}
+\alias{explain.ifnot}
+\alias{fakecallstring}
+\alias{fastFindInterval}
+\alias{fave.order}
+\alias{fillseq}
+\alias{findfirstfactor}
+\alias{firstfactor}
+\alias{flat.deparse}
+\alias{fontify}
+\alias{forbidNA}
+\alias{geomseq}
+\alias{getdataobjects}
+\alias{good.names}
+\alias{graphicsPars}
+\alias{gsubdot}
+\alias{indexCartesian}
+\alias{inject.expr}
+\alias{insertinlist}
+\alias{is.blank}
+\alias{is.parseable}
+\alias{lty2char}
+\alias{make.parseable}
+\alias{mapstrings}
+\alias{matchNameOrPosition}
+\alias{matcolall}
+\alias{matcolany}
+\alias{matcolsum}
+\alias{matrixsample}
+\alias{matrowall}
+\alias{matrowany}
+\alias{matrowsum}
+\alias{natozero}
+\alias{niceround}
+\alias{NNdist2segments}
+\alias{numalign}
+\alias{nzpaste}
+\alias{orderstats}
+\alias{orderwhich}
+\alias{\%orifnull\%} %DoNotExport
+%NAMESPACE export("%orifnull%")
+\alias{padtowidth}
+\alias{passthrough}
+\alias{paste.expr}
+\alias{pasteFormula}
+\alias{pasteN}
+\alias{prettydiscrete}
+\alias{prettyinside}
+\alias{prolongseq}
+\alias{ratiotweak}
+\alias{romansort}
+\alias{samefunction}
+\alias{sensiblevarname}
+\alias{short.deparse}
+\alias{singlestring}
+\alias{startinrange}
+\alias{strsplitretain}
+\alias{substringcount}
+\alias{there.is.no.try}
+\alias{trap.extra.arguments}
+\alias{truncline}
+\alias{uptrimat}
+\alias{validposint}
+\alias{variablesintext}
+\alias{warn.ignored.args}
+%%
+\alias{inside.triangle}
+\alias{overlap.trapezium}
+%%
+\alias{xysegXcircle}
+\alias{xysegMcircle}
+\alias{xysegPcircle}
+%%
+\alias{matchIntegerDataFrames}
+\description{
+ Internal utility functions of the \code{spatstat.utils} package.
+}
+\usage{
+adjustthinrange(ur,vstep,vr)
+apply23sum(x)
+as2vector(x)
+asNumericMatrix(x)
+assignDFcolumn(x, name, value, \dots)
+badprobability(x, NAvalue)
+blockdiagarray(\dots)
+blockdiagmatrix(\dots)
+check.finite(x, context, xname, fatal, usergiven)
+choptext(\dots, prefix, indent)
+choptextline(st, w, prefix, indent)
+complaining(whinge, fatal, value)
+distpl(p, l)
+distppl(p, l)
+distppll(p, l, mintype, method, listit)
+distppllmin(p, l, big)
+dont.complain.about(\dots)
+dotexpr.to.call(expr, dot, evaluator)
+dropifsingle(x)
+dround(x)
+ensure2vector(x)
+ensure3Darray(x)
+equispaced(z, reltol)
+eratosthenes(nmax, startset)
+evenly.spaced(x, tol)
+exhibitStringList(prefix, strings)
+exceedsMaxArraySize(\dots)
+explain.ifnot(expr, context)
+fakecallstring(fname, parlist)
+fastFindInterval(x, b, labels, reltol, dig.lab)
+fave.order(x)
+fillseq(x, step)
+findfirstfactor(x)
+firstfactor(x)
+flat.deparse(x)
+fontify(x, font)
+forbidNA(x, context, xname, fatal, usergiven)
+geomseq(from, to, length.out)
+getdataobjects(nama, envir, datalist, fatal)
+good.names(nama, defaults, suffices)
+graphicsPars(key)
+gsubdot(replacement, x)
+indexCartesian(nn)
+inject.expr(base, expr)
+insertinlist(x, i, y)
+is.blank(s)
+is.parseable(x)
+lty2char(i)
+make.parseable(x)
+mapstrings(x, map)
+matchNameOrPosition(expected, avail)
+matcolall(x)
+matcolany(x)
+matcolsum(x)
+matrixsample(mat, newdim, phase, scale, na.value)
+matrowall(x)
+matrowany(x)
+matrowsum(x)
+natozero(x)
+niceround(x, m)
+NNdist2segments(xp, yp, x0, y0, x1, y1, bigvalue)
+numalign(i, nmax, zero)
+nzpaste(\dots, sep, collapse)
+orderstats(x, k, decreasing)
+orderwhich(x, k, decreasing)
+a \%orifnull\% b
+padtowidth(a, b, justify)
+passthrough(.Fun, \dots, .Fname)
+paste.expr(x)
+pasteFormula(f)
+pasteN(\dots)
+prettydiscrete(x, n)
+prettyinside(x, \dots)
+prolongseq(x, newrange, step)
+ratiotweak(a, b, overzero, zerozero)
+romansort(x)
+samefunction(f, g)
+sensiblevarname(guess, fallback, maxlen)
+short.deparse(x, maxlen)
+singlestring(s, coll)
+startinrange(x0, dx, r)
+strsplitretain(x, split)
+substringcount(x,y)
+there.is.no.try(\dots)
+trap.extra.arguments(\dots, .Context, .Fatal)
+truncline(x, nc)
+uptrimat(x)
+validposint(n, caller, fatal)
+variablesintext(x)
+warn.ignored.args(\dots, context)
+%% xypolygon %%%
+inside.triangle(x, y, xx, yy)
+overlap.trapezium(xa, ya, xb, yb, verb = FALSE)
+%%
+xysegXcircle(xcentres, ycentres, radii, x0, y0, x1, y1, check)
+xysegMcircle(xcentres, ycentres, radmat, x0, y0, x1, y1, check)
+xysegPcircle(xc, yc, rc, x0, y0, x1, y1, check)
+%%
+matchIntegerDataFrames(X,Y,sort)
+}
+\details{
+ These internal \pkg{spatstat.utils} functions are not usually called
+ directly by the user. Their names and capabilities may change
+ without warning from one version of \pkg{spatstat.utils} to the next.
+}
+\keyword{internal}
diff --git a/man/spatstat.utils-package.Rd b/man/spatstat.utils-package.Rd
new file mode 100644
index 0000000..a84afad
--- /dev/null
+++ b/man/spatstat.utils-package.Rd
@@ -0,0 +1,121 @@
+\name{spatstat.utils-package}
+\alias{spatstat.utils-package}
+\alias{spatstat.utils}
+\docType{package}
+\title{The spatstat.utils Package}
+\description{
+ The \pkg{spatstat.utils} package contains low-level utilities,
+ written for the \pkg{spatstat} package, which may be useful in
+ other packages as well.
+}
+\details{
+ The functions in \pkg{spatstat.utils} were originally written
+ as internal, undocumented, utility functions in the \pkg{spatstat} package.
+
+ Many of these functions could be useful to other programmers,
+ so we have made them available in a separate package
+ \pkg{spatstat.utils} and provided documentation.
+
+ The functionality contained in \pkg{spatstat.utils} includes:
+ \describe{
+ \item{Factorisation of integers}{
+ Find prime numbers (\code{\link{primesbelow}}),
+ factorise a composite number into its prime factors
+ (\code{\link{primefactors}}),
+ determine whether a number is prime (\code{\link{is.prime}})
+ or whether two numbers are relatively prime
+ (\code{\link{relatively.prime}}),
+ and find the least common multiple
+ or greatest common divisor of two numbers
+ (\code{\link{least.common.multiple}},
+ \code{\link{greatest.common.divisor}}).
+ }
+ \item{Faster versions of basic \R tools}{
+ Faster versions of some basic \R tools and idioms are provided.
+ These are only faster in particular cases, but if you know that
+ your data have a particular form, the acceleration can be substantial.
+ See \code{\link{ifelseAB}}, \code{\link{fave.order}},
+ \code{\link{revcumsum}}, \code{\link{tapplysum}}.
+ }
+ \item{Grammar}{
+ Use the correct word in English to refer to
+ an ordinal number (\code{\link{ordinal}},
+ \code{\link{ordinalsuffix}})
+ and the correct indefinite article
+ (\code{\link{articlebeforenumber}}).
+ }
+ \item{Tools for generating printed output}{
+ The function \code{\link{splat}} is a replacement for
+ \code{cat(paste(\dots))} which ensures that output stays
+ inside the declared text margin (\code{\link{getOption}("width")})
+ and can also perform automatic indentation.
+ There are useful functions to add or remove parentheses
+ (\code{\link{paren}}, \code{\link{unparen}}) and
+ to make comma-separated lists (\code{\link{commasep}}).
+ }
+ \item{Handling intervals (ranges) of real numbers}{
+ Simple functions handle an interval (range) of numerical
+ values: \code{\link{check.range}},
+ \code{\link{intersect.ranges}},
+ \code{\link{inside.range}},
+ \code{\link{check.in.range}},
+ \code{\link{prange}}.
+ }
+ \item{Handling a formula}{
+ Tools for handling a formula in the \R language include
+ \code{\link{lhs.of.formula}},
+ \code{\link{rhs.of.formula}},
+ \code{\link{variablesinformula}},
+ \code{\link{termsinformula}},
+ \code{\link{offsetsinformula}},
+ \code{\link{can.be.formula}} and
+ \code{\link{identical.formulae}}.
+ }
+ \item{Polynomials}{
+ There are tools for creating and manipulating symbolic expressions
+ for polynomials, as they might appear in a formula
+ (\code{\link{sympoly}},
+ \code{\link{expand.polynom}}).
+ }
+ \item{Validating arguments}{
+ There are many tools for validating an argument and
+ generating a comprehensible error or warning message
+ if the argument is not valid:
+ \code{\link{check.1.integer}},
+ \code{\link{check.nvector}},
+ \code{\link{check.named.vector}}.
+ }
+ \item{Passing arguments}{
+ There are many tools for calling a function while passing
+ only some of the arguments in a supplied list of arguments:
+ \code{\link{do.call.matched}}, \code{\link{do.call.without}},
+ \code{\link{resolve.defaults}}.
+ }
+ \item{Traced optimization}{
+ \code{\link{optimizeWithTrace}} is a simple wrapper for
+ the one-dimensional optimization routine
+ \code{\link[stats]{optimize}}. It stores the values of the
+ function argument each time it is called, stores the resulting
+ function values, and returns them along with the optimal value.
+ }
+ \item{Workarounds}{
+ There are workarounds for known bugs or undesirable features in
+ other software.
+ \code{\link{spatstatLocator}} is a replacement for
+ \code{\link[graphics]{locator}} which works around a bug in the
+ \code{RStudio} graphics interface.
+ \code{\link{cat.factor}} concatenates several factors,
+ merging the levels, to produce a new factor.
+ }
+ }
+}
+\section{Licence}{
+ This library and its documentation are usable under the terms of the
+ \dQuote{GNU General Public License},
+ a copy of which is distributed with \R.
+}
+\author{
+ \spatstatAuthors.
+}
+\keyword{spatial}
+\keyword{package}
diff --git a/man/spatstatLocator.Rd b/man/spatstatLocator.Rd
new file mode 100644
index 0000000..715a124
--- /dev/null
+++ b/man/spatstatLocator.Rd
@@ -0,0 +1,49 @@
+\name{spatstatLocator}
+\alias{spatstatLocator}
+\title{
+ Graphical Input
+}
+\description{
+ This is an alternative to the \code{\link[graphics]{locator}}
+ function. It contains a workaround for a bug that occurs in
+ \code{RStudio}.
+}
+\usage{
+spatstatLocator(n, type = c("p", "l", "o", "n"), \dots)
+}
+\arguments{
+ \item{n}{
+ Optional. Maximum number of points to locate.
+ }
+ \item{type}{
+ Character specifying how to plot the locations.
+ If \code{"p"} or \code{"o"} the points are plotted;
+ if \code{"l"} or \code{"o"} they are joined by lines.
+ }
+ \item{\dots}{
+ Additional graphics parameters used to plot the locations.
+ }
+}
+\details{
+ This is a replacement/workaround for the \code{\link{locator}}
+ function in some versions of \pkg{RStudio} which do not seem to
+ recognise the option \code{type="p"}.
+
+ See \code{\link[graphics]{locator}} for a description of the
+ behaviour.
+}
+\value{
+ A list containing components \code{x} and \code{y} which are vectors
+ giving the coordinates of the identified points in the
+ user coordinate system, i.e., the one specified by \code{par("usr")}.
+}
+\author{
+ \spatstatAuthors.
+}
+\seealso{
+ \code{\link[graphics]{locator}}
+}
+\examples{
+ if(interactive()) locator(1, type="p")
+}
+\keyword{iplot}
diff --git a/man/splat.Rd b/man/splat.Rd
new file mode 100644
index 0000000..cdad0d1
--- /dev/null
+++ b/man/splat.Rd
@@ -0,0 +1,62 @@
+\name{splat}
+\alias{splat}
+\title{
+ Print Text Within Margins
+}
+\description{
+ Prints a given character string or strings
+ inside the text margin specified by \code{options("width")}.
+ Indents the text if required.
+}
+\usage{
+splat(\dots, indent = 0)
+}
+\arguments{
+ \item{\dots}{
+ Character strings, or other arguments acceptable to
+ \code{\link[base]{paste}}.
+ }
+ \item{indent}{
+ Optional. Indentation of the text.
+ Either an integer specifying the number of character positions
+ by which the text should be indented, or a character string whose
+ length determines the indentation.
+ }
+}
+\details{
+ \code{splat} stands for \sQuote{split cat}.
+
+ The command \code{splat(\dots)} is like \code{cat(paste(\dots))}
+ except that the output will be split into lines that can be
+ printed within the current text margin
+ specified by \code{\link[base]{getOption}("width")}.
+
+ The arguments \code{\dots} are first combined into a character vector
+ using \code{\link[base]{paste}}. Then they are split into words
+ separated by white space. A newline will be inserted whenever the next
+ word does not fit in the available text area.
+ (Words will not be broken, so the text margin could be exceeded
+ if any word is longer than \code{\link[base]{getOption}("width")}).
+
+ If any argument is a vector, each
+ element of the vector is treated as a separate line.
+ Existing newline characters in \code{\dots} are also respected.
+}
+\value{
+Null.
+}
+\author{
+ \spatstatAuthors.
+}
+\examples{
+ op <- options(width=20)
+ splat("There is more than one way to skin a cat.")
+ splat("There is more than one", "way to skin a cat.", indent=5)
+
+ options(width=10)
+ splat("The value of pi is", pi)
+ splat("The value of pi is", signif(pi))
+ options(op)
+}
+\keyword{print}
+\keyword{utilities}
diff --git a/man/tapplysum.Rd b/man/tapplysum.Rd
new file mode 100644
index 0000000..4640f39
--- /dev/null
+++ b/man/tapplysum.Rd
@@ -0,0 +1,60 @@
+\name{tapplysum}
+\alias{tapplysum}
+\title{
+ Sum By Factor Level
+}
+\description{
+ A faster equivalent of \code{tapply(FUN=sum)}.
+}
+\usage{
+ tapplysum(x, flist, do.names = FALSE, na.rm = TRUE)
+}
+\arguments{
+ \item{x}{
+ Numeric vector.
+ }
+ \item{flist}{
+ A list of \code{factor}s of the same length as \code{x}.
+ }
+ \item{do.names}{
+ Logical value indicating whether to attach names to the result.
+ }
+ \item{na.rm}{
+ Logical value indicating whether to remove \code{NA} values
+ before computing the sums.
+ }
+}
+\details{
+ This function is designed to be a faster alternative to
+ the idiom \code{y <- \link[base]{tapply}(x, flist, sum); y[is.na(y)] <- 0}.
+ The result \code{y} is a vector, matrix or array of
+ dimension equal to the number of factors in \code{flist}. Each
+ position in \code{y} represents one of the possible combinations
+ of the factor levels. The resulting value in this position
+ is the sum of all entries of \code{x}
+ where the factors in \code{flist} take this particular combination of
+ values. The sum is zero if this combination does not occur.
+
+ Currently this is implemented for the cases where \code{flist}
+ has length 2 or 3 (resulting in a matrix or 3D array, respectively).
+ For other cases we fall back on \code{\link[base]{tapply}}.
+}
+\value{
+ A numeric vector, matrix or array.
+}
+\author{
+ \adrian and Tilman Davies.
+}
+\seealso{
+ \code{\link[base]{tapply}}, \code{\link[base]{table}}
+}
+\examples{
+ x <- 1:12
+ a <- factor(rep(LETTERS[1:2], each=6))
+ b <- factor(rep(letters[1:4], times=3))
+ ff <- list(a, b)
+ tapply(x, ff, sum)
+ tapplysum(x, ff, do.names=TRUE)
+}
+\keyword{arith}
+\keyword{utilities}
diff --git a/man/termsinformula.Rd b/man/termsinformula.Rd
new file mode 100644
index 0000000..9eac67f
--- /dev/null
+++ b/man/termsinformula.Rd
@@ -0,0 +1,102 @@
+\name{termsinformula}
+\alias{termsinformula}
+\alias{offsetsinformula}
+\alias{variablesinformula}
+\alias{lhs.of.formula}
+\alias{lhs.of.formula<-}
+\alias{rhs.of.formula}
+\alias{rhs.of.formula<-}
+\alias{can.be.formula}
+\alias{identical.formulae}
+\title{
+ Manipulate Formulae
+}
+\description{
+ Operations for manipulating formulae.
+}
+\usage{
+termsinformula(x)
+variablesinformula(x)
+offsetsinformula(x)
+lhs.of.formula(x)
+rhs.of.formula(x, tilde=TRUE)
+lhs.of.formula(x) <- value
+rhs.of.formula(x) <- value
+can.be.formula(x)
+identical.formulae(x,y)
+}
+\arguments{
+ \item{x,y}{
+ Formulae, or character strings representing formulae.
+ }
+ \item{tilde}{
+ Logical value indicating whether to retain the tilde.
+ }
+ \item{value}{
+ Symbol or expression in the \R language. See Examples.
+ }
+}
+\details{
+ \code{variablesinformula(x)}
+ returns a character vector of the names
+ of all variables which appear in the formula \code{x}.
+
+ \code{termsinformula(x)} returns a character vector of all
+ terms in the formula \code{x} (after expansion of interaction terms).
+
+ \code{offsetsinformula(x)} returns a character vector of all
+ offset terms in the formula.
+
+ \code{rhs.of.formula(x)} returns the right-hand side of the formula
+ as another formula (that is, it removes the left-hand side) provided
+ \code{tilde=TRUE} (the default). If \code{tilde=FALSE}, then the
+ right-hand side is returned as a language object.
+
+ \code{lhs.of.formula(x)} returns the left-hand side of the formula
+ as a symbol or language object, or \code{NULL} if the formula has no
+ left-hand side.
+
+ \code{lhs.of.formula(x) <- value}
+ and \code{rhs.of.formula(x) <- value}
+ change the formula \code{x} by replacing the left or right hand side
+ of the formula by \code{value}.
+
+ \code{can.be.formula(x)} returns \code{TRUE} if \code{x} is a formula
+ or a character string that can be parsed as a formula, and returns
+ \code{FALSE} otherwise.
+
+ \code{identical.formulae(x,y)} returns \code{TRUE} if \code{x} and
+ \code{y} are identical formulae (ignoring their environments).
+}
+\value{
+ \code{variablesinformula},
+ \code{termsinformula} and
+ \code{offsetsinformula} return a character vector.
+
+ \code{rhs.of.formula} returns a formula.
+ \code{lhs.of.formula} returns a symbol or language object, or \code{NULL}.
+
+ \code{can.be.formula} and \code{identical.formulae} return
+ a logical value.
+}
+\author{
+ \spatstatAuthors.
+}
+\examples{
+ f <- (y ~ x + z*w + offset(h))
+ lhs.of.formula(f)
+ rhs.of.formula(f)
+ variablesinformula(f)
+ termsinformula(f)
+ offsetsinformula(f)
+ g <- f
+ environment(g) <- new.env()
+ identical(f,g)
+ identical.formulae(f,g)
+ lhs.of.formula(f) <- quote(mork) # or as.name("mork")
+ f
+ rhs.of.formula(f) <- quote(x+y+z) # or parse(text="x+y+z")[[1]]
+ f
+}
+\keyword{utilities}
+
diff --git a/man/verbalogic.Rd b/man/verbalogic.Rd
new file mode 100644
index 0000000..8e223c2
--- /dev/null
+++ b/man/verbalogic.Rd
@@ -0,0 +1,67 @@
+\name{verbalogic}
+\alias{verbalogic}
+\title{
+ Verbal Logic
+}
+\description{
+ Perform the specified logical operation
+ on the character vector \code{x}, recognising the
+ special strings \code{"TRUE"} and \code{"FALSE"}
+ and treating other strings as logical variables.
+}
+\usage{
+verbalogic(x, op = "and")
+}
+\arguments{
+ \item{x}{
+ Character vector.
+ }
+ \item{op}{
+ Logical operation: one of the character strings
+ \code{"and"}, \code{"or"} or \code{"not"}.
+ }
+}
+\details{
+ This function performs simple logical operations
+ on character strings that represent human-readable statements.
+
+ The character vector \code{x} may contain any strings:
+ the special strings \code{"TRUE"} and \code{"FALSE"} are
+ treated as the logical values \code{TRUE} and \code{FALSE},
+ while all other strings are treated as if they were
+ logical variables.
+
+ If \code{op="and"}, the result is a single string,
+ logically equivalent to \code{x[1] && x[2] && ... && x[n]}.
+ First, any entries of \code{x} equal to \code{"TRUE"} are removed.
+ The result is \code{"FALSE"} if any of the entries of \code{x}
+ is \code{"FALSE"}; otherwise it is equivalent to
+ \code{paste(x, collapse=" and ")}.
+
+ If \code{op="or"}, the result is a single string, logically equivalent to
+ \code{x[1] || x[2] || ... || x[n]}.
+ First, any entries of \code{x} equal to \code{"FALSE"} are removed.
+ The result is \code{"TRUE"} if any of the entries of \code{x}
+ is \code{"TRUE"}; otherwise it is equivalent to
+ \code{paste(x, collapse=" or ")}.
+
+ If \code{op="not"}, the result is a character vector \code{y}
+ such that \code{y[i]} is the logical negation of \code{x[i]}.
+
+ The code does not understand English grammar and cannot expand logical
+ expressions.
+}
+\value{
+ A character string.
+}
+\author{
+ \adrian.
+}
+\examples{
+ x <- c("The sky is blue", "my name is not Einstein", "FALSE")
+ verbalogic(x, "and")
+ verbalogic(x, "or")
+ verbalogic(x, "not")
+}
+\keyword{logic}
+\keyword{manip}
diff --git a/src/chunkloop.h b/src/chunkloop.h
new file mode 100644
index 0000000..f0adceb
--- /dev/null
+++ b/src/chunkloop.h
@@ -0,0 +1,37 @@
+/*
+ chunkloop.h
+
+ Divide a loop into chunks
+
+ Convenient for divide-and-recombine,
+ and reducing calls to R_CheckUserInterrupt, etc.
+
+ $Revision: 1.2 $ $Date: 2013/05/27 02:09:10 $
+
+*/
+
+#define OUTERCHUNKLOOP(IVAR, LOOPLENGTH, ICHUNK, CHUNKSIZE) \
+ IVAR = 0; \
+ ICHUNK = 0; \
+ while(IVAR < LOOPLENGTH)
+
+#define INNERCHUNKLOOP(IVAR, LOOPLENGTH, ICHUNK, CHUNKSIZE) \
+ ICHUNK += CHUNKSIZE; \
+ if(ICHUNK > LOOPLENGTH) ICHUNK = LOOPLENGTH; \
+ for(; IVAR < ICHUNK; IVAR++)
+
+#define XOUTERCHUNKLOOP(IVAR, ISTART, IEND, ICHUNK, CHUNKSIZE) \
+ IVAR = ISTART; \
+ ICHUNK = 0; \
+ while(IVAR <= IEND)
+
+#define XINNERCHUNKLOOP(IVAR, ISTART, IEND, ICHUNK, CHUNKSIZE) \
+ ICHUNK += CHUNKSIZE; \
+ if(ICHUNK > IEND) ICHUNK = IEND; \
+ for(; IVAR <= IEND; IVAR++)
+
+#define CHUNKLOOP_H
+
+
+
+
diff --git a/src/circxseg.c b/src/circxseg.c
new file mode 100644
index 0000000..43f215e
--- /dev/null
+++ b/src/circxseg.c
@@ -0,0 +1,930 @@
+/*
+
+ circxseg.c
+
+ Intersections between circles and line segments
+
+ circXseg: centres * radii * segments
+
+ circMseg: matrix of radii with rows corresponding to centres
+
+ circPseg: parallel vectors of centres and radii
+
+ $Revision: 1.8 $ $Date: 2017/01/21 10:50:32 $
+
+ */
+
+#include <R.h>
+#include <Rdefines.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include <math.h>
+
+#undef BUGGER
+
+/* circXseg: consider every combination of centre, radius, segment */
+
+SEXP circXseg(SEXP XC, /* circle centres */
+ SEXP YC,
+ SEXP R, /* radii */
+ SEXP X0, /* segments */
+ SEXP Y0,
+ SEXP X1,
+ SEXP Y1
+ )
+{
+ double *xc, *yc, *r, *x0, *y0, *x1, *y1;
+ int Nc, Ns, Nr, Ne, NeMax, newmax;
+
+ /* outputs */
+ int *ie, *je, *ke; /* provenance of each intersection */
+ double *xe, *ye, *sinalpha; /* cut points and angles */
+ SEXP out, iout, jout, kout, xout, yout, sinout;
+ int *ip, *jp, *kp;
+ double *xp, *yp, *sp;
+
+ /* internal */
+ int i, j, k, m;
+ double xci, yci, rk, x0c, y0c, dx, dy, A, B, C, Det, sqrtDet, sina;
+ double u, u1, u2, slope, intcept, xcut, ycut, xnorm, ynorm, hx, hy;
+ double twoA, fourA, Bsquared, Cdist2;
+
+ PROTECT(XC = AS_NUMERIC(XC));
+ PROTECT(YC = AS_NUMERIC(YC));
+ PROTECT(R = AS_NUMERIC(R));
+ PROTECT(X0 = AS_NUMERIC(X0));
+ PROTECT(Y0 = AS_NUMERIC(Y0));
+ PROTECT(X1 = AS_NUMERIC(X1));
+ PROTECT(Y1 = AS_NUMERIC(Y1));
+ /* That's 7 protected */
+
+ /* get pointers */
+ xc = NUMERIC_POINTER(XC);
+ yc = NUMERIC_POINTER(YC);
+ r = NUMERIC_POINTER(R);
+ x0 = NUMERIC_POINTER(X0);
+ y0 = NUMERIC_POINTER(Y0);
+ x1 = NUMERIC_POINTER(X1);
+ y1 = NUMERIC_POINTER(Y1);
+
+ /* determine lengths of vectors */
+ Nc = LENGTH(XC);
+ Nr = LENGTH(R);
+ Ns = LENGTH(X0);
+
+ /* Guess amount of storage required */
+ NeMax = 4 * (Ns + Nr + Nc);
+ if(NeMax < 2048) NeMax = 2048;
+ ie = (int *) R_alloc(NeMax, sizeof(int));
+ je = (int *) R_alloc(NeMax, sizeof(int));
+ ke = (int *) R_alloc(NeMax, sizeof(int));
+ xe = (double *) R_alloc(NeMax, sizeof(double));
+ ye = (double *) R_alloc(NeMax, sizeof(double));
+ sinalpha = (double *) R_alloc(NeMax, sizeof(double));
+
+ /* initialise output */
+ Ne = 0;
+
+ if(Nc > 0 && Ns > 0 && Nr > 0) {
+ /* loop over circle centres */
+ for(i = 0; i < Nc; i++) {
+#ifdef BUGGER
+ Rprintf("Circle %d\n", i);
+#endif
+ R_CheckUserInterrupt();
+ xci = xc[i];
+ yci = yc[i];
+ /* loop over segments */
+ for(j = 0; j < Ns; j++) {
+#ifdef BUGGER
+ Rprintf("\tSegment %d\n", j);
+#endif
+ dx = x1[j] - x0[j];
+ dy = y1[j] - y0[j];
+ x0c = x0[j] - xci;
+ y0c = y0[j] - yci;
+ /* find intersections between circle and infinite line */
+ A = dx * dx + dy * dy;
+ B = 2 * (dx * x0c + dy * y0c);
+ twoA = 2.0 * A;
+ fourA = 4.0 * A;
+ Bsquared = B * B;
+ Cdist2 = x0c * x0c + y0c * y0c;
+ /* loop over radii */
+ for(k = 0; k < Nr; k++) {
+#ifdef BUGGER
+ Rprintf("\t\tRadius %d\n", k);
+#endif
+ rk = r[k];
+ C = Cdist2 - rk * rk;
+ Det = Bsquared - fourA * C;
+ if(Det > 0.0) {
+ /* two intersection points */
+ sqrtDet = sqrt(Det);
+ u1 = (-B - sqrtDet)/twoA;
+ u2 = (-B + sqrtDet)/twoA;
+ if(u1 > 0.0 && u1 < 1.0) {
+ /* first intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u1 * dx;
+ ycut = y0c + u1 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u1 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rk;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ if(u2 > 0.0 && u2 < 1.0) {
+ /* second intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u2 * dx;
+ ycut = y0c + u2 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u2 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rk;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ } else if(Det == 0.0) {
+ /* tangent point */
+ u = -B/twoA;
+ if(u > 0.0 && u < 1.0) {
+ /* tangent point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u * dx;
+ ycut = y0c + u * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rk;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* pack up */
+ PROTECT(out = NEW_LIST(6));
+ PROTECT(iout = NEW_INTEGER(Ne));
+ PROTECT(jout = NEW_INTEGER(Ne));
+ PROTECT(kout = NEW_INTEGER(Ne));
+ PROTECT(xout = NEW_NUMERIC(Ne));
+ PROTECT(yout = NEW_NUMERIC(Ne));
+ PROTECT(sinout = NEW_NUMERIC(Ne));
+ /* 7 + 1 + 6 = 14 protected */
+ ip = INTEGER_POINTER(iout);
+ jp = INTEGER_POINTER(jout);
+ kp = INTEGER_POINTER(kout);
+ xp = NUMERIC_POINTER(xout);
+ yp = NUMERIC_POINTER(yout);
+ sp = NUMERIC_POINTER(sinout);
+ for(m = 0; m < Ne; m++) {
+ ip[m] = ie[m];
+ jp[m] = je[m];
+ kp[m] = ke[m];
+ xp[m] = xe[m];
+ yp[m] = ye[m];
+ sp[m] = sinalpha[m];
+ }
+ SET_VECTOR_ELT(out, 0, xout);
+ SET_VECTOR_ELT(out, 1, yout);
+ SET_VECTOR_ELT(out, 2, iout);
+ SET_VECTOR_ELT(out, 3, jout);
+ SET_VECTOR_ELT(out, 4, kout);
+ SET_VECTOR_ELT(out, 5, sinout);
+ UNPROTECT(14);
+ return(out);
+}
+
+/* circMseg: matrix of radii with rows corresponding to centres */
+
+SEXP circMseg(SEXP XC, /* circle centres */
+ SEXP YC,
+ SEXP R, /* radii */
+ SEXP X0, /* segments */
+ SEXP Y0,
+ SEXP X1,
+ SEXP Y1
+ )
+{
+ double *xc, *yc, *r, *x0, *y0, *x1, *y1;
+ int Nc, Ns, Nr, Ne, NeMax, newmax;
+
+ /* outputs */
+ int *ie, *je, *ke; /* provenance of each intersection */
+ double *xe, *ye, *sinalpha; /* cut points and angles */
+ SEXP out, iout, jout, kout, xout, yout, sinout;
+ int *ip, *jp, *kp;
+ double *xp, *yp, *sp;
+
+ /* internal */
+ int i, j, k, m;
+ double xci, yci, rik, x0c, y0c, dx, dy, A, B, C, Det, sqrtDet, sina;
+ double u, u1, u2, slope, intcept, xcut, ycut, xnorm, ynorm, hx, hy;
+ double twoA, fourA, Bsquared, Cdist2;
+
+ PROTECT(XC = AS_NUMERIC(XC));
+ PROTECT(YC = AS_NUMERIC(YC));
+ PROTECT(R = AS_NUMERIC(R));
+ PROTECT(X0 = AS_NUMERIC(X0));
+ PROTECT(Y0 = AS_NUMERIC(Y0));
+ PROTECT(X1 = AS_NUMERIC(X1));
+ PROTECT(Y1 = AS_NUMERIC(Y1));
+ /* That's 7 protected */
+
+ /* get pointers */
+ xc = NUMERIC_POINTER(XC);
+ yc = NUMERIC_POINTER(YC);
+ r = NUMERIC_POINTER(R);
+ x0 = NUMERIC_POINTER(X0);
+ y0 = NUMERIC_POINTER(Y0);
+ x1 = NUMERIC_POINTER(X1);
+ y1 = NUMERIC_POINTER(Y1);
+
+ /* determine lengths of vectors */
+ Nc = LENGTH(XC);
+ Nr = LENGTH(R)/Nc; /* n.b. */
+ Ns = LENGTH(X0);
+
+ /* Guess amount of storage required */
+ NeMax = 4 * Nr * Nc;
+ if(NeMax < 2048) NeMax = 2048;
+ ie = (int *) R_alloc(NeMax, sizeof(int));
+ je = (int *) R_alloc(NeMax, sizeof(int));
+ ke = (int *) R_alloc(NeMax, sizeof(int));
+ xe = (double *) R_alloc(NeMax, sizeof(double));
+ ye = (double *) R_alloc(NeMax, sizeof(double));
+ sinalpha = (double *) R_alloc(NeMax, sizeof(double));
+
+ /* initialise output */
+ Ne = 0;
+
+ if(Nc > 0 && Ns > 0 && Nr > 0) {
+ /* loop over circle centres */
+ for(i = 0; i < Nc; i++) {
+#ifdef BUGGER
+ Rprintf("Circle %d\n", i);
+#endif
+ R_CheckUserInterrupt();
+ xci = xc[i];
+ yci = yc[i];
+ /* loop over segments */
+ for(j = 0; j < Ns; j++) {
+#ifdef BUGGER
+ Rprintf("\tSegment %d\n", j);
+#endif
+ dx = x1[j] - x0[j];
+ dy = y1[j] - y0[j];
+ x0c = x0[j] - xci;
+ y0c = y0[j] - yci;
+ /* find intersections between circle and infinite line */
+ A = dx * dx + dy * dy;
+ B = 2 * (dx * x0c + dy * y0c);
+ twoA = 2.0 * A;
+ fourA = 4.0 * A;
+ Bsquared = B * B;
+ Cdist2 = x0c * x0c + y0c * y0c;
+ /* loop over radii */
+ for(k = 0; k < Nr; k++) {
+#ifdef BUGGER
+ Rprintf("\t\tRadius [%d, %d]\n", i, k);
+#endif
+ rik = r[i + k*Nc];
+ C = Cdist2 - rik * rik;
+ Det = Bsquared - fourA * C;
+ if(Det > 0.0) {
+ /* two intersection points */
+ sqrtDet = sqrt(Det);
+ u1 = (-B - sqrtDet)/twoA;
+ u2 = (-B + sqrtDet)/twoA;
+ if(u1 > 0.0 && u1 < 1.0) {
+ /* first intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u1 * dx;
+ ycut = y0c + u1 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u1 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rik;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ if(u2 > 0.0 && u2 < 1.0) {
+ /* second intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u2 * dx;
+ ycut = y0c + u2 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u2 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rik;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ } else if(Det == 0.0) {
+ /* tangent point */
+ u = -B/twoA;
+ if(u > 0.0 && u < 1.0) {
+ /* tangent point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u * dx;
+ ycut = y0c + u * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rik;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ke[Ne] = k + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ ke = (int *) S_realloc((char *) ke,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* pack up */
+ PROTECT(out = NEW_LIST(6));
+ PROTECT(iout = NEW_INTEGER(Ne));
+ PROTECT(jout = NEW_INTEGER(Ne));
+ PROTECT(kout = NEW_INTEGER(Ne));
+ PROTECT(xout = NEW_NUMERIC(Ne));
+ PROTECT(yout = NEW_NUMERIC(Ne));
+ PROTECT(sinout = NEW_NUMERIC(Ne));
+ /* 7 + 1 + 6 = 14 protected */
+ ip = INTEGER_POINTER(iout);
+ jp = INTEGER_POINTER(jout);
+ kp = INTEGER_POINTER(kout);
+ xp = NUMERIC_POINTER(xout);
+ yp = NUMERIC_POINTER(yout);
+ sp = NUMERIC_POINTER(sinout);
+ for(m = 0; m < Ne; m++) {
+ ip[m] = ie[m];
+ jp[m] = je[m];
+ kp[m] = ke[m];
+ xp[m] = xe[m];
+ yp[m] = ye[m];
+ sp[m] = sinalpha[m];
+ }
+ SET_VECTOR_ELT(out, 0, xout);
+ SET_VECTOR_ELT(out, 1, yout);
+ SET_VECTOR_ELT(out, 2, iout);
+ SET_VECTOR_ELT(out, 3, jout);
+ SET_VECTOR_ELT(out, 4, kout);
+ SET_VECTOR_ELT(out, 5, sinout);
+ UNPROTECT(14);
+ return(out);
+}
+
+/* circPseg: centres and radii matched ('vectors') */
+
+SEXP circPseg(SEXP XC, /* circles (x0, y0, r) */
+ SEXP YC,
+ SEXP RC,
+ SEXP X0, /* segments */
+ SEXP Y0,
+ SEXP X1,
+ SEXP Y1
+ )
+{
+ double *xc, *yc, *rc, *x0, *y0, *x1, *y1;
+ int Nc, Ns, Ne, NeMax, newmax;
+
+ /* outputs */
+ int *ie, *je; /* provenance of each intersection */
+ double *xe, *ye, *sinalpha; /* cut points and angles */
+ SEXP out, iout, jout, xout, yout, sinout;
+ int *ip, *jp;
+ double *xp, *yp, *sp;
+
+ /* internal */
+ int i, j, m;
+ double xci, yci, rci, x0c, y0c, dx, dy, A, B, C, Det, sqrtDet, sina;
+ double u, u1, u2, slope, intcept, xcut, ycut, xnorm, ynorm, hx, hy;
+ double twoA, rci2;
+
+ PROTECT(XC = AS_NUMERIC(XC));
+ PROTECT(YC = AS_NUMERIC(YC));
+ PROTECT(RC = AS_NUMERIC(RC));
+ PROTECT(X0 = AS_NUMERIC(X0));
+ PROTECT(Y0 = AS_NUMERIC(Y0));
+ PROTECT(X1 = AS_NUMERIC(X1));
+ PROTECT(Y1 = AS_NUMERIC(Y1));
+ /* That's 7 protected */
+
+ /* get pointers */
+ xc = NUMERIC_POINTER(XC);
+ yc = NUMERIC_POINTER(YC);
+ rc = NUMERIC_POINTER(RC);
+ x0 = NUMERIC_POINTER(X0);
+ y0 = NUMERIC_POINTER(Y0);
+ x1 = NUMERIC_POINTER(X1);
+ y1 = NUMERIC_POINTER(Y1);
+
+ /* determine lengths of vectors */
+ Nc = LENGTH(XC);
+ Ns = LENGTH(X0);
+
+ /* Guess amount of storage required */
+ NeMax = 4 * (Ns + Nc);
+ if(NeMax < 2048) NeMax = 2048;
+ ie = (int *) R_alloc(NeMax, sizeof(int));
+ je = (int *) R_alloc(NeMax, sizeof(int));
+ xe = (double *) R_alloc(NeMax, sizeof(double));
+ ye = (double *) R_alloc(NeMax, sizeof(double));
+ sinalpha = (double *) R_alloc(NeMax, sizeof(double));
+
+ /* initialise output */
+ Ne = 0;
+
+ if(Nc > 0 && Ns > 0) {
+ /* loop over circles */
+ for(i = 0; i < Nc; i++) {
+
+#ifdef BUGGER
+ Rprintf("Circle %d\n", i);
+#endif
+ R_CheckUserInterrupt();
+
+ xci = xc[i];
+ yci = yc[i];
+ rci = rc[i];
+
+ rci2 = rci * rci;
+
+ /* loop over segments */
+ for(j = 0; j < Ns; j++) {
+
+#ifdef BUGGER
+ Rprintf("\tSegment %d\n", j);
+#endif
+
+ dx = x1[j] - x0[j];
+ dy = y1[j] - y0[j];
+ x0c = x0[j] - xci;
+ y0c = y0[j] - yci;
+ /* find intersections between circle and infinite line */
+ A = dx * dx + dy * dy;
+ B = 2 * (dx * x0c + dy * y0c);
+ C = x0c * x0c + y0c * y0c - rci2;
+ Det = B * B - 4.0 * A * C;
+ twoA = 2.0 * A;
+ if(Det > 0.0) {
+ /* two intersection points */
+ sqrtDet = sqrt(Det);
+ u1 = (-B - sqrtDet)/twoA;
+ u2 = (-B + sqrtDet)/twoA;
+ if(u1 > 0.0 && u1 < 1.0) {
+ /* first intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u1 * dx;
+ ycut = y0c + u1 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u1 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rci;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ if(u2 > 0.0 && u2 < 1.0) {
+ /* second intersection point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u2 * dx;
+ ycut = y0c + u2 * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u2 * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rci;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ } else if(Det == 0.0) {
+ /* tangent point */
+ u = -B/twoA;
+ if(u > 0.0 && u < 1.0) {
+ /* tangent point is inside segment */
+ if(dx != 0) {
+ /* sloping line */
+ slope = dy/dx;
+ intcept = y0c - slope * x0c;
+ xcut = x0c + u * dx;
+ ycut = y0c + u * dy;
+ ynorm = intcept/(slope * slope + 1.0);
+ xnorm = - ynorm * slope;
+ } else {
+ /* vertical line */
+ xcut = x0c;
+ ycut = y0c + u * dy;
+ xnorm = x0c;
+ ynorm = 0.0;
+ }
+ hx = xcut - xnorm;
+ hy = ycut - ynorm;
+
+ sina = sqrt(hx * hx + hy * hy)/rci;
+ if(sina > 1.0) sina = 1.0; else if(sina < -1.0) sina = -1.0;
+
+ /* add to output */
+#ifdef BUGGER
+ Rprintf("\t\t\tAdding..\n");
+#endif
+ sinalpha[Ne] = sina;
+ xe[Ne] = xcut + xci;
+ ye[Ne] = ycut + yci;
+ ie[Ne] = i + 1;
+ je[Ne] = j + 1;
+ ++Ne;
+ if(Ne >= NeMax) {
+ /* storage overflow; reallocate */
+#ifdef BUGGER
+ Rprintf("\t\t\tOVERFLOW..\n");
+#endif
+ newmax = 2 * NeMax;
+ xe = (double *) S_realloc((char *) xe,
+ newmax, NeMax, sizeof(double));
+ ye = (double *) S_realloc((char *) ye,
+ newmax, NeMax, sizeof(double));
+ ie = (int *) S_realloc((char *) ie,
+ newmax, NeMax, sizeof(int));
+ je = (int *) S_realloc((char *) je,
+ newmax, NeMax, sizeof(int));
+ sinalpha = (double *) S_realloc((char *) sinalpha,
+ newmax, NeMax, sizeof(double));
+ NeMax = newmax;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /* pack up */
+ PROTECT(out = NEW_LIST(5));
+ PROTECT(iout = NEW_INTEGER(Ne));
+ PROTECT(jout = NEW_INTEGER(Ne));
+ PROTECT(xout = NEW_NUMERIC(Ne));
+ PROTECT(yout = NEW_NUMERIC(Ne));
+ PROTECT(sinout = NEW_NUMERIC(Ne));
+ /* 7 + 1 + 5 = 13 protected */
+ ip = INTEGER_POINTER(iout);
+ jp = INTEGER_POINTER(jout);
+ xp = NUMERIC_POINTER(xout);
+ yp = NUMERIC_POINTER(yout);
+ sp = NUMERIC_POINTER(sinout);
+ for(m = 0; m < Ne; m++) {
+ ip[m] = ie[m];
+ jp[m] = je[m];
+ xp[m] = xe[m];
+ yp[m] = ye[m];
+ sp[m] = sinalpha[m];
+ }
+ SET_VECTOR_ELT(out, 0, xout);
+ SET_VECTOR_ELT(out, 1, yout);
+ SET_VECTOR_ELT(out, 2, iout);
+ SET_VECTOR_ELT(out, 3, jout);
+ SET_VECTOR_ELT(out, 4, sinout);
+ UNPROTECT(13);
+ return(out);
+}
diff --git a/src/distseg.c b/src/distseg.c
new file mode 100755
index 0000000..21b7847
--- /dev/null
+++ b/src/distseg.c
@@ -0,0 +1,175 @@
+/*
+ distseg.c
+
+ Distances from point pattern to line segment pattern
+ Distance transform of a line segment pattern
+
+ nndist2segs: minimum distance from point to any line segment
+ prdist2segs: pairwise distances from each point to each line segment
+
+ $Revision: 1.9 $ $Date: 2012/03/27 05:38:51 $
+
+ Author: Adrian Baddeley
+
+*/
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include <math.h>
+
+#include "chunkloop.h"
+
+void
+nndist2segs(xp, yp, npoints, x0, y0, x1, y1, nsegments, epsilon, dist2, index)
+ /* input */
+ double *xp, *yp; /* point/pixel coordinates */
+ int *npoints;
+ double *x0, *y0, *x1, *y1; /* line segment endpoints */
+ int *nsegments;
+ double *epsilon; /* tolerance for short segments */
+ /* output */
+ double *dist2; /* squared distance from pixel
+ to nearest line segment
+ INITIALISED TO LARGE VALUE
+ */
+ int *index; /* which line segment is closest */
+{
+ int i,j, np, nseg, maxchunk;
+ double dx,dy,leng,co,si; /* parameters of segment */
+ double xdif0,ydif0,xdif1,ydif1,xpr,ypr; /* vectors */
+ double dsq0,dsq1,dsq,dsqperp; /* squared distances */
+ double eps;
+
+ np = *npoints;
+ nseg = *nsegments;
+ eps = *epsilon;
+
+ OUTERCHUNKLOOP(j, nseg, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(j, nseg, maxchunk, 16384) {
+ dx = x1[j] - x0[j];
+ dy = y1[j] - y0[j];
+ leng = hypot(dx, dy);
+ if(leng > eps) {
+ /* normal case */
+ co = dx/leng;
+ si = dy/leng;
+ for(i = 0; i < np; i++) {
+ /* vectors from pixel to segment endpoints */
+ xdif0 = xp[i] - x0[j];
+ ydif0 = yp[i] - y0[j];
+ xdif1 = xp[i] - x1[j];
+ ydif1 = yp[i] - y1[j];
+ /* squared distances to segment endpoints */
+ dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
+ dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
+ dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
+ /* rotate pixel around 1st endpoint of segment
+ so that line segment lies in x axis */
+ xpr = xdif0 * co + ydif0 * si;
+ ypr = -xdif0 * si + ydif0 * co;
+ /* perpendicular distance applies only in perpendicular region */
+ if(xpr >= 0.0 && xpr <= leng) {
+ dsqperp = ypr * ypr;
+ if(dsqperp < dsq) dsq = dsqperp;
+ }
+ if(dist2[i] > dsq) {
+ dist2[i] = dsq;
+ index[i] = j;
+ }
+ }
+ } else {
+ /* short segment - use endpoints only */
+ for(i = 0; i < np; i++) {
+ /* vectors from pixel to segment endpoints */
+ xdif0 = xp[i] - x0[j];
+ ydif0 = yp[i] - y0[j];
+ xdif1 = xp[i] - x1[j];
+ ydif1 = yp[i] - y1[j];
+ /* squared distances to segment endpoints */
+ dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
+ dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
+ dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
+ if(dist2[i] > dsq) {
+ dist2[i] = dsq;
+ index[i] = j;
+ }
+ }
+ }
+ }
+ }
+}
+
+void
+prdist2segs(xp, yp, npoints, x0, y0, x1, y1, nsegments, epsilon, dist2)
+ /* input */
+ double *xp, *yp; /* point/pixel coordinates */
+ int *npoints;
+ double *x0, *y0, *x1, *y1; /* line segment endpoints */
+ int *nsegments;
+ double *epsilon; /* tolerance for short segments */
+ /* output */
+ double *dist2; /* squared distances from each pixel
+ to each line segment */
+{
+ int i,j, np, nseg, maxchunk;
+ double dx,dy,leng,co,si; /* parameters of segment */
+ double xdif0,ydif0,xdif1,ydif1,xpr,ypr; /* vectors */
+ double dsq0,dsq1,dsq,dsqperp; /* squared distances */
+ double eps;
+
+ np = *npoints;
+ nseg = *nsegments;
+ eps = *epsilon;
+
+ OUTERCHUNKLOOP(j, nseg, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(j, nseg, maxchunk, 16384) {
+ dx = x1[j] - x0[j];
+ dy = y1[j] - y0[j];
+ leng = hypot(dx, dy);
+ if(leng > eps) {
+ /* normal case */
+ co = dx/leng;
+ si = dy/leng;
+ for(i = 0; i < np; i++) {
+ /* vectors from pixel to segment endpoints */
+ xdif0 = xp[i] - x0[j];
+ ydif0 = yp[i] - y0[j];
+ xdif1 = xp[i] - x1[j];
+ ydif1 = yp[i] - y1[j];
+ /* squared distances to segment endpoints */
+ dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
+ dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
+ dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
+ /* rotate pixel around 1st endpoint of segment
+ so that line segment lies in x axis */
+ xpr = xdif0 * co + ydif0 * si;
+ ypr = -xdif0 * si + ydif0 * co;
+ /* perpendicular distance applies only in perpendicular region */
+ if(xpr >= 0.0 && xpr <= leng) {
+ dsqperp = ypr * ypr;
+ if(dsqperp < dsq) dsq = dsqperp;
+ }
+ dist2[i + j * np] = dsq;
+ }
+ } else {
+ /* short segment */
+ for(i = 0; i < np; i++) {
+ /* vectors from pixel to segment endpoints */
+ xdif0 = xp[i] - x0[j];
+ ydif0 = yp[i] - y0[j];
+ xdif1 = xp[i] - x1[j];
+ ydif1 = yp[i] - y1[j];
+ /* squared distances to segment endpoints */
+ dsq0 = xdif0 * xdif0 + ydif0 * ydif0;
+ dsq1 = xdif1 * xdif1 + ydif1 * ydif1;
+ dsq = (dsq0 < dsq1) ? dsq0 : dsq1;
+ dist2[i + j * np] = dsq;
+ }
+ }
+ }
+ }
+}
+
diff --git a/src/fastinterv.c b/src/fastinterv.c
new file mode 100644
index 0000000..bda079f
--- /dev/null
+++ b/src/fastinterv.c
@@ -0,0 +1,37 @@
+#include <R.h>
+#include <Rdefines.h>
+
+/*
+
+ fastinterv.c
+
+ Fast version of findInterval when breaks are known to be evenly spaced
+ and are known to embrace the data.
+
+ $Revision: 1.2 $ $Date: 2015/10/19 11:09:18 $
+
+*/
+
+void fastinterv(x, n, brange, nintervals, y)
+ double *x, *brange;
+ int *n, *nintervals;
+ int *y;
+{
+ double bmin, bmax, db;
+ int i, j, m, N;
+
+ m = *nintervals;
+ N = *n;
+
+ bmin = brange[0];
+ bmax = brange[1];
+ db = (bmax - bmin)/m;
+
+ for(i = 0; i < N; i++) {
+ j = (int) ceil((x[i] - bmin)/db);
+ if(j <= 0) { j = 1; } else if(j > m) { j = m; }
+ y[i] = j;
+ }
+}
+
+
diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..eeb6ea8
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,47 @@
+
+/*
+ Native symbol registration table for spatstat.utils package
+
+ Automatically generated - do not edit this file!
+
+*/
+
+#include "proto.h"
+#include <R.h>
+#include <Rinternals.h>
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+/*
+ See proto.h for declarations for the native routines registered below.
+*/
+
+static const R_CMethodDef CEntries[] = {
+ {"Cmatchxy", (DL_FUNC) &Cmatchxy, 7},
+ {"drevcumsum", (DL_FUNC) &drevcumsum, 2},
+ {"fastinterv", (DL_FUNC) &fastinterv, 5},
+ {"inxyp", (DL_FUNC) &inxyp, 8},
+ {"irevcumsum", (DL_FUNC) &irevcumsum, 2},
+ {"nndist2segs", (DL_FUNC) &nndist2segs, 11},
+ {"ply2sum", (DL_FUNC) &ply2sum, 8},
+ {"ply3sum", (DL_FUNC) &ply3sum, 10},
+ {"prdist2segs", (DL_FUNC) &prdist2segs, 10},
+ {"primefax", (DL_FUNC) &primefax, 3},
+ {"CUmatch2int", (DL_FUNC) &CUmatch2int, 7},
+ {"CSmatch2int", (DL_FUNC) &CSmatch2int, 7},
+ {"CUmatch3int", (DL_FUNC) &CUmatch3int, 9},
+ {"CSmatch3int", (DL_FUNC) &CSmatch3int, 9},
+ {NULL, NULL, 0}
+};
+
+static const R_CallMethodDef CallEntries[] = {
+ {"circMseg", (DL_FUNC) &circMseg, 7},
+ {"circXseg", (DL_FUNC) &circXseg, 7},
+ {NULL, NULL, 0}
+};
+
+void R_init_spatstat_utils(DllInfo *dll)
+{
+ R_registerRoutines(dll, CEntries, CallEntries, NULL, NULL);
+ R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/src/inxyp.c b/src/inxyp.c
new file mode 100755
index 0000000..7ab3ea5
--- /dev/null
+++ b/src/inxyp.c
@@ -0,0 +1,75 @@
+/*
+ inxyp.c
+
+ Point-in-polygon test
+
+ NB: relative to other versions, 'score' is multiplied by 2
+ (and is an integer)
+
+ $Revision: 1.7 $ $Date: 2013/09/18 04:20:13 $
+
+ */
+
+#include <R_ext/Utils.h>
+#include "chunkloop.h"
+
+void inxyp(x,y,xp,yp,npts,nedges,score,onbndry)
+ /* inputs */
+ double *x, *y; /* points to be tested */
+ int *npts;
+ double *xp, *yp; /* polygon vertices */
+ int *nedges;
+ /* outputs */
+ int *score;
+ int *onbndry;
+{
+ int i, j, Npts, Nedges, Ne1, contrib, maxchunk;
+ double x0, y0, x1, y1, dx, dy, xj, yj, xcrit, ycrit;
+
+ Npts = *npts;
+ Nedges = *nedges;
+ Ne1 = Nedges - 1;
+
+ x0 = xp[Ne1];
+ y0 = yp[Ne1];
+
+ OUTERCHUNKLOOP(i, Nedges, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Nedges, maxchunk, 16384) {
+ /* visit edge (x0,y0) -> (x1,y1) */
+ x1 = xp[i];
+ y1 = yp[i];
+ dx = x1 - x0;
+ dy = y1 - y0;
+ for(j = 0; j < Npts; j++) {
+ xj = x[j];
+ yj = y[j];
+ xcrit = (xj - x0) * (xj - x1);
+ if(xcrit <= 0) {
+ if(xcrit == 0) {
+ contrib = 1;
+ } else {
+ contrib = 2;
+ }
+ ycrit = yj * dx - xj * dy + x0 * dy - y0 * dx;
+ if(dx < 0) {
+ if(ycrit >= 0)
+ score[j] += contrib;
+ onbndry[j] = onbndry[j] | (ycrit == 0);
+ } else if(dx > 0) {
+ if(ycrit < 0)
+ score[j] -= contrib;
+ onbndry[j] = onbndry[j] | (ycrit == 0);
+ } else {
+ if(xj == x0)
+ ycrit = (yj - y0) * (yj - y1);
+ onbndry[j] = onbndry[j] | (ycrit <= 0);
+ }
+ }
+ }
+ /* next edge */
+ x0 = x1;
+ y0 = y1;
+ }
+ }
+}
diff --git a/src/matchindices.c b/src/matchindices.c
new file mode 100644
index 0000000..ed93bd2
--- /dev/null
+++ b/src/matchindices.c
@@ -0,0 +1,197 @@
+/*
+
+ matchindices.c
+
+ $Revision$ $Date$
+
+ CSmatch2int Find matches between two sorted lists of pairs of integers
+ CSmatch3int Find matches between two sorted lists of triples of integers
+
+ CUmatch2int Find matches between two UNsorted lists of pairs of integers
+ CUmatch3int Find matches between two UNsorted lists of triples of integers
+
+*/
+
+#include <math.h>
+#include <R_ext/Utils.h>
+#include "chunkloop.h"
+
+/* ................ unsorted ............................ */
+/* ........ Behaviour equivalent to match() ............. */
+
+
+/*
+
+ CUmatch2int
+
+ Find matches between two unsorted lists of pairs of integers
+
+ */
+
+void CUmatch2int(na, xa, ya, nb, xb, yb, match)
+ /* inputs */
+ int *na, *nb;
+ int *xa, *ya;
+ int *xb, *yb;
+ /* output */
+ int *match;
+ /* match[i] = j+1 if xb[j], yb[j] matches xa[i], ya[i] */
+ /* match[i] = 0 if no such point matches xa[i], ya[i] */
+{
+ int i, j, Na, Nb, maxchunk;
+ int xai, yai;
+
+ Na = *na;
+ Nb = *nb;
+
+ OUTERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ xai = xa[i];
+ yai = ya[i];
+ match[i] = 0;
+ for(j = 0; j < Nb; j++) {
+ if(xb[j] == xai && yb[j] == yai) {
+ match[i] = j+1;
+ break;
+ }
+ }
+ }
+ }
+}
+
+
+/*
+
+ CUmatch3int
+
+ Find matches between two unsorted lists of triples of integers
+
+ */
+
+void CUmatch3int(na, xa, ya, za, nb, xb, yb, zb, match)
+ /* inputs */
+ int *na, *nb;
+ int *xa, *ya, *za;
+ int *xb, *yb, *zb;
+ /* output */
+ int *match;
+ /* match[i] = j+1 if xb[j], yb[j], zb[j] matches xa[i], ya[i], za[i] */
+ /* match[i] = 0 if no such point matches xa[i], ya[i], za[i] */
+{
+ int i, j, Na, Nb, maxchunk;
+ int xai, yai, zai;
+
+ Na = *na;
+ Nb = *nb;
+
+ j = 0;
+
+ OUTERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ xai = xa[i];
+ yai = ya[i];
+ zai = za[i];
+ match[i] = 0;
+ for(j = 0; j < Nb; j++) {
+ if(xb[j] == xai && yb[j] == yai && zb[j] == zai) {
+ match[i] = j+1;
+ break;
+ }
+ }
+ }
+ }
+}
+
+
+/* ................ sorted ............................ */
+
+/*
+
+ CSmatch2int
+
+ Find matches between two lists of pairs of integers
+
+ Each list sorted by order(x,y)
+
+ */
+
+void CSmatch2int(na, xa, ya, nb, xb, yb, match)
+ /* inputs */
+ int *na, *nb;
+ int *xa, *ya; /* sorted into increasing order of (xa, ya) */
+ int *xb, *yb; /* sorted into increasing order of (xb, yb) */
+ /* output */
+ int *match;
+ /* match[i] = j+1 if xb[j], yb[j] matches xa[i], ya[i] */
+ /* match[i] = 0 if no such point matches xa[i], ya[i] */
+{
+ int i, j, Na, Nb, maxchunk;
+ int xai, yai;
+
+ Na = *na;
+ Nb = *nb;
+
+ j = 0;
+
+ OUTERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ xai = xa[i];
+ yai = ya[i];
+ match[i] = 0;
+ while(j < Nb && xb[j] < xai) ++j;
+ while(j < Nb && xb[j] == xai && yb[j] < yai) ++j;
+ if(j < Nb && xb[j] == xai && yb[j] == yai)
+ match[i] = j+1;
+ if(j >= Nb) return;
+ }
+ }
+}
+
+
+/*
+
+ CSmatch3int
+
+ Find matches between two lists of triples of integers
+
+ Each list sorted by order(x,y,z)
+
+ */
+
+void CSmatch3int(na, xa, ya, za, nb, xb, yb, zb, match)
+ /* inputs */
+ int *na, *nb;
+ int *xa, *ya, *za; /* sorted into increasing order of (xa, ya, za) */
+ int *xb, *yb, *zb; /* sorted into increasing order of (xb, yb, zb) */
+ /* output */
+ int *match;
+ /* match[i] = j+1 if xb[j], yb[j], zb[j] matches xa[i], ya[i], za[i] */
+ /* match[i] = 0 if no such point matches xa[i], ya[i], za[i] */
+{
+ int i, j, Na, Nb, maxchunk;
+ int xai, yai, zai;
+
+ Na = *na;
+ Nb = *nb;
+
+ j = 0;
+
+ OUTERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ xai = xa[i];
+ yai = ya[i];
+ zai = za[i];
+ match[i] = 0;
+ while(j < Nb && xb[j] < xai) ++j;
+ while(j < Nb && xb[j] == xai && yb[j] < yai) ++j;
+ while(j < Nb && xb[j] == xai && yb[j] == yai && zb[j] < zai) ++j;
+ if(j < Nb && xb[j] == xai && yb[j] == yai && zb[j] == zai)
+ match[i] = j+1;
+ if(j >= Nb) return;
+ }
+ }
+}
diff --git a/src/matchpoints.c b/src/matchpoints.c
new file mode 100644
index 0000000..ba82bd9
--- /dev/null
+++ b/src/matchpoints.c
@@ -0,0 +1,52 @@
+/*
+
+ matchpoints.c
+
+ $Revision: 1.1 $ $Date: 2017/01/08 00:32:43 $
+
+ Cmatchxy Find matches between two sets of points
+
+*/
+
+#include <math.h>
+#include <R_ext/Utils.h>
+#include "chunkloop.h"
+
+/*
+
+ Cmatchxy
+
+ Find matches between two lists of points
+
+ */
+
+void Cmatchxy(na, xa, ya, nb, xb, yb, match)
+ /* inputs */
+ int *na, *nb;
+ double *xa, *ya, *xb, *yb;
+ /* output */
+ int *match;
+ /* match[i] = j+1 if xb[j], yb[j] matches xa[i], ya[i] */
+ /* match[i] = 0 if no such point matches xa[i], ya[i] */
+{
+ int i, j, Na, Nb, maxchunk;
+ double xai, yai;
+
+ Na = *na;
+ Nb = *nb;
+
+ OUTERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ R_CheckUserInterrupt();
+ INNERCHUNKLOOP(i, Na, maxchunk, 16384) {
+ xai = xa[i];
+ yai = ya[i];
+ match[i] = 0;
+ for (j=0; j < Nb; j++) {
+ if(xai == xb[j] && yai == yb[j]) {
+ match[i] = j+1;
+ break;
+ }
+ }
+ }
+ }
+}
diff --git a/src/ply.c b/src/ply.c
new file mode 100644
index 0000000..496bf0a
--- /dev/null
+++ b/src/ply.c
@@ -0,0 +1,31 @@
+/*
+ ply.c
+
+ Faster versions of tapply(..., FUN=sum)
+ assuming indices are sorted.
+
+ Code template is in 'ply.h'
+
+ Adrian Baddeley and Tilman Davies
+
+ $Revision: 1.2 $ $Date: 2016/08/15 02:34:03 $
+
+*/
+
+#include <R.h>
+#include <R_ext/Utils.h>
+#include <Rmath.h>
+
+#define FNAME ply3sum
+#define NDIM 3
+#include "ply.h"
+#undef FNAME
+#undef NDIM
+
+#define FNAME ply2sum
+#define NDIM 2
+#include "ply.h"
+#undef FNAME
+#undef NDIM
+
+
diff --git a/src/ply.h b/src/ply.h
new file mode 100644
index 0000000..c116d24
--- /dev/null
+++ b/src/ply.h
@@ -0,0 +1,85 @@
+/*
+ ply.h
+
+ Template for functions in ply.c
+ This file is #included several times
+
+ Macros used:
+ FNAME Name of C routine
+ NDIM Number of dimensions of result (2 or 3)
+
+ Adrian Baddeley and Tilman Davies
+
+ $Revision: 1.1 $ $Date: 2016/08/15 02:29:15 $
+
+*/
+
+
+void FNAME(nin,
+ xin,
+ iin,
+ jin,
+#if (NDIM > 2)
+ kin,
+#endif
+ nout,
+ xout,
+ iout,
+ jout
+#if (NDIM > 2)
+ , kout
+#endif
+)
+ int *nin, *nout;
+ double *xin, *xout;
+ int *iin, *jin, *iout, *jout;
+#if (NDIM > 2)
+ int *kin, *kout;
+#endif
+{
+ int Nin, l, m, icur, jcur;
+#if (NDIM > 2)
+ int kcur;
+#endif
+ double xsum;
+ Nin = *nin;
+ if(Nin == 0) {
+ *nout = 0;
+ return;
+ }
+ /* initialise first cell using first entry */
+ m = 0;
+ iout[0] = icur = iin[0];
+ jout[0] = jcur = jin[0];
+#if (NDIM > 2)
+ kout[0] = kcur = kin[0];
+#endif
+ xout[0] = xsum = xin[0];
+ /* process subsequent entries */
+ if(Nin > 1) {
+ for(l = 1; l < Nin; l++) {
+ if(iin[l] == icur && jin[l] == jcur
+#if (NDIM > 2)
+ && kin[l] == kcur
+#endif
+ ) {
+ /* increment current sum */
+ xsum += xin[l];
+ } else {
+ /* write cell result */
+ xout[m] = xsum;
+ /* initialise next cell */
+ ++m;
+ iout[m] = icur = iin[l];
+ jout[m] = jcur = jin[l];
+#if (NDIM > 2)
+ kout[m] = kcur = kin[l];
+#endif
+ xsum = xin[l];
+ }
+ /* write last cell */
+ xout[m] = xsum;
+ }
+ }
+ *nout = m + 1;
+}
diff --git a/src/primefax.c b/src/primefax.c
new file mode 100644
index 0000000..c96b3c3
--- /dev/null
+++ b/src/primefax.c
@@ -0,0 +1,175 @@
+/*
+
+ primefax.c
+
+ Prime numbers and prime factorisation
+
+ $Revision: 1.2 $ $Date: 2016/12/31 08:40:29 $
+
+*/
+
+#include <R.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+
+
+int primetable[] = {
+ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
+ 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
+ 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191,
+ 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
+ 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
+ 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
+ 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
+ 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
+ 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
+ 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757,
+ 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853,
+ 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
+ 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021,
+ 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,
+ 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
+ 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259,
+ 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
+ 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433,
+ 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
+ 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
+ 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
+ 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
+ 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831,
+ 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913,
+ 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
+ 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087,
+ 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
+ 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269,
+ 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347,
+ 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,
+ 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
+ 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
+ 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
+ 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767,
+ 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
+ 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953,
+ 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041,
+ 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163,
+ 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251,
+ 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329,
+ 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
+ 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517,
+ 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
+ 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673,
+ 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767,
+ 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853,
+ 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931,
+ 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027,
+ 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
+ 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229,
+ 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
+ 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421,
+ 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513,
+ 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603,
+ 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691,
+ 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793,
+ 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,
+ 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987,
+ 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077,
+ 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171,
+ 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
+ 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393,
+ 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,
+ 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557,
+ 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,
+ 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741,
+ 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839,
+ 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903,
+ 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043,
+ 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131,
+ 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
+ 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311,
+ 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,
+ 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521,
+ 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607,
+ 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703,
+ 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803,
+ 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899,
+ 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983,
+ 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079,
+ 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
+ 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307,
+ 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433,
+ 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523,
+ 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589,
+ 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687,
+ 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
+ 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883,
+ 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009,
+ 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101,
+ 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191,
+ 0 /* end marker */
+};
+
+#define PMAX 8192
+
+#undef BUGGY
+
+void primefax(n, factors, nfactors)
+ int *n;
+ int *factors;
+ int *nfactors;
+{
+ int m, p, dmax, k, d;
+ int *ptr;
+
+ m = *n;
+ k = 0;
+ /* upper limit on divisors */
+ dmax = (int) ceil(sqrt((double) m));
+
+#ifdef BUGGY
+ Rprintf("n = %d, dmax=%d\n", m, dmax);
+#endif
+
+ /* search for prime divisors in table of primes */
+ for(ptr = primetable; *ptr != 0; ++ptr) {
+ p = *ptr;
+#ifdef BUGGY
+ Rprintf("m = %d, p = %d\n", m, p);
+#endif
+ while(m % p == 0) {
+#ifdef BUGGY
+ Rprintf("\tdivides!\n");
+#endif
+ factors[k] = p;
+ ++k;
+ m = m/p;
+ }
+ if(p > m || p > dmax) break;
+ }
+
+ if(m > 1 && dmax > PMAX) {
+ /* search for divisors above PMAX */
+#ifdef BUGGY
+ Rprintf("searching beyond table..\n");
+#endif
+ for(d = PMAX; d * d <= m; ++d) {
+ while(m % d == 0) {
+ factors[k] = d;
+ k++;
+ m = m/d;
+ }
+ if(d > dmax) break;
+ }
+ }
+
+ if(m == 1) {
+ /* n has been factorised */
+ *nfactors = k;
+ return;
+ }
+
+ /* m is prime */
+ factors[k] = m;
+ ++k;
+ *nfactors = k;
+ return;
+}
diff --git a/src/proto.h b/src/proto.h
new file mode 100644
index 0000000..2d55c8d
--- /dev/null
+++ b/src/proto.h
@@ -0,0 +1,37 @@
+#include <R.h>
+#include <Rinternals.h>
+
+/*
+ Prototype declarations for all native routines in spatstat.utils package
+
+ Automatically generated - do not edit!
+
+*/
+
+/*
+
+ Functions invoked by .C
+
+*/
+
+void primefax(int *, int *, int *);
+void ply2sum(int *, double *, int *, int *, int *, double *, int *, int *);
+void ply3sum(int *, double *, int *, int *, int *, int *, double *, int *, int *, int *);
+void CSmatch2int(int *, int *, int *, int *, int *, int *, int *);
+void CUmatch2int(int *, int *, int *, int *, int *, int *, int *);
+void CSmatch3int(int *, int *, int *, int *, int *, int *, int *, int *, int *);
+void CUmatch3int(int *, int *, int *, int *, int *, int *, int *, int *, int *);
+void irevcumsum(int *, int *); void drevcumsum(double *, int *);
+void fastinterv(double *, int *, double *, int *, int *);
+void Cmatchxy(int *, double *, double *, int *, double *, double *, int *);
+void inxyp(double *, double *, double *, double *, int *, int *, int *, int *);
+void prdist2segs(double *, double *, int *, double *, double *, double *, double *, int *, double *, double *);
+void nndist2segs(double *, double *, int *, double *, double *, double *, double *, int *, double *, double *, int *);
+/*
+
+ Functions invoked by .Call
+
+*/
+SEXP circXseg(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+SEXP circMseg(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+SEXP circXseg(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
diff --git a/src/revcum.c b/src/revcum.c
new file mode 100755
index 0000000..c96addf
--- /dev/null
+++ b/src/revcum.c
@@ -0,0 +1,41 @@
+/*
+
+ revcum.c
+
+ $Revision: 1.3 $ $Date: 2016/12/30 07:28:13 $
+
+ Reverse cumulative sums
+
+*/
+
+void drevcumsum(double *x, int *nx) {
+ int i;
+ double sumx;
+ double *xp;
+
+ i = *nx - 1;
+ xp = x + i;
+ sumx = *xp;
+ while(i > 0) {
+ --i;
+ --xp;
+ sumx += *xp;
+ *xp = sumx;
+ }
+}
+
+void irevcumsum(int *x, int *nx) {
+ int i;
+ int sumx;
+ int *xp;
+
+ i = *nx - 1;
+ xp = x + i;
+ sumx = *xp;
+ while(i > 0) {
+ --i;
+ --xp;
+ sumx += *xp;
+ *xp = sumx;
+ }
+}
diff --git a/tests/numerical.R b/tests/numerical.R
new file mode 100644
index 0000000..d0eedee
--- /dev/null
+++ b/tests/numerical.R
@@ -0,0 +1,25 @@
+#' spatstat.utils/tests/numerical.R
+#' Tests of numerical code
+
+require(spatstat.utils)
+
+#' validity of 'tapplysum'
+aa <- factor(letters[1:3])
+bb <- factor(letters[1:4])[c(1,2,2)]
+xx <- round(runif(3), 3)
+yy <- tapplysum(xx, list(A=aa, B=bb), do.names=TRUE)
+zz <- tapply(xx, list(A=aa, B=bb), sum)
+zz[is.na(zz)] <- 0
+if(any(yy != zz))
+ stop("tapplysum does not agree with tapply(, sum)")
+#' tapplysum with zero-length data
+tapplysum(xx[FALSE], list(A=aa[FALSE], B=bb[FALSE]), do.names=TRUE)
+
+#' validity of matchIntegerDataFrames
+
+A <- data.frame(a=sample(1:5), b=sample(1:5, replace=TRUE))
+B <- data.frame(u=sample(1:3), w=3:1)
+acode <- paste(A[,1], A[,2])
+bcode <- paste(B[,1], B[,2])
+stopifnot(identical(matchIntegerDataFrames(A,B), match(acode,bcode)))
+
diff --git a/tests/segments.R b/tests/segments.R
new file mode 100644
index 0000000..53ba229
--- /dev/null
+++ b/tests/segments.R
@@ -0,0 +1,13 @@
+#' spatstat.utils/tests/segments.R
+
+require(spatstat.utils)
+
+#' test of distppll pointed out by Ang Qi Wei
+
+p <- matrix(c(1.5, 0), 1, 2)
+l <- matrix(c(0,0,1,0,1,0,2,0), 2, 4, byrow=T)
+a <- distppll(p, l, mintype=2, method="interpreted")
+d <- distppll(p, l, mintype=2, method="C")
+if(a$min.which != d$min.which)
+ stop("conflict between C and interpreted code in distppll")
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-spatstat.utils.git
More information about the debian-science-commits
mailing list