[r-cran-polycub] 11/18: New upstream version 0.6.0
Andreas Tille
tille at debian.org
Fri Oct 20 14:21:49 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-polycub.
commit 8fa29d0d1575dda4be35d929ed1578603ee5e2a9
Author: Andreas Tille <tille at debian.org>
Date: Fri Oct 20 16:06:20 2017 +0200
New upstream version 0.6.0
---
DESCRIPTION | 42 ++++
MD5 | 53 +++++
NAMESPACE | 52 +++++
R/circleCub.R | 52 +++++
R/coerce-gpc-methods.R | 116 +++++++++++
R/coerce-sp-methods.R | 93 +++++++++
R/plotpolyf.R | 86 +++++++++
R/polyCub.R | 50 +++++
R/polyCub.SV.R | 352 ++++++++++++++++++++++++++++++++++
R/polyCub.exact.Gauss.R | 217 +++++++++++++++++++++
R/polyCub.iso.R | 252 ++++++++++++++++++++++++
R/polyCub.midpoint.R | 80 ++++++++
R/sysdata.rda | Bin 0 -> 18852 bytes
R/tools.R | 98 ++++++++++
R/xylist.R | 124 ++++++++++++
R/zzz.R | 115 +++++++++++
README.md | 36 ++++
debian/README.source | 8 -
debian/README.test | 21 --
debian/changelog | 36 ----
debian/compat | 1 -
debian/control | 33 ----
debian/copyright | 29 ---
debian/docs | 3 -
debian/patches/series | 1 -
debian/patches/skip_gpclib_test.patch | 25 ---
debian/rules | 4 -
debian/source/format | 1 -
debian/tests/control | 3 -
debian/tests/run-unit-test | 12 --
debian/watch | 2 -
inst/CITATION | 22 +++
inst/NEWS.Rd | 248 ++++++++++++++++++++++++
inst/examples/plotpolyf.R | 11 ++
inst/examples/polyCub.R | 55 ++++++
inst/examples/polyCub.iso.R | 22 +++
inst/include/polyCubAPI.h | 37 ++++
man/checkintrfr.Rd | 45 +++++
man/circleCub.Gauss.Rd | 51 +++++
man/coerce-gpc-methods.Rd | 54 ++++++
man/coerce-sp-methods.Rd | 43 +++++
man/dotprod.Rd | 19 ++
man/gpclibPermit.Rd | 18 ++
man/isClosed.Rd | 21 ++
man/isScalar.Rd | 18 ++
man/makegrid.Rd | 23 +++
man/plot_polyregion.Rd | 23 +++
man/plotpolyf.Rd | 68 +++++++
man/polyCub-package.Rd | 73 +++++++
man/polyCub.Rd | 103 ++++++++++
man/polyCub.SV.Rd | 97 ++++++++++
man/polyCub.exact.Gauss.Rd | 76 ++++++++
man/polyCub.iso.Rd | 134 +++++++++++++
man/polyCub.midpoint.Rd | 61 ++++++
man/polygauss.Rd | 56 ++++++
man/vecnorm.Rd | 19 ++
man/xylist.Rd | 87 +++++++++
src/init.c | 29 +++
src/polyCub.SV.c | 85 ++++++++
src/polyCub.SV.h | 18 ++
src/polyCub.iso.c | 140 ++++++++++++++
src/polyCub.iso.h | 21 ++
tests/test-all.R | 3 +
tests/testthat/polyiso_powerlaw.c | 46 +++++
tests/testthat/test-NWGL.R | 9 +
tests/testthat/test-polyCub.R | 59 ++++++
tests/testthat/test-polyiso.R | 103 ++++++++++
tests/testthat/test-regression.R | 12 ++
68 files changed, 3827 insertions(+), 179 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..6610b67
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,42 @@
+Package: polyCub
+Title: Cubature over Polygonal Domains
+Version: 0.6.0
+Date: 2017-05-24
+Authors at R: c(
+ person("Sebastian", "Meyer",
+ email = "seb.meyer at fau.de",
+ role = c("aut","cre","trl")),
+ person("Leonhard", "Held",
+ email = "Leonhard.Held at uzh.ch",
+ role = "ths"),
+ person("Michael", "Hoehle",
+ email = "hoehle at math.su.se",
+ role = "ths")
+ )
+Description: The following methods for cubature (numerical integration)
+ over polygonal domains are currently implemented:
+ the two-dimensional midpoint rule as a simple wrapper around
+ as.im.function() from package 'spatstat' (Baddeley and Turner, 2005),
+ the product Gauss cubature by Sommariva and Vianello (2007),
+ an adaptive cubature for isotropic functions via line integrate()
+ along the boundary (Meyer and Held, 2014),
+ and quasi-exact methods specific to the integration of the
+ bivariate Gaussian density over polygonal and circular domains
+ (based on formulae from the Abramowitz and Stegun (1972) handbook).
+ For cubature over simple hypercubes, the packages 'cubature' and
+ 'R2Cuba' are more appropriate.
+License: GPL-2
+URL: https://github.com/WastlM/polyCub
+BugReports: https://github.com/WastlM/polyCub/issues
+Depends: R (>= 2.15.0), methods, sp
+Imports: grDevices, graphics, stats, spatstat
+Suggests: lattice, testthat, mvtnorm, statmod, rgeos, gpclib
+RoxygenNote: 6.0.1
+NeedsCompilation: yes
+Packaged: 2017-05-26 09:23:56 UTC; smeyer
+Author: Sebastian Meyer [aut, cre, trl],
+ Leonhard Held [ths],
+ Michael Hoehle [ths]
+Maintainer: Sebastian Meyer <seb.meyer at fau.de>
+Repository: CRAN
+Date/Publication: 2017-05-26 10:05:58 UTC
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..2f936bf
--- /dev/null
+++ b/MD5
@@ -0,0 +1,53 @@
+d4cd408a6dda4d4fd832fad0adcf4ebc *DESCRIPTION
+ebce3d96440b5368513e2d434be16a24 *NAMESPACE
+0a60a7ee72d69596c5c42d0e8070e42b *R/circleCub.R
+6c81138181a8c99ed924e298dec063d2 *R/coerce-gpc-methods.R
+866d2c71b8fdf86f02a1253cade4dbf4 *R/coerce-sp-methods.R
+671d2d34604542e536dea2027789228b *R/plotpolyf.R
+00f8503508c8fe19abb7e261193adff8 *R/polyCub.R
+d26a5f10f78a832a3dfd00d4e2b8e043 *R/polyCub.SV.R
+e352b88e0860a77b7df505008d722400 *R/polyCub.exact.Gauss.R
+6f6287cf77cd2e15bafbf5916dddeba1 *R/polyCub.iso.R
+33b9219de529db39375a1e04e593cea4 *R/polyCub.midpoint.R
+75a1bb73253360d3ffdbb005b55f083f *R/sysdata.rda
+482735f9ca96e14ff8963614b524d63e *R/tools.R
+331b633713ef98289a6f462cc968c97c *R/xylist.R
+a9f5580a4daa787b5c36804c61336084 *R/zzz.R
+cae095e970caf0de622c4d4fc65fe38e *README.md
+915021f57341348a5f9348605269001e *inst/CITATION
+ade60c55007075ce440ccdd6ec24323b *inst/NEWS.Rd
+a2e8c8c02633c62187d4b62c80ed7bde *inst/examples/plotpolyf.R
+71c354b04903463684d7ee0c21afc397 *inst/examples/polyCub.R
+fec41191e00776d1cb4615e89d6e3612 *inst/examples/polyCub.iso.R
+fa7799fccf63a1561bb2ed114a2f07db *inst/include/polyCubAPI.h
+48f7672ddb81a77f3793f2b26f4f9cbe *man/checkintrfr.Rd
+80c46e828a042699ab4e7500baf363f5 *man/circleCub.Gauss.Rd
+14855b08e09754316834a834f4dc2d22 *man/coerce-gpc-methods.Rd
+2b9885fd97c096fd4f553c35a4bb119f *man/coerce-sp-methods.Rd
+323f56fe3449f60180325c6263454377 *man/dotprod.Rd
+c38994de2c8de6f2372c42eb4ce915a4 *man/gpclibPermit.Rd
+db65b0e51cfcb8f175e6a32dc5c495e6 *man/isClosed.Rd
+b31e8ec92525bba018a7368f2d270ce2 *man/isScalar.Rd
+31f03498d486eeaef7ba89fec69a3415 *man/makegrid.Rd
+6f71baeb60ea79b309dd533de9bc3f0f *man/plot_polyregion.Rd
+4b162374f6661084873f7566033465be *man/plotpolyf.Rd
+414657ab72d2336196a2ae3dde718c2a *man/polyCub-package.Rd
+46a1bcad91a7041eb29af54604ddd2d7 *man/polyCub.Rd
+cebb825e063d60ddfa10ad67d39d18ea *man/polyCub.SV.Rd
+79c8a789e697e7afc5b148dd32d2c99a *man/polyCub.exact.Gauss.Rd
+74681dcd35af645008bb1e07239283a4 *man/polyCub.iso.Rd
+33f30984750652e060b59fd0a9fd54b1 *man/polyCub.midpoint.Rd
+05a598d51ea290a90cf6ba4cfab5554f *man/polygauss.Rd
+3e57be2294c537ed441d6273872bd716 *man/vecnorm.Rd
+f2aa40b3e8ed1379ae7a2d51fd7d2eb9 *man/xylist.Rd
+17ed5a8294f4ffd649be6e4d710e9907 *src/init.c
+a18f7041fa9c917423f88b60fd000384 *src/polyCub.SV.c
+c88a07495cbb628e0bee68be96c46b5f *src/polyCub.SV.h
+da53b0e4064c5cf355e7a3eb37fe3819 *src/polyCub.iso.c
+9b7e40e6ba862b9d54817793c3827ec4 *src/polyCub.iso.h
+3e4e9e53ad2f425939503fd77e030cd2 *tests/test-all.R
+169ef46029292ff8bfae48ba27834b8c *tests/testthat/polyiso_powerlaw.c
+b56f6ab7f21b2441ccd03873d4a1605c *tests/testthat/test-NWGL.R
+4076eacdcf7edcf0e376ecec57521d3f *tests/testthat/test-polyCub.R
+4222d51ae8cc74424b8fbf735c1bb55c *tests/testthat/test-polyiso.R
+f60d877c5efbec6d36927cfd1fd2a4d5 *tests/testthat/test-regression.R
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..2703a8d
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,52 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(as.owin,Polygon)
+S3method(as.owin,Polygons)
+S3method(as.owin,SpatialPolygons)
+S3method(as.owin,gpc.poly)
+S3method(xylist,Polygon)
+S3method(xylist,Polygons)
+S3method(xylist,SpatialPolygons)
+S3method(xylist,default)
+S3method(xylist,gpc.poly)
+S3method(xylist,owin)
+export(.polyCub.iso)
+export(checkintrfr)
+export(circleCub.Gauss)
+export(gpc2owin)
+export(gpclibPermit)
+export(gpclibPermitStatus)
+export(owin2gpc)
+export(plotpolyf)
+export(polyCub)
+export(polyCub.SV)
+export(polyCub.exact.Gauss)
+export(polyCub.iso)
+export(polyCub.midpoint)
+export(xylist)
+exportMethods(coerce)
+import(methods)
+import(sp)
+importFrom(grDevices,extendrange)
+importFrom(grDevices,gray)
+importFrom(grDevices,heat.colors)
+importFrom(grDevices,xy.coords)
+importFrom(graphics,image)
+importFrom(graphics,lines)
+importFrom(graphics,points)
+importFrom(graphics,polygon)
+importFrom(spatstat,as.im.function)
+importFrom(spatstat,as.owin)
+importFrom(spatstat,as.polygonal)
+importFrom(spatstat,integral.im)
+importFrom(spatstat,is.polygonal)
+importFrom(spatstat,owin)
+importFrom(spatstat,plot.im)
+importFrom(spatstat,plot.owin)
+importFrom(spatstat,summary.owin)
+importFrom(stats,cov2cor)
+importFrom(stats,dist)
+importFrom(stats,integrate)
+importFrom(stats,pchisq)
+importFrom(stats,pnorm)
+useDynLib(polyCub, .registration = TRUE)
diff --git a/R/circleCub.R b/R/circleCub.R
new file mode 100644
index 0000000..448ca96
--- /dev/null
+++ b/R/circleCub.R
@@ -0,0 +1,52 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2013-2014 Sebastian Meyer
+### Time-stamp: <[circleCub.R] by SM Die 06/05/2014 10:02 (CEST)>
+###
+### Special cases of cubature over circular domains (center, r)
+################################################################################
+
+
+##' Integration of the Isotropic Gaussian Density over Circular Domains
+##'
+##' This function calculates the integral of the bivariate, isotropic Gaussian
+##' density (i.e. \eqn{\Sigma} = \code{sd^2*diag(2)}) over circular domains via
+##' the cumulative distribution function of the (non-central) Chi-Squared
+##' distribution (\code{pchisq}), cp. Formula 26.3.24 in Abramowitz and Stegun
+##' (1972).
+##'
+##' @references
+##' Abramowitz, M. and Stegun, I. A. (1972).
+##' Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+##' Tables. New York: Dover Publications.
+##' @param center numeric vector of length 2 (center of the circle).
+##' @param r numeric (radius of the circle). Several radii may be supplied.
+##' @param mean numeric vector of length 2
+##' (mean of the bivariate Gaussian density).
+##' @param sd numeric (common standard deviation of the isotropic
+##' Gaussian density in both dimensions).
+##' @return The integral value (one for each supplied radius).
+##' @note The non-centrality parameter of the evaluated chi-squared distribution
+##' equals the squared distance between the \code{mean} and the
+##' \code{center}. If this becomes too large, the result becomes inaccurate, see
+##' \code{\link{pchisq}}.
+##' @keywords math spatial
+##' @importFrom stats pchisq
+##' @export
+##' @examples
+##' circleCub.Gauss(center=c(1,2), r=3, mean=c(4,5), sd=6)
+##'
+##' if (requireNamespace("mvtnorm") && gpclibPermit()) {
+##' ## compare with cubature over a polygonal approximation of a circle
+##' disc.poly <- spatstat::disc(radius=3, centre=c(1,2), npoly=32)
+##' polyCub.exact.Gauss(disc.poly, mean=c(4,5), Sigma=6^2*diag(2))
+##' }
+
+circleCub.Gauss <- function (center, r, mean, sd)
+{
+ stopifnot(isScalar(sd), length(center) == 2, length(mean) == 2)
+ pchisq((r/sd)^2, df=2, ncp=sum(((center-mean)/sd)^2))
+}
diff --git a/R/coerce-gpc-methods.R b/R/coerce-gpc-methods.R
new file mode 100644
index 0000000..fe5e23c
--- /dev/null
+++ b/R/coerce-gpc-methods.R
@@ -0,0 +1,116 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2012-2015 Sebastian Meyer
+### Time-stamp: <[coerce-gpc-methods.R] 2015-02-25 21:07 (CET) by SM>
+################################################################################
+
+
+##' Conversion between polygonal \code{"owin"} and \code{"gpc.poly"}
+##'
+##' Package \pkg{polyCub} implements converters between the classes
+##' \code{"\link[=owin.object]{owin}"} of package \pkg{spatstat} and
+##' \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} of package \pkg{rgeos}
+##' (originally from \pkg{gpclib}).
+##' Support for the \code{"gpc.poly"} class was dropped from
+##' \pkg{spatstat} as of version 1.34-0.
+##'
+##' @param object an object of class \code{"gpc.poly"} or \code{"owin"},
+##' respectively.
+##' @return The converted polygon of class \code{"gpc.poly"} or \code{"owin"},
+##' respectively. If neither package \pkg{rgeos} nor \pkg{gpclib} are available,
+##' \code{owin2gpc} will just return the \code{pts} slot of the
+##' \code{"gpc.poly"} (no formal class) with a warning.
+##' @author Sebastian Meyer
+##' @note The converter \code{owin2gpc} requires the package \pkg{rgeos} (or
+##' \pkg{gpclib}) for the formal class definition of a \code{"gpc.poly"}.
+##' It will produce vertices ordered according to the \pkg{sp} convention,
+##' i.e. clockwise for normal boundaries and anticlockwise for holes, where,
+##' however, the first vertex is \emph{not} repeated!
+##' @seealso \code{\link{xylist}}, and the package \pkg{rgeos} for
+##' conversions of \code{"gpc.poly"} objects from and to \pkg{sp}'s
+##' \code{"\linkS4class{SpatialPolygons}"} class.
+##' @name coerce-gpc-methods
+##' @rdname coerce-gpc-methods
+##' @keywords spatial methods
+##' @importFrom spatstat as.polygonal summary.owin
+##' @import methods
+##' @export
+owin2gpc <- function (object)
+{
+ object <- as.polygonal(object)
+
+ ## FIXME: spatstat functions to extract the areas and hole flags
+ ## of the individual polygons in a polygonal "owin" would be nice
+ holes <- summary.owin(object)$areas < 0
+
+ ## reverse vertices and set hole flags
+ pts <- mapply(
+ FUN = function (poly, hole) {
+ list(x = rev(poly$x), y = rev(poly$y), hole = hole)
+ ## or hole = area.owin(owin(poly = poly, check = FALSE)) < 0,
+ ## but spatstat::is.hole.xypolygon is marked as an internal function
+ },
+ poly = object$bdry, hole = holes,
+ SIMPLIFY = FALSE, USE.NAMES = FALSE)
+
+ ## formal class
+ if (know_gpc.poly()) {
+ new("gpc.poly", pts = pts)
+ } else {
+ warning("formal class \"gpc.poly\" not available")
+ pts
+ }
+}
+
+##' @inheritParams owin2gpc
+##' @param ... further arguments passed to \code{\link{owin}}.
+##' @rdname coerce-gpc-methods
+##' @importFrom spatstat owin summary.owin
+##' @export
+gpc2owin <- function (object, ...)
+{
+ ## first convert to an "owin" without checking areas etc.
+ ## to determine the hole status according to vertex order (area)
+ res <- owin(poly = object at pts, check = FALSE)
+ holes_owin <- summary.owin(res)$areas < 0
+ ## Note: cannot rely on spatstat::Area.xypolygon since it is marked internal
+
+ ## now fix the vertex order
+ bdry <- mapply(
+ FUN = function (poly, owinhole) {
+ if (poly$hole != owinhole) {
+ poly$x <- rev(poly$x)
+ poly$y <- rev(poly$y)
+ }
+ poly
+ },
+ poly = object at pts, owinhole = holes_owin,
+ SIMPLIFY = FALSE, USE.NAMES = FALSE)
+
+ ## now really convert to owin with appropriate vertex order
+ owin(poly = bdry, ...)
+}
+
+##' @inheritParams gpc2owin
+##' @param W an object of class \code{"gpc.poly"}.
+##' @rdname coerce-gpc-methods
+##' @importFrom spatstat as.owin
+##' @method as.owin gpc.poly
+##' @export
+as.owin.gpc.poly <- function (W, ...)
+{
+ gpc2owin(W, ...)
+}
+
+
+## check for the formal class "gpc.poly" (loading rgeos or gpclib if necessary)
+##' @import methods
+know_gpc.poly <- function ()
+{
+ isClass("gpc.poly") ||
+ suppressWarnings(requireNamespace("rgeos", quietly=TRUE) ||
+ requireNamespace("gpclib", quietly=TRUE))
+}
diff --git a/R/coerce-sp-methods.R b/R/coerce-sp-methods.R
new file mode 100644
index 0000000..e3040ea
--- /dev/null
+++ b/R/coerce-sp-methods.R
@@ -0,0 +1,93 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2012-2013, 2015 Sebastian Meyer
+### Time-stamp: <[coerce-sp-methods.R] 2015-02-25 22:43 (CET) by SM>
+################################################################################
+
+
+##' Coerce \code{"SpatialPolygons"} to \code{"owin"}
+##'
+##' Package \pkg{polyCub} implements \code{coerce}-methods
+##' (\code{as(object, Class)}) to convert \code{"\linkS4class{SpatialPolygons}"}
+##' (or \code{"\linkS4class{Polygons}"} or \code{"\linkS4class{Polygon}"})
+##' to \code{"\link[=owin.object]{owin}"}.
+##' They are also registered as \code{\link{as.owin}}-methods to support
+##' \code{\link{polyCub.midpoint}}.
+##' Note that the \pkg{maptools} package contains an alternative implementation
+##' of coercion from \code{"SpatialPolygons"} to \code{"owin"} (and reverse),
+##' and \R will use the S4 \code{coerce}-method that was loaded last,
+##' and prefer the \code{as.owin.SpatialPolygons} S3-method exported from
+##' \pkg{maptools} if attached.
+##' @author Sebastian Meyer
+##' @keywords spatial methods
+##' @name coerce-sp-methods
+##' @rdname coerce-sp-methods
+##' @exportMethod coerce
+NULL
+
+##' @param W an object of class \code{"SpatialPolygons"},
+##' \code{"Polygons"}, or \code{"Polygon"}.
+##' @param ... further arguments passed to \code{\link{owin}}.
+##' @rdname coerce-sp-methods
+##' @importFrom spatstat owin
+##' @method as.owin SpatialPolygons
+##' @export
+as.owin.SpatialPolygons <- function (W, ...)
+ owin(poly = xylist.SpatialPolygons(W), ...)
+
+##' @rdname coerce-sp-methods
+##' @importFrom spatstat owin
+##' @method as.owin Polygons
+##' @export
+as.owin.Polygons <- function (W, ...)
+ owin(poly = xylist.Polygons(W), ...)
+
+##' @rdname coerce-sp-methods
+##' @importFrom spatstat owin
+##' @method as.owin Polygon
+##' @export
+as.owin.Polygon <- function (W, ...)
+ owin(poly = xylist.Polygon(W), ...)
+
+
+## Register "owin" as class in S4 so we can define methods for it
+##setClass("owin")
+## -> no need to register "owin", since we depend on sp which does it !
+## Otherwise we would get the following warning upon package installation:
+## Warning in .simpleDuplicateClass(def, prev) :
+## the specification for class "owin" in package 'polyCub' seems
+## equivalent to one from package 'sp' and is not turning on
+## duplicate class definitions for this class
+## Using setOldClass("owin") is incompatible with package "maptools", which
+## does setClass("owin") _and_ exports this class! Specifically, loading
+## library("polyCub"); library("maptools"); library("gpclib")
+## in this order would not work (no idea why) throwing:
+## Error : package slot missing from signature for generic 'plot'
+## and classes gpc.poly, ANY
+## cannot use with duplicate class names (the package may need to be
+## re-installed)
+## Error: package/namespace load failed for 'gpclib'
+
+##' @name coerce,SpatialPolygons,owin-method
+##' @rdname coerce-sp-methods
+setAs(from = "SpatialPolygons", to = "owin",
+ def = function (from) as.owin.SpatialPolygons(from))
+
+##' @name coerce,Polygons,owin-method
+##' @rdname coerce-sp-methods
+setAs(from = "Polygons", to = "owin",
+ def = function (from) as.owin.Polygons(from))
+
+##' @name coerce,Polygon,owin-method
+##' @rdname coerce-sp-methods
+setAs(from = "Polygon", to = "owin",
+ def = function (from) as.owin.Polygon(from))
+
+
+##' @name coerce,Polygon,Polygons-method
+##' @rdname coerce-sp-methods
+setAs(from = "Polygon", to = "Polygons",
+ def = function (from) Polygons(list(from), "Polygon"))
diff --git a/R/plotpolyf.R b/R/plotpolyf.R
new file mode 100644
index 0000000..a725752
--- /dev/null
+++ b/R/plotpolyf.R
@@ -0,0 +1,86 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2013-2014 Sebastian Meyer
+### Time-stamp: <[plotpolyf.R] 2014-09-27 14:39 (CEST) by SM>
+###
+### Plot polygonal domain with image of bivariate function
+################################################################################
+
+
+##' Plot Polygonal Domain on Image of Bivariate Function
+##'
+##' Produces a combined plot of a polygonal domain and an image of a bivariate
+##' function, using either \code{\link[lattice:levelplot]{lattice::levelplot}}
+##' or \code{\link{image}}.
+##'
+##' @param polyregion a polygonal domain.
+##' The following classes are supported: \code{"\link[spatstat]{owin}"},
+##' \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"},
+##' \code{"\linkS4class{SpatialPolygons}"}, \code{"\linkS4class{Polygons}"},
+##' and \code{"\linkS4class{Polygon}"}
+##' (for these we have an internal \code{\link{xylist}} method).
+##' @param f a two-dimensional real function.
+##' As its first argument it must take a coordinate matrix, i.e., a
+##' numeric matrix with two columns, and it must return a numeric vector of
+##' length the number of coordinates.
+##' @param ... further arguments for \code{f}.
+##' @param npixel numeric vector of length 1 or 2 setting the number of pixels
+##' in each dimension.
+##' @param cuts number of cut points in the \eqn{z} dimension.
+##' The range of function values will be divided into \code{cuts+1} levels.
+##' @param col colour vector used for the function levels.
+##' @param lwd line width of the polygon edges.
+##' @param xlim,ylim numeric vectors of length 2 setting the axis limits.
+##' \code{NULL} means using the bounding box of \code{polyregion}.
+##' @param use.lattice logical indicating if \pkg{lattice} graphics
+##' (\code{\link[lattice]{levelplot}}) should be used.
+##' @param print.args a list of arguments passed to \code{\link{print.trellis}}
+##' for plotting the produced \code{\link[=trellis.object]{"trellis"}} object
+##' (given \code{use.lattice = TRUE}). The latter will be returned without
+##' explicit \code{print}ing if \code{print.args} is not a list.
+##' @author Sebastian Meyer
+##' @keywords hplot
+##' @example inst/examples/plotpolyf.R
+##' @importFrom grDevices extendrange heat.colors
+##' @importFrom graphics image
+##' @export
+
+plotpolyf <- function (polyregion, f, ...,
+ npixel=100, cuts=15, col=rev(heat.colors(cuts+1)), lwd=3,
+ xlim=NULL, ylim=NULL, use.lattice=TRUE, print.args=list())
+{
+ polys <- xylist(polyregion)
+ npixel <- rep(npixel, length.out=2)
+
+ ## make two-dimensional grid
+ if (is.null(xlim))
+ xlim <- extendrange(unlist(lapply(polys, "[[", "x"), use.names=FALSE))
+ if (is.null(ylim))
+ ylim <- extendrange(unlist(lapply(polys, "[[", "y"), use.names=FALSE))
+ xgrid <- makegrid(xlim, npixel[1])
+ ygrid <- makegrid(ylim, npixel[2])
+ xygrid <- expand.grid(x=xgrid, y=ygrid, KEEP.OUT.ATTRS=FALSE)
+
+ ## compute function values on the grid
+ xygrid$fval <- f(as.matrix(xygrid, rownames.force = FALSE), ...)
+
+ ## plot
+ if (use.lattice && requireNamespace("lattice")) {
+ mypanel <- function(...) {
+ lattice::panel.levelplot(...)
+ lapply(polys, function(xy) lattice::panel.polygon(xy, lwd=lwd))
+ }
+ trobj <- lattice::levelplot(fval ~ x*y, data=xygrid, aspect="iso",
+ cuts=cuts, col.regions=col, panel=mypanel)
+ if (is.list(print.args)) {
+ do.call("print", c(alist(x=trobj), print.args))
+ } else trobj
+ } else {
+ image(xgrid, ygrid, matrix(xygrid$fval, npixel[1], npixel[2]), col=col,
+ xlab="x", ylab="y", asp=1)
+ plot_polyregion(polyregion, lwd=lwd, add=TRUE)
+ }
+}
diff --git a/R/polyCub.R b/R/polyCub.R
new file mode 100644
index 0000000..966518c
--- /dev/null
+++ b/R/polyCub.R
@@ -0,0 +1,50 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2009-2013 Sebastian Meyer
+### Time-stamp: <[polyCub.R] by SM Sam 06/07/2013 12:52 (CEST)>
+################################################################################
+
+
+#' Wrapper Function for the Various Cubature Methods
+#'
+#' Instead of calling one of the specific cubature methods of this package, the
+#' wrapper function \code{polyCub} may be used together with the \code{method}
+#' argument.
+#'
+#' @param polyregion a polygonal integration domain.
+#' The supported classes depend on the specific method, however, the
+#' \code{"\link[spatstat]{owin}"} class from package \pkg{spatstat} works for
+#' all methods, as well should a \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"}
+#' polygon (but see the comments in \code{help("\link{coerce-methods}")}).
+#' @param f two-dimensional function to be integrated.
+#' As its first argument the function must take a coordinate matrix, i.e. a
+#' numeric matrix with two columns. For the \code{"exact.Gauss"} \code{method},
+#' \code{f} is ignored since it is specific to the bivariate normal density.
+#' @param method choose one of the implemented cubature methods (partial
+#' argument matching is applied), see \code{help("\link{polyCub-package}")}
+#' for an overview. Defaults to using the product Gauss cubature
+#' implemented in \code{\link{polyCub.SV}}.
+#' @param ... arguments of \code{f} or of the specific \code{method}.
+#' @param plot logical indicating if an illustrative plot of the numerical
+#' integration should be produced.
+#' @return The approximated integral of \code{f} over \code{polyregion}.
+#' @example inst/examples/polyCub.R
+#' @keywords math spatial
+#' @family polyCub-methods
+#' @export
+
+polyCub <- function (polyregion, f,
+ method = c("SV", "midpoint", "iso", "exact.Gauss"), ...,
+ plot = FALSE)
+{
+ method <- match.arg(method)
+ cl <- match.call()
+ cl$method <- NULL
+ cl[[1]] <- as.name(paste("polyCub", method, sep="."))
+ if (method == "exact.Gauss") cl$f <- NULL
+ int <- eval(cl, parent.frame())
+ int #structure(int, method = method)
+}
diff --git a/R/polyCub.SV.R b/R/polyCub.SV.R
new file mode 100644
index 0000000..a399d05
--- /dev/null
+++ b/R/polyCub.SV.R
@@ -0,0 +1,352 @@
+################################################################################
+### polyCub.SV: Product Gauss Cubature over Polygonal Domains
+###
+### Copyright (C) 2009-2014,2017 Sebastian Meyer
+###
+### This file is part of the R package "polyCub",
+### free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.R-project.org/Licenses/.
+################################################################################
+
+
+##' Product Gauss Cubature over Polygonal Domains
+##'
+##' Product Gauss cubature over polygons as proposed by
+##' Sommariva and Vianello (2007).
+##'
+##' @inheritParams plotpolyf
+##' @param f a two-dimensional real function (or \code{NULL} to only compute
+##' nodes and weights).
+##' As its first argument it must take a coordinate matrix, i.e., a
+##' numeric matrix with two columns, and it must return a numeric vector of
+##' length the number of coordinates.
+##' @param nGQ degree of the one-dimensional Gauss-Legendre quadrature rule
+##' (default: 20) as implemented in function \code{\link[statmod]{gauss.quad}}
+##' of package \pkg{statmod}. Nodes and weights up to \code{nGQ=60} are cached
+##' in \pkg{polyCub}, for larger degrees \pkg{statmod} is required.
+##' @param alpha base-line of the (rotated) polygon at \eqn{x = \alpha} (see
+##' Sommariva and Vianello (2007) for an explication). If \code{NULL} (default),
+##' the midpoint of the x-range of each polygon is chosen if no \code{rotation}
+##' is performed, and otherwise the \eqn{x}-coordinate of the rotated point
+##' \code{"P"} (see \code{rotation}). If \code{f} has its maximum value at the
+##' origin \eqn{(0,0)}, e.g., the bivariate Gaussian density with zero mean,
+##' \code{alpha = 0} is a reasonable choice.
+##' @param rotation logical (default: \code{FALSE}) or a list of points
+##' \code{"P"} and \code{"Q"} describing the preferred direction. If
+##' \code{TRUE}, the polygon is rotated according to the vertices \code{"P"} and
+##' \code{"Q"}, which are farthest apart (see Sommariva and Vianello, 2007). For
+##' convex polygons, this rotation guarantees that all nodes fall inside the
+##' polygon.
+##' @param engine character string specifying the implementation to use.
+##' Up to \pkg{polyCub} version 0.4-3, the two-dimensional nodes and weights
+##' were computed by \R functions and these are still available by setting
+##' \code{engine = "R"}.
+##' The new C-implementation is now the default (\code{engine = "C"}) and
+##' requires approximately 30\% less computation time.\cr
+##' The special setting \code{engine = "C+reduce"} will discard redundant nodes
+##' at (0,0) with zero weight resulting from edges on the base-line
+##' \eqn{x = \alpha} or orthogonal to it.
+##' This extra cleaning is only worth its cost for computationally intensive
+##' functions \code{f} over polygons which really have some edges on the
+##' baseline or parallel to the x-axis. Note that the old \R
+##' implementation does not have such unset zero nodes and weights.
+##' @param plot logical indicating if an illustrative plot of the numerical
+##' integration should be produced.
+##' @return The approximated value of the integral of \code{f} over
+##' \code{polyregion}.\cr
+##' In the case \code{f = NULL}, only the computed nodes and weights are
+##' returned in a list of length the number of polygons of \code{polyregion},
+##' where each component is a list with \code{nodes} (a numeric matrix with
+##' two columns), \code{weights} (a numeric vector of length
+##' \code{nrow(nodes)}), the rotation \code{angle}, and \code{alpha}.
+##' @author Sebastian Meyer\cr
+##' The product Gauss cubature is based on the
+##' original \acronym{MATLAB} implementation \code{polygauss} by Sommariva and
+##' Vianello (2007), which is available under the GNU GPL (>=2) license from
+##' \url{http://www.math.unipd.it/~alvise/software.html}.
+##' @references
+##' Sommariva, A. and Vianello, M. (2007).
+##' Product Gauss cubature over polygons based on Green's integration formula.
+##' \emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+##' @keywords math spatial
+##' @family polyCub-methods
+##' @importFrom graphics points
+##' @examples # see example(polyCub)
+##' @export
+
+polyCub.SV <- function (polyregion, f, ...,
+ nGQ = 20, alpha = NULL, rotation = FALSE, engine = "C",
+ plot = FALSE)
+{
+ polys <- xylist(polyregion) # transform to something like "owin$bdry"
+ # which means anticlockwise vertex order with
+ # first vertex not repeated
+ stopifnot(isScalar(nGQ), nGQ > 0,
+ is.null(alpha) || (isScalar(alpha) && !is.na(alpha)))
+
+ ## COMPUTE NODES AND WEIGHTS OF 1D GAUSS QUADRATURE RULE.
+ ## DEGREE "N" (as requested) (ORDER GAUSS PRIMITIVE)
+ nw_N <- gauss.quad(nGQ)
+ ## DEGREE "M" = N+1 (ORDER GAUSS INTEGRATION)
+ nw_M <- gauss.quad(nGQ + 1)
+
+ ## Special case f=NULL: compute and return nodes and weights only
+ if (is.null(f)) {
+ return(lapply(X = polys, FUN = polygauss, nw_MN = c(nw_M, nw_N),
+ alpha = alpha, rotation = rotation, engine = engine))
+ }
+
+ ## Cubature over every single polygon of the "polys" list
+ f <- match.fun(f)
+ int1 <- function (poly) {
+ nw <- polygauss(poly, c(nw_M, nw_N), alpha, rotation, engine)
+ fvals <- f(nw$nodes, ...)
+ cubature_val <- sum(nw$weights * fvals)
+ ## if (!isTRUE(all.equal(0, cubature_val))) {
+ ## if ((1 - 2 * as.numeric(poly$hole)) * sign(cubature_val) == -1)
+ ## warning("wrong sign if positive integral")
+ ## }
+ cubature_val
+ }
+ respolys <- vapply(X=polys, FUN=int1, FUN.VALUE=0, USE.NAMES=FALSE)
+ int <- sum(respolys)
+
+### ILLUSTRATION ###
+ if (plot) {
+ plotpolyf(polys, f, ..., use.lattice=FALSE)
+ for (i in seq_along(polys)) {
+ nw <- polygauss(polys[[i]], c(nw_M, nw_N), alpha, rotation, engine)
+ points(nw$nodes, cex=0.6, pch = i) #, col=1+(nw$weights<=0)
+ }
+ }
+###################
+
+ int
+}
+
+## this wrapper provides a partially memoized version of
+## unname(statmod::gauss.quad(n, kind="legendre"))
+gauss.quad <- function (n)
+{
+ if (n <= 61) { # results cached in R/sysdata.rda
+ .NWGL[[n]]
+ } else if (requireNamespace("statmod")) {
+ unname(statmod::gauss.quad(n = n, kind = "legendre"))
+ } else {
+ stop("package ", sQuote("statmod"), " is required for nGQ > 60")
+ }
+}
+
+
+##' Calculate 2D Nodes and Weights of the Product Gauss Cubature
+##'
+##' @param xy list with elements \code{"x"} and \code{"y"} containing the
+##' polygon vertices in \emph{anticlockwise} order (otherwise the result of the
+##' cubature will have a negative sign) with first vertex not repeated at the
+##' end (like \code{owin.object$bdry}).
+##' @param nw_MN unnamed list of nodes and weights of one-dimensional Gauss
+##' quadrature rules of degrees \eqn{N} and \eqn{M=N+1} (as returned by
+##' \code{\link[statmod]{gauss.quad}}): \code{list(s_M, w_M, s_N, w_N)}.
+##' @inheritParams polyCub.SV
+##' @references
+##' Sommariva, A. and Vianello, M. (2007):
+##' Product Gauss cubature over polygons based on Green's integration formula.
+##' \emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+##' @keywords internal
+##' @useDynLib polyCub, .registration = TRUE
+
+polygauss <- function (xy, nw_MN, alpha = NULL, rotation = FALSE, engine = "C")
+{
+ ## POLYGON ROTATION
+
+ xyrot <- if (identical(FALSE, rotation)) {
+ if (is.null(alpha)) { # choose midpoint of x-range
+ xrange <- range(xy[["x"]])
+ alpha <- (xrange[1L] + xrange[2L]) / 2
+ }
+ angle <- 0
+ xy[c("x", "y")]
+ } else {
+ ## convert to coordinate matrix
+ xy <- cbind(xy[["x"]], xy[["y"]], deparse.level=0)
+ ## determine P and Q
+ if (identical(TRUE, rotation)) { # automatic choice of rotation angle
+ ## such that for a convex polygon all nodes fall inside the polygon
+ QP <- vertexpairmaxdist(xy)
+ Q <- QP[1L,,drop=TRUE]
+ P <- QP[2L,,drop=TRUE]
+ } else if (is.list(rotation)) { # predefined rotation
+ P <- rotation$P
+ Q <- rotation$Q
+ stopifnot(is.vector(P, mode="numeric") && length(P) == 2L,
+ is.vector(Q, mode="numeric") && length(Q) == 2L)
+ stopifnot(any(P != Q))
+ rotation <- TRUE
+ } else {
+ stop("'rotation' must be logical or a list of points \"P\" and \"Q\"")
+ }
+ rotmat <- rotmatPQ(P,Q)
+ angle <- attr(rotmat, "angle")
+ if (is.null(alpha)) {
+ Prot <- rotmat %*% P
+ alpha <- Prot[1]
+ }
+ xyrot <- xy %*% t(rotmat) # = t(rotmat %*% t(xy))
+ ## convert back to list
+ list(x = xyrot[,1L,drop=TRUE], y = xyrot[,2L,drop=TRUE])
+ }
+
+ ## number of vertices
+ L <- length(xyrot[[1L]])
+
+
+ ## COMPUTE 2D NODES AND WEIGHTS.
+
+ if (engine == "R") {
+
+ toIdx <- c(seq.int(2, L), 1L)
+ nwlist <- mapply(.polygauss.side,
+ xyrot[[1L]], xyrot[[2L]],
+ xyrot[[1L]][toIdx], xyrot[[2L]][toIdx],
+ MoreArgs = c(nw_MN, alpha),
+ SIMPLIFY = FALSE, USE.NAMES = FALSE)
+
+ nodes <- c(lapply(nwlist, "[[", 1L),
+ lapply(nwlist, "[[", 2L),
+ recursive=TRUE)
+ dim(nodes) <- c(length(nodes)/2, 2L)
+ weights <- unlist(lapply(nwlist, "[[", 3L),
+ recursive=FALSE, use.names=FALSE)
+
+ } else { # use C-implementation
+
+ ## degrees of cubature and vector template for results
+ M <- length(nw_MN[[1L]])
+ N <- length(nw_MN[[3L]])
+ zerovec <- double(L*M*N)
+
+ ## rock'n'roll
+ nwlist <- .C(C_polygauss,
+ as.double(xyrot[[1L]]), as.double(xyrot[[2L]]),
+ as.double(nw_MN[[1L]]), as.double(nw_MN[[2L]]),
+ as.double(nw_MN[[3L]]), as.double(nw_MN[[4L]]),
+ as.double(alpha),
+ as.integer(L), as.integer(M), as.integer(N),
+ x = zerovec, y = zerovec, w = zerovec)[c("x", "y", "w")]
+
+ nodes <- cbind(nwlist[[1L]], nwlist[[2L]], deparse.level=0)
+ weights <- nwlist[[3L]]
+
+ ## remove unset nodes from edges on baseline or orthogonal to it
+ ## (note that the R implementation does not return such redundant nodes)
+ if (engine == "C+reduce" && any(unset <- weights == 0)) {
+ nodes <- nodes[!unset,]
+ weights <- weights[!unset]
+ }
+
+ }
+
+ ## back-transform rotated nodes by t(t(rotmat) %*% t(nodes))
+ ## (inverse of rotation matrix is its transpose)
+ list(nodes = if (rotation) nodes %*% rotmat else nodes,
+ weights = weights, angle = angle, alpha = alpha)
+}
+
+
+## The working horse .polygauss.side below is an R translation
+## of the original MATLAB implementation by Sommariva and Vianello (2007).
+
+.polygauss.side <- function (x1, y1, x2, y2, s_loc, w_loc, s_N, w_N, alpha)
+{
+ if ((x1 == alpha && x2 == alpha) || (y2 == y1))
+ ## side lies on base-line or is orthogonal to it -> skip
+ return(NULL)
+
+ if (x2 == x1) { # side is parallel to base-line => degree N
+ s_loc <- s_N
+ w_loc <- w_N
+ }
+
+ half_pt_x <- (x1+x2)/2
+ half_length_x <- (x2-x1)/2
+
+ half_pt_y <- (y1+y2)/2
+ half_length_y <- (y2-y1)/2
+
+ ## GAUSSIAN POINTS ON THE SIDE.
+ x_gauss_side <- half_pt_x + half_length_x * s_loc
+ y_gauss_side <- half_pt_y + half_length_y * s_loc
+
+ scaling_fact_minus <- (x_gauss_side - alpha) / 2
+
+ ## construct nodes and weights: x and y coordinates ARE STORED IN MATRICES.
+ ## A COUPLE WITH THE SAME INDEX IS A POINT, i.e. P_i=(x(k),y(k)).
+ ## Return in an unnamed list of nodes_x, nodes_y, weights
+ ## (there is no need for c(nodes_x) and c(weights))
+ list(alpha + tcrossprod(scaling_fact_minus, s_N + 1), # degree_loc x N
+ rep.int(y_gauss_side, length(s_N)), # length: degree_loc*N
+ tcrossprod(half_length_y*scaling_fact_minus*w_loc, w_N)) # degree_loc x N
+}
+
+## NOTE: The above .polygauss.side() function is already efficient R code.
+## Passing via C only at this deep level (see below) turned out to be
+## slower than staying with R! However, stepping into C already for
+## looping over the edges in polygauss() improves the speed.
+## ## @useDynLib polyCub C_polygauss_side
+## .polygauss.side <- function (x1, y1, x2, y2, s_M, w_M, s_N, w_N, alpha)
+## {
+## if ((x1 == alpha && x2 == alpha) || (y2 == y1))
+## ## side lies on base-line or is orthogonal to it -> skip
+## return(NULL)
+##
+## parallel2baseline <- x2 == x1 # side is parallel to base-line => degree N
+## M <- length(s_M)
+## N <- length(s_N)
+## loc <- if (parallel2baseline) N else M
+## zerovec <- double(loc * N)
+## .C(C_polygauss_side,
+## as.double(x1), as.double(y1), as.double(x2), as.double(y2),
+## as.double(if (parallel2baseline) s_N else s_M),
+## as.double(if (parallel2baseline) w_N else w_M),
+## as.double(s_N), as.double(w_N), as.double(alpha),
+## as.integer(loc), as.integer(N),
+## x = zerovec, y = zerovec, w = zerovec)[c("x", "y", "w")]
+## }
+
+
+##' @importFrom stats dist
+vertexpairmaxdist <- function (xy)
+{
+ ## compute euclidean distance matrix
+ distances <- dist(xy)
+ size <- attr(distances, "Size")
+
+ ## select two points with maximum distance
+ maxdistidx <- which.max(distances)
+ lowertri <- seq_along(distances) == maxdistidx
+ mat <- matrix(FALSE, size, size)
+ mat[lower.tri(mat)] <- lowertri
+ QPidx <- which(mat, arr.ind=TRUE, useNames=FALSE)[1L,]
+ xy[QPidx,]
+}
+
+rotmatPQ <- function (P, Q)
+{
+ direction_axis <- (Q-P) / sqrt(sum((Q-P)^2))
+
+ ## determine rotation angle
+ rot_angle_x <- acos(direction_axis[1L])
+ rot_angle_y <- acos(direction_axis[2L])
+
+ rot_angle <- if (rot_angle_y <= pi/2) {
+ if (rot_angle_x <= pi/2) -rot_angle_y else rot_angle_y
+ } else {
+ if (rot_angle_x <= pi/2) pi-rot_angle_y else rot_angle_y
+ }
+ ## cat(sprintf(' [ANGLE CLOCKWISE (IN DEGREES)]: %5.5f\n', rot_angle*180/pi))
+
+ ## rotation matrix
+ rot_matrix <- diag(cos(rot_angle), nrow=2L)
+ rot_matrix[2:3] <- c(-1,1) * sin(rot_angle) # clockwise rotation
+ structure(rot_matrix, angle=rot_angle)
+}
diff --git a/R/polyCub.exact.Gauss.R b/R/polyCub.exact.Gauss.R
new file mode 100644
index 0000000..62b7fee
--- /dev/null
+++ b/R/polyCub.exact.Gauss.R
@@ -0,0 +1,217 @@
+################################################################################
+### polyCub.exact.Gauss: Quasi-Exact Cubature of the Bivariate Normal Density
+###
+### Copyright (C) 2009-2016 Sebastian Meyer
+###
+### This file is part of the R package "polyCub",
+### free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+################################################################################
+
+
+#' Quasi-Exact Cubature of the Bivariate Normal Density
+#'
+#' Integration is based on triangulation of the (transformed) polygonal domain
+#' and formulae from the
+#' Abramowitz and Stegun (1972) handbook (Section 26.9, Example 9, pp. 956f.).
+#' This method is quite cumbersome because the A&S formula is only for triangles
+#' where one vertex is the origin (0,0). For each triangle of the
+#' \code{\link[gpclib]{tristrip}} we have to check in which of the 6 outer
+#' regions of the triangle the origin (0,0) lies and adapt the signs in the
+#' formula appropriately: \eqn{(AOB+BOC-AOC)} or \eqn{(AOB-AOC-BOC)} or
+#' \eqn{(AOB+AOC-BOC)} or \eqn{(AOC+BOC-AOB)} or \ldots.
+#' However, the most time consuming step is the
+#' evaluation of \code{\link[mvtnorm]{pmvnorm}}.
+#'
+#' @note The package \pkg{gpclib} (which is required to produce the
+#' \code{tristrip}, since this is not yet implemented in \pkg{rgeos})
+#' has a restricted license (commercial use prohibited).
+#' It has to be accepted explicitly via
+#' \code{\link{gpclibPermit}()} prior to using \code{polyCub.exact.Gauss}.
+#'
+#' @param polyregion a \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} polygon or
+#' something that can be coerced to this class, e.g., an \code{"owin"} polygon
+#' (converted via \code{\link{owin2gpc}} and -- given \pkg{rgeos} is available
+#' -- \code{"SpatialPolygons"} also work.
+#' @param mean,Sigma mean and covariance matrix of the bivariate normal density
+#' to be integrated.
+#' @param plot logical indicating if an illustrative plot of the numerical
+#' integration should be produced. Note that the \code{polyregion} will be
+#' transformed (shifted and scaled).
+#' @return The integral of the bivariate normal density over \code{polyregion}.
+#' Two attributes are appended to the integral value:
+#' \item{nEval}{
+#' number of triangles over which the standard bivariate normal density had to
+#' be integrated, i.e. number of calls to \code{\link[mvtnorm]{pmvnorm}} and
+#' \code{\link[stats]{pnorm}}, the former of which being the most time-consuming
+#' operation.
+#' }
+#' \item{error}{
+#' Approximate absolute integration error stemming from the error introduced by
+#' the \code{nEval} \code{\link[mvtnorm]{pmvnorm}} evaluations.
+#' For this reason, the cubature method is in fact only
+#' quasi-exact (as is the \code{pmvnorm} function).
+#' }
+#' @references
+#' Abramowitz, M. and Stegun, I. A. (1972).
+#' Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+#' Tables. New York: Dover Publications.
+#' @keywords math spatial
+#' @seealso \code{\link{circleCub.Gauss}} for quasi-exact cubature of the
+#' isotropic Gaussian density over a circular domain.
+#' @family polyCub-methods
+#' @examples # see example(polyCub)
+#' @import methods
+#' @import sp
+#' @importFrom stats cov2cor
+#' @importFrom spatstat is.polygonal
+#' @importFrom graphics lines
+#' @export
+## NOTE: we don't import graphics::plot since it is already imported via sp
+
+polyCub.exact.Gauss <- function (polyregion, mean = c(0,0), Sigma = diag(2),
+ plot = FALSE)
+{
+ gpclibCheck(fatal=TRUE)
+ if (is.polygonal(polyregion)) {
+ polyregion <- owin2gpc(polyregion)
+ } else if (!inherits(polyregion, "gpc.poly")) {
+ if (inherits(polyregion, "SpatialPolygons") &&
+ !requireNamespace("rgeos")) {
+ stop("package ", sQuote("rgeos"), " is required to handle ",
+ "\"SpatialPolygons\" input")
+ }
+ polyregion <- as(polyregion, "gpc.poly")
+ }
+
+ ## coordinate transformation so that the standard bivariat normal density
+ ## can be used in integrations (cf. formula 26.3.22)
+ polyregion at pts <- transform_pts(polyregion at pts, mean = mean, Sigma = Sigma)
+
+ ## triangulation: tristrip() returns a list where each element is a
+ ## coordinate matrix of vertices of triangles
+ triangleSets <- gpclib::tristrip(polyregion)
+
+### ILLUSTRATION ###
+ if (plot) {
+ plot(polyregion, poly.args=list(lwd=2), ann=FALSE)
+ lapply(triangleSets, lines, lty=2)
+ }
+####################
+
+ integrals <- vapply(X = triangleSets, FUN = function (triangles) {
+ int <- 0
+ error <- 0
+ nTriangles <- nrow(triangles) - 2L
+ for (i in seq_len(nTriangles)) {
+ res <- .intTriangleAS(triangles[i+(0:2),])
+ int <- int + res
+ error <- error + attr(res, "error")
+ }
+ c(int, nTriangles, error)
+ }, FUN.VALUE = numeric(3L), USE.NAMES = FALSE)
+ int <- sum(integrals[1,])
+
+ ## number of .V() evaluations (if there were no degenerate triangles)
+ attr(int, "nEval") <- 6 * sum(integrals[2,])
+ ## approximate absolute integration error
+ attr(int, "error") <- sum(integrals[3,])
+ return(int)
+}
+
+
+
+###########################
+### Auxiliary Functions ###
+###########################
+
+## transform coordinates according to Formula 26.3.22
+transform_pts <- function (pts, mean, Sigma)
+{
+ mx <- mean[1L]
+ my <- mean[2L]
+ rho <- cov2cor(Sigma)[1L,2L]
+ sdx <- sqrt(Sigma[1L,1L])
+ sdy <- sqrt(Sigma[2L,2L])
+ lapply(pts, function (poly) {
+ x0 <- (poly[["x"]] - mx) / sdx
+ y0 <- (poly[["y"]] - my) / sdy
+ list(x = (x0 + y0) / sqrt(2 + 2*rho),
+ y = (y0 - x0) / sqrt(2 - 2*rho),
+ hole = poly[["hole"]])
+ })
+}
+
+## calculates the integral of the standard bivariat normal over a triangle ABC
+.intTriangleAS <- function (xy)
+{
+ if (anyDuplicated(xy)) # degenerate triangle
+ return(structure(0, error = 0))
+ A <- xy[1,]
+ B <- xy[2,]
+ C <- xy[3,]
+ intAOB <- .intTriangleAS0(A, B)
+ intBOC <- .intTriangleAS0(B, C)
+ intAOC <- .intTriangleAS0(A, C)
+
+ # determine signs of integrals
+ signAOB <- -1 + 2*.pointsOnSameSide(A,B,C)
+ signBOC <- -1 + 2*.pointsOnSameSide(B,C,A)
+ signAOC <- -1 + 2*.pointsOnSameSide(A,C,B)
+
+ int <- signAOB*intAOB + signBOC*intBOC + signAOC*intAOC
+ attr(int, "error") <- attr(intAOB, "error") +
+ attr(intBOC, "error") + attr(intAOC, "error")
+ return(int)
+}
+
+## calculates the integral of the standard bivariat normal over a triangle A0B
+.intTriangleAS0 <- function (A, B)
+{
+ BmA <- B - A
+ d <- sqrt(sum(BmA^2))
+ h <- abs(B[2L]*A[1L] - A[2L]*B[1L]) / d # distance of AB to the origin
+ if (d == 0 || h == 0) # degenerate triangle: A == B or 0, A, B on a line
+ return(structure(0, error = 0))
+
+ k1 <- dotprod(A, BmA) / d
+ k2 <- dotprod(B, BmA) / d
+ V2 <- .V(h, abs(k2))
+ V1 <- .V(h, abs(k1))
+ res <- if (sign(k1) == sign(k2)) {
+ ## A and B are on the same side of the normal line through 0
+ abs(V2 - V1)
+ } else {
+ V2 + V1
+ }
+ attr(res, "error") <- attr(V1, "error") + attr(V2, "error")
+ return(res)
+}
+
+## checks if point1 and point2 lie on the same side of a line through
+## linepoint1 and linepoint2
+.pointsOnSameSide <- function (linepoint1, linepoint2, point1, point2 = c(0,0))
+{
+ n <- c(-1,1) * rev(linepoint2-linepoint1) # normal vector
+ S <- dotprod(point1-linepoint1,n) * dotprod(point2-linepoint1,n)
+ return(S > 0)
+}
+
+## calculates the integral of the standard bivariat normal
+## over a triangle bounded by y=0, y=ax, x=h (cf. formula 26.3.23)
+##' @importFrom stats pnorm
+.V <- function(h,k) {
+ if (k == 0) # degenerate triangle
+ return(structure(0, error = 0))
+ a <- k/h
+ rho <- -a/sqrt(1+a^2)
+ # V = 0.25 + L(h,0,rho) - L(0,0,rho) - Q(h) / 2
+ # L(0,0,rho) = 0.25 + asin(rho) / (2*pi)
+ # V = L(h,0,rho) - asin(rho)/(2*pi) - Q(h) / 2
+ Lh0rho <- mvtnorm::pmvnorm(
+ lower = c(h,0), upper = c(Inf,Inf),
+ mean = c(0,0), corr = matrix(c(1,rho,rho,1), 2L, 2L)
+ )
+ Qh <- pnorm(h, mean = 0, sd = 1, lower.tail = FALSE)
+ return(Lh0rho - asin(rho)/2/pi - Qh/2)
+}
diff --git a/R/polyCub.iso.R b/R/polyCub.iso.R
new file mode 100644
index 0000000..88bc64d
--- /dev/null
+++ b/R/polyCub.iso.R
@@ -0,0 +1,252 @@
+################################################################################
+### polyCub.iso: Cubature of Isotropic Functions over Polygonal Domains
+###
+### Copyright (C) 2013-2017 Sebastian Meyer
+###
+### This file is part of the R package "polyCub",
+### free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+################################################################################
+
+
+#' Cubature of Isotropic Functions over Polygonal Domains
+#'
+#' Conducts numerical integration of a two-dimensional isotropic function
+#' \eqn{f(x,y) = f_r(||(x,y)-\boldsymbol{\mu}||)}{f(x,y) = f_r(||(x,y)-\mu||)},
+#' with \eqn{\mu} being the center of isotropy, over a polygonal domain.
+#' It internally solves a line integral along the polygon boundary using
+#' \code{\link{integrate}} where the integrand requires the antiderivative of
+#' \eqn{r f_r(r)}), which ideally is analytically available and supplied to the
+#' function as argument \code{intrfr}.
+#' The two-dimensional integration problem thereby reduces to an efficient
+#' adaptive quadrature in one dimension.
+#' See Meyer and Held (2014, Supplement B, Section 2.4) for mathematical
+#' details.
+#'
+#' @inheritParams plotpolyf
+#' @param intrfr analytical antiderivative of \eqn{r f_r(r)} from 0 to \code{R}
+#' (first argument, not necessarily named \code{"R"}, must be vectorized).
+#' If missing, \code{intrfr} is approximated numerically using
+#' \code{\link{integrate}} configured with \code{control}.
+#' @param ... further arguments for \code{f} or \code{intrfr}.
+#' @param center numeric vector of length 2, the center of isotropy.
+#' @param control list of arguments passed to \code{\link{integrate}}, the
+#' quadrature rule used for the line integral along the polygon boundary.
+#' @param check.intrfr logical (or numeric vector) indicating if
+#' (for which \code{r}'s) the supplied \code{intrfr} function should be
+#' checked against a numeric approximation. This check requires \code{f}
+#' to be specified. If \code{TRUE}, the set of test
+#' \code{r}'s defaults to a \code{\link{seq}} of length 20 from 1 to
+#' the maximum absolute x or y coordinate of any edge of the \code{polyregion}.
+#' @param plot logical indicating if an image of the function should be plotted
+#' together with the polygonal domain, i.e.,
+#' \code{\link{plotpolyf}(polyregion, f, \dots)}.
+#' @return The approximate integral of the isotropic function
+#' \code{f} over \code{polyregion}.\cr
+#' If the \code{intrfr} function is provided (which is assumed to be exact), an
+#' upper bound for the absolute integration error is appended to the result as
+#' attribute \code{"abs.error"}. It equals the sum of the absolute errors
+#' reported by all \code{\link{integrate}} calls
+#' (there is one for each edge of \code{polyregion}).
+#' @author Sebastian Meyer
+#'
+#' The basic mathematical formulation of this efficient integration for radially
+#' symmetric functions was ascertained with great support by
+#' Emil Hedevang (2013), Dept. of Mathematics, Aarhus University, Denmark.
+#' @references
+#' Hedevang, E. (2013). Personal communication at the Summer School on Topics in
+#' Space-Time Modeling and Inference (May 2013, Aalborg, Denmark).
+#'
+#' Meyer, S. and Held, L. (2014).
+#' Power-law models for infectious disease spread.
+#' \emph{The Annals of Applied Statistics}, \bold{8} (3), 1612-1639.\cr
+#' DOI-Link: \url{http://dx.doi.org/10.1214/14-AOAS743},
+#' \href{http://arxiv.org/abs/1308.5115}{arXiv:1308.5115}
+#' @seealso
+#' \code{system.file("include", "polyCubAPI.h", package = "polyCub")}
+#' for a full C-implementation of this cubature method (for a \emph{single}
+#' polygon). The corresponding C-routine \code{polyCub_iso} can be used by
+#' other \R packages, notably \pkg{surveillance}, via \samp{LinkingTo: polyCub}
+#' (in the \file{DESCRIPTION}) and \samp{#include <polyCubAPI.h>} (in suitable
+#' \file{/src} files). Note that the \code{intrfr} function must then also be
+#' supplied as a C-routine. An example can be found in the package tests.
+#' @keywords math spatial
+#' @family polyCub-methods
+#' @example inst/examples/polyCub.iso.R
+#' @importFrom stats integrate
+#' @export
+
+polyCub.iso <- function (polyregion, f, intrfr, ..., center,
+ control = list(), check.intrfr = FALSE, plot = FALSE)
+{
+ polys <- xylist(polyregion) # transform to something like "owin$bdry"
+ # which means anticlockwise vertex order with
+ # first vertex not repeated
+ getError <- !missing(intrfr) # can't estimate error of double approximation
+ center <- as.vector(center, mode = "numeric")
+ stopifnot(length(center) == 2L, is.finite(center))
+
+ ## check 'intrfr' function
+ rs <- if (isTRUE(check.intrfr)) {
+ seq(1, max(abs(unlist(lapply(polys, "[", c("x","y"))))), length.out=20L)
+ } else if (identical(check.intrfr, FALSE)) {
+ numeric(0L)
+ } else {
+ check.intrfr
+ }
+ intrfr <- checkintrfr(intrfr, f, ..., center=center, control=control, rs=rs)
+
+ ## plot polygon and function image
+ if (plot) plotpolyf(polys, f, ...)
+
+ ## do the cubature over all polygons of the 'polys' list
+ .polyCub.iso(polys, intrfr, ..., center=center,
+ control=control, .witherror=getError)
+}
+
+
+##' Check the Integral of \eqn{r f_r(r)}
+##'
+##' This function is auxiliary to \code{\link{polyCub.iso}}.
+##' The (analytical) integral of \eqn{r f_r(r)} from 0 to \eqn{R} is checked
+##' against a numeric approximation using \code{\link{integrate}} for various
+##' values of the upper bound \eqn{R}. A warning is issued if inconsistencies
+##' are found.
+##'
+##' @inheritParams polyCub.iso
+##' @param rs numeric vector of upper bounds for which to check the validity of
+##' \code{intrfr}. If it has length 0, no checks are performed.
+##' @param tolerance of \code{\link{all.equal.numeric}} when comparing
+##' \code{intrfr} results with numerical integration. Defaults to the
+##' relative tolerance used for \code{integrate}.
+##' @return The \code{intrfr} function. If it was not supplied, its quadrature
+##' version using \code{integrate} is returned.
+##' @importFrom stats integrate
+##' @export
+checkintrfr <- function (intrfr, f, ..., center, control = list(),
+ rs = numeric(0L), tolerance = control$rel.tol)
+{
+ doCheck <- length(rs) > 0L
+ if (!missing(f)) {
+ f <- match.fun(f)
+ rfr <- function (r, ...)
+ r * f(cbind(center[1L]+r, center[2L], deparse.level=0L), ...)
+ quadrfr1 <- function (R, ...) integrate(rfr, 0, R, ...)$value
+ if (length(control))
+ body(quadrfr1)[[2L]] <- as.call(c(as.list(body(quadrfr1)[[2L]]),
+ control))
+ quadrfr <- function (R, ...)
+ vapply(X = R, FUN = quadrfr1, FUN.VALUE = 0, ..., USE.NAMES = FALSE)
+ if (missing(intrfr)) {
+ return(quadrfr)
+ } else if (doCheck) {
+ cat("Checking 'intrfr' against a numeric approximation ... ")
+ stopifnot(is.vector(rs, mode="numeric"))
+ if (is.null(tolerance))
+ tolerance <- eval(formals(integrate)$rel.tol)
+ ana <- intrfr(rs, ...)
+ num <- quadrfr(rs, ...)
+ if (!isTRUE(comp <- all.equal(num, ana, tolerance=tolerance))) {
+ cat("\n->", comp, "\n")
+ warning("'intrfr' might be incorrect: ", comp)
+ } else cat("OK\n")
+ }
+ } else if (doCheck) {
+ stop("numerical verification of 'intrfr' requires 'f'")
+ }
+
+ match.fun(intrfr)
+}
+
+
+##' \code{.polyCub.iso} is a \dQuote{bare-bone} version of \code{polyCub.iso}.
+##' @rdname polyCub.iso
+##' @param polys something like \code{owin$bdry}, but see \code{\link{xylist}}.
+##' @param .witherror logical indicating if an upper bound for the absolute
+##' integration error should be attached as an attribute to the result?
+##' @export
+.polyCub.iso <- function (polys, intrfr, ..., center,
+ control = list(), .witherror = FALSE)
+{
+ ints <- lapply(polys, polyCub1.iso,
+ intrfr, ..., center=center,
+ control=control, .witherror=.witherror)
+ if (.witherror) {
+ res <- sum(vapply(X=ints, FUN="[", FUN.VALUE=0, 1L, USE.NAMES=FALSE))
+ attr(res, "abs.error") <-
+ sum(vapply(X=ints, FUN="[", FUN.VALUE=0, 2L, USE.NAMES=FALSE))
+ res
+ } else {
+ sum(unlist(ints, recursive=FALSE, use.names=FALSE))
+ }
+}
+
+## cubature method for a single polygon
+polyCub1.iso <- function (poly, intrfr, ..., center,
+ control = list(), .witherror = TRUE)
+{
+ xy <- cbind(poly[["x"]], poly[["y"]], deparse.level=0L)
+ nedges <- nrow(xy)
+ intedges <- erredges <- numeric(nedges)
+ for (i in seq_len(nedges)) {
+ v0 <- xy[i, ] - center
+ v1 <- xy[if (i==nedges) 1L else i+1L, ] - center
+ int <- lineInt(v0, v1, intrfr, ..., control=control)
+ intedges[i] <- int$value
+ erredges[i] <- int$abs.error
+ }
+ int <- sum(intedges)
+ ## if (!is.null(poly$hole) && !isTRUE(all.equal(0, int))) {
+ ## if ((1 - 2 * as.numeric(poly$hole)) * sign(int) == -1)
+ ## warning("wrong sign if positive integral")
+ ## }
+ if (.witherror) {
+ c(int, sum(erredges))
+ } else {
+ int
+ }
+}
+
+## line integral for one edge
+##' @importFrom stats integrate
+lineInt <- function (v0, v1, intrfr, ..., control = list())
+{
+ d <- v1 - v0
+ num <- v1[2L]*v0[1L] - v1[1L]*v0[2L] # = d[2]*p[,1] - d[1]*p[,2]
+ # for any point p on the edge
+ if (num == 0) { # i.e., if 'center' is part of this polygon edge
+ return(list(value = 0, abs.error = 0))
+ }
+ integrand <- function (t) {
+ ## get the points on the edge corresponding to t
+ p <- cbind(v0[1L] + t*d[1L], v0[2L] + t*d[2L], deparse.level=0L)
+ norm2 <- .rowSums(p^2, length(t), 2L)
+ ints <- intrfr(sqrt(norm2), ...)
+ ##ints[is.infinite(ints)] <- 1e300
+ num * ints / norm2
+ }
+ if (length(control)) { # use slower do.call()-construct
+ do.call("integrate", c(list(integrand, 0, 1), control))
+ } else {
+ integrate(integrand, 0, 1)
+ }
+}
+
+## equally fast method _only_ for convex polygonal domains including the origin
+## (formula obtained via polar coordinate representation)
+lineInt2 <- function (v0, v1, intrfr, ..., control = list())
+{
+ d <- v1 - v0
+ ld <- vecnorm(d)
+ l0 <- vecnorm(v0)
+ l1 <- vecnorm(v1)
+ dp <- dotprod(v0,v1)
+ theta <- acos((l0 - dp/l0) / ld)
+ num <- sin(theta) * l0
+ phispan <- acos(dp / l0 / l1)
+ integrand <- function (phi, ...) {
+ r <- num / sin(theta+phi)
+ intrfr(r, ...)
+ }
+ do.call("integrate", c(list(integrand, 0, phispan, ...), control))
+}
diff --git a/R/polyCub.midpoint.R b/R/polyCub.midpoint.R
new file mode 100644
index 0000000..0d207ad
--- /dev/null
+++ b/R/polyCub.midpoint.R
@@ -0,0 +1,80 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2009-2015 Sebastian Meyer
+### Time-stamp: <[polyCub.midpoint.R] 2015-02-25 20:53 (CET) by SM>
+################################################################################
+
+
+#' Two-Dimensional Midpoint Rule
+#'
+#' The surface is converted to a binary pixel image
+#' using the \code{\link[spatstat]{as.im.function}} method from package
+#' \pkg{spatstat} (Baddeley and Turner, 2005).
+#' The integral under the surface is then approximated as the
+#' sum over (pixel area * f(pixel midpoint)).
+#'
+#' @inheritParams plotpolyf
+#' @param polyregion a polygonal integration domain.
+#' It can be any object coercible to the \pkg{spatstat} class
+#' \code{"\link[spatstat]{owin}"} via a corresponding
+#' \code{\link[spatstat]{as.owin}}-method.
+#' Note that this includes polygons of the classes \code{"gpc.poly"} and
+#' \code{"\linkS4class{SpatialPolygons}"}, because \pkg{polyCub} defines
+#' methods \code{\link{as.owin.gpc.poly}} and
+#' \code{\link{as.owin.SpatialPolygons}}, respectively.
+#' @param eps width and height of the pixels (squares),
+#' see \code{\link[spatstat]{as.mask}}.
+#' @param dimyx number of subdivisions in each dimension,
+#' see \code{\link[spatstat]{as.mask}}.
+#' @param plot logical indicating if an illustrative plot of the numerical
+#' integration should be produced.
+#' @return The approximated value of the integral of \code{f} over
+#' \code{polyregion}.
+#' @references
+#' Baddeley, A. and Turner, R. (2005).
+#' \pkg{spatstat}: an \R package for analyzing spatial point patterns.
+#' \emph{Journal of Statistical Software}, \bold{12} (6), 1-42.
+#' @keywords math spatial
+#' @family polyCub-methods
+#' @import sp
+#' @importFrom spatstat as.im.function plot.im integral.im
+#' @importFrom grDevices gray
+#' @examples # see example(polyCub)
+#' @export
+## NOTE: we don't import graphics::plot since it is already imported via sp
+
+polyCub.midpoint <- function (polyregion, f, ...,
+ eps = NULL, dimyx = NULL, plot = FALSE)
+{
+ ## as.im needs seperate x and y arguments
+ fxy <- function (x, y, ...) f(cbind(x,y), ...)
+
+ ## calculate pixel values of fxy
+ IM <- tryCatch(
+ as.im.function(X=fxy, W=polyregion, ..., eps=eps, dimyx=dimyx),
+ error = function (e) {
+ ## if eps was to small such that the dimensions of the image would
+ ## be too big then the operation matrix(TRUE, nr, nc) throws an
+ ## error. (try e.g. devnull <- matrix(TRUE, 1e6,1e6))
+ ## unfortunately, it is not clear what we should do in this
+ ## case... => stop
+ stop("inapplicable choice of bandwidth (eps=", format(eps),
+ ") in midpoint rule:\n", e)
+ })
+
+### ILLUSTRATION ###
+ if (plot) {
+ plot.im(IM, axes=TRUE, col=gray(31:4/35), main="")
+ ## add evaluation points
+ #with(IM, points(expand.grid(xcol, yrow), col=!is.na(v), cex=0.5))
+ plot(polyregion, add=TRUE, poly.args=list(lwd=2), lwd=2)
+ ##<- two 'lwd'-specifications such that it works with owin and gpc.poly
+ }
+####################
+
+ ## return the approximated integral
+ integral.im(IM)
+}
diff --git a/R/sysdata.rda b/R/sysdata.rda
new file mode 100644
index 0000000..64ddd89
Binary files /dev/null and b/R/sysdata.rda differ
diff --git a/R/tools.R b/R/tools.R
new file mode 100644
index 0000000..0667d04
--- /dev/null
+++ b/R/tools.R
@@ -0,0 +1,98 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2009-2015 Sebastian Meyer
+### Time-stamp: <[tools.R] 2015-01-05 23:20 (CET) by SM>
+###
+### Tiny toolbox of internal function
+################################################################################
+
+
+##' Check if Polygon is Closed
+##'
+##' Check if the first and last coordinates of a coordinate matrix are
+##' identical.
+##' @param coords numeric coordinate matrix. It is interpreted by
+##' \code{\link{xy.coords}}.
+##' @return logical
+##' @keywords spatial internal
+##' @importFrom grDevices xy.coords
+isClosed <- function (coords)
+{
+ xycoords <- xy.coords(coords)[c("x","y")]
+ n <- length(xycoords$x)
+ return(identical(xycoords$x[1], xycoords$x[n]) &&
+ identical(xycoords$y[1], xycoords$y[n]))
+}
+
+
+##' Dot/Scalar Product of Two Vectors
+##'
+##' This is nothing else than \code{sum(x*y)}.
+##' @param x,y numeric vectors (of compatible lengths).
+##' @return \code{sum(x*y)}
+##' @keywords math internal
+dotprod <- function (x,y) sum(x*y)
+
+##' Euclidean Vector Norm (Length)
+##'
+##' This is nothing else than \code{sqrt(sum(x^2))}.
+##' @param x numeric vector.
+##' @return \code{sqrt(sum(x^2))}
+##' @keywords math internal
+vecnorm <- function (x) sqrt(sum(x^2))
+
+##' Checks if Argument is Scalar
+##'
+##' Check if the argument is scalar, i.e. a numeric vector of length 1.
+##' @param x any object
+##' @return logical
+##' @keywords internal
+isScalar <- function (x) {
+ length(x) == 1L && is.vector(x, mode = "numeric")
+}
+
+
+##' Plots a Polygonal Domain (of Various Classes)
+##'
+##' @inheritParams plotpolyf
+##' @param add logical. Add to existing plot?
+##' @import methods
+##' @import sp
+##' @importFrom graphics polygon
+##' @importFrom spatstat plot.owin
+## CAVE: need to import plot.owin for compatibility with spatstat <1.33-0,
+## since plot.owin was not registered as an S3-method for plot
+## NOTE: we don't import graphics::plot since it is already imported via sp
+plot_polyregion <- function (polyregion, lwd=2, add=FALSE)
+{
+ if (is.vector(polyregion, mode="list")) { # internal xylist object
+ stopifnot(add)
+ lapply(polyregion, polygon, lwd=lwd)
+ invisible()
+ } else if (inherits(polyregion, "gpc.poly")) {
+ plot(polyregion, poly.args=list(lwd=lwd), ann=FALSE, add=add)
+ } else {
+ if (inherits(polyregion, "Polygon"))
+ polyregion <- Polygons(list(polyregion), "ID")
+ if (inherits(polyregion, "Polygons"))
+ polyregion <- SpatialPolygons(list(polyregion))
+ ## plot call which works for "SpatialPolygons" and "owin"
+ plot(polyregion, lwd=lwd, axes=TRUE, main="", add=add)
+ }
+}
+
+
+##' Constructs Equally-Spaced Grid
+##'
+##' Construct an equally-spaced grid given a range and the number of cut points
+##' (one more than the number of resulting bins).
+##' This is nothing else than \code{seq(range[1], range[2], length.out=n)}.
+##' @param range numeric vector of length 2.
+##' @param n length of the desired grid, i.e. number of bins + 1.
+##' @return the desired grid, a numeric vector of length \code{n} covering
+##' \code{range}.
+##' @keywords internal
+makegrid <- function(range, n) seq(range[1], range[2], length.out=n)
diff --git a/R/xylist.R b/R/xylist.R
new file mode 100644
index 0000000..9598373
--- /dev/null
+++ b/R/xylist.R
@@ -0,0 +1,124 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2012-2014 Sebastian Meyer
+### Time-stamp: <[xylist.R] 2015-02-25 21:18 (CET) by SM>
+################################################################################
+
+
+##' Convert Various Polygon Classes to a Simple List of Vertices
+##'
+##' Different packages concerned with spatial data use different polygon
+##' specifications, which sometimes becomes very confusing (see Details below).
+##' To be compatible with the various polygon classes, package \pkg{polyCub}
+##' uses an S3 class \code{"xylist"}, which represents
+##' polygons by their core feature only, a list of lists of vertex coordinates
+##' (see the "Value" section below).
+##' The generic function \code{xylist} can deal with the
+##' following polygon classes:
+##' \itemize{
+##' \item{\code{"\link[=owin.object]{owin}"} from package \pkg{spatstat}}
+##' \item{\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} from package
+##' \pkg{rgeos} (or \pkg{gpclib})}
+##' \item{\code{"\linkS4class{Polygons}"} from package \pkg{sp}
+##' (as well as \code{"\linkS4class{Polygon}"} and
+##' \code{"\linkS4class{SpatialPolygons}"})}
+##' }
+##' The (somehow useless) default \code{xylist}-method
+##' does not perform any transformation but only ensures that the polygons are
+##' not closed (first vertex not repeated).
+##'
+##' Different packages concerned with spatial data use different polygon
+##' specifications with respect to:
+##' \itemize{
+##' \item{do we repeat the first vertex?}
+##' \item{which direction represents holes?}
+##' }
+##' Package overview:
+##' \describe{
+##' \item{\pkg{sp}:}{\emph{Repeat} first vertex at the end (closed),
+##' anticlockwise = hole, clockwise = normal boundary}
+##' \item{\pkg{spatstat}:}{do \emph{not repeat} first vertex,
+##' anticlockwise = normal boundary, clockwise = hole. This convention is also
+##' used in \code{xylist}.}
+##' \item{\pkg{gpclib}:}{Unfortunately, there seems to be no convention
+##' for the specification of polygons of class \code{"gpc.poly"}.}
+##' }
+##'
+##' @param object an object of one of the supported spatial classes.
+##' @param ... (unused) argument of the generic.
+##' @return Applying \code{xylist} to a polygon object, one gets a simple list,
+##' where each component (polygon) is a list of \code{"x"} and \code{"y"}
+##' coordinates. These represent vertex coordinates following \pkg{spatstat}'s
+##' \code{"owin"} convention (anticlockwise order without repeating any vertex).
+##' The opposite vertex order can be retained for the \pkg{sp}-classes
+##' by the non-default use with \code{reverse=FALSE}.
+##' @author Sebastian Meyer
+##' @name xylist
+##' @keywords spatial methods
+##' @export
+xylist <- function (object, ...) UseMethod("xylist")
+
+##' @rdname xylist
+##' @importFrom spatstat as.polygonal
+##' @export
+xylist.owin <- function (object, ...)
+{
+ as.polygonal(object)$bdry
+}
+
+##' @rdname xylist
+##' @export
+xylist.gpc.poly <- function (object, ...)
+{
+ xylist.owin(gpc2owin(object, check = FALSE))
+}
+
+##' @rdname xylist
+##' @inheritParams xylist.Polygons
+##' @export
+xylist.SpatialPolygons <- function (object, reverse = TRUE, ...)
+{
+ unlist(lapply(object at polygons, xylist.Polygons, reverse=reverse, ...),
+ recursive=FALSE, use.names=FALSE)
+}
+
+##' @rdname xylist
+##' @param reverse logical (\code{TRUE}) indicating if the vertex order of the
+##' \pkg{sp} classes should be reversed to get the \code{xylist}/\code{owin}
+##' convention.
+##' @import sp
+##' @export
+xylist.Polygons <- function (object, reverse = TRUE, ...)
+{
+ lapply(object at Polygons, function (sr) {
+ coords <- coordinates(sr)
+ n <- nrow(coords) - 1L # number of vertices
+ idxs <- if (reverse) seq.int(n,1) else seq_len(n)
+ list(x = coords[idxs,1L], y = coords[idxs,2L])
+ #area = sr at area, hole = sr at hole
+ })
+}
+
+##' @rdname xylist
+##' @import methods
+##' @export
+xylist.Polygon <- function (object, reverse = TRUE, ...)
+ xylist.Polygons(as(object,"Polygons"), reverse=reverse, ...)
+
+##' @rdname xylist
+##' @importFrom grDevices xy.coords
+##' @export
+xylist.default <- function (object, ...) {
+ lapply(object, function (xy) {
+ poly <- xy.coords(xy)[c("x","y")]
+ if (isClosed(poly)) {
+ sel <- seq_len(length(poly$x) - 1L)
+ poly$x <- poly$x[sel]
+ poly$y <- poly$y[sel]
+ }
+ poly
+ })
+}
diff --git a/R/zzz.R b/R/zzz.R
new file mode 100644
index 0000000..fadaba5
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,115 @@
+################################################################################
+### Part of the R package "polyCub".
+### Free software under the terms of the GNU General Public License, version 2,
+### a copy of which is available at http://www.r-project.org/Licenses/.
+###
+### Copyright (C) 2009-2014 Sebastian Meyer
+### Time-stamp: <[zzz.R] 2014-10-24 11:11 (CEST) by SM>
+###
+### Package administration
+################################################################################
+
+
+#' Cubature over Polygonal Domains
+#'
+#' The \R package \pkg{polyCub} provides methods for \strong{cubature}
+#' (numerical integration) \strong{over polygonal domains}.
+#' The function \code{\link{polyCub}()} is the main entry point of the package.
+#' It is a wrapper around the specific cubature methods listed below.
+#'
+#' \describe{
+#' \item{\code{\link{polyCub.midpoint}}:}{
+#' Two-dimensional midpoint rule.
+#' Polygons are converted to binary pixel images
+#' using the \code{\link[spatstat]{as.im.function}} method from package
+#' \pkg{spatstat} (Baddeley and Turner, 2005).
+#' The integral is then obtained as the sum over
+#' (pixel area * f(pixel midpoint)).
+#' }
+#' \item{\code{\link{polyCub.SV}}:}{
+#' Product Gauss cubature as proposed by Sommariva and Vianello (2007).
+#' }
+#' \item{\code{\link{polyCub.iso}}:}{
+#' Efficient adaptive cubature for \emph{isotropic} functions via line
+#' \code{\link{integrate}()} along the polygon boundary, see Meyer and Held
+#' (2014, Supplement B, Section 2.4).
+#' }
+#' \item{\code{\link{polyCub.exact.Gauss}}:}{
+#' Quasi-exact method specific to the integration of the \emph{bivariate Gaussian
+#' density} over polygonal domains. It is based on formulae from Chapter 26 of
+#' the Abramowitz and Stegun (1972) handbook, i.e. triangulation of the
+#' polygonal domain (using \code{\link[gpclib]{tristrip}} of package
+#' \pkg{gpclib}) and appropriate evaluations of
+#' \code{\link[mvtnorm]{pmvnorm}} from package \pkg{mvtnorm}.
+#' Note that there is also a function \code{\link{circleCub.Gauss}}
+#' to perform integration of the \emph{isotropic} Gaussian density over
+#' \emph{circular domains}.
+#' }
+#' }
+#' See Section 3.2 of Meyer (2010) for a more detailed description and benchmark
+#' experiment of some of the above cubature methods (and others).
+#'
+#' @references
+#' Abramowitz, M. and Stegun, I. A. (1972).
+#' Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+#' Tables. New York: Dover Publications.
+#'
+#' Baddeley, A. and Turner, R. (2005).
+#' \pkg{spatstat}: an \R package for analyzing spatial point patterns.
+#' \emph{Journal of Statistical Software}, \bold{12} (6), 1-42.
+#'
+#' Meyer, S. (2010).
+#' Spatio-Temporal Infectious Disease Epidemiology based on Point Processes.
+#' Master's Thesis, LMU Munich.
+#' Available as \url{http://epub.ub.uni-muenchen.de/11703/}.
+#'
+#' Meyer, S. and Held, L. (2014).
+#' Power-law models for infectious disease spread.
+#' \emph{The Annals of Applied Statistics}, \bold{8} (3), 1612-1639.\cr
+#' DOI-Link: \url{http://dx.doi.org/10.1214/14-AOAS743},
+#' \href{http://arxiv.org/abs/1308.5115}{arXiv:1308.5115}
+#'
+#' Sommariva, A. and Vianello, M. (2007).
+#' Product Gauss cubature over polygons based on Green's integration formula.
+#' \emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+#' @docType package
+#' @name polyCub-package
+#' @seealso The packages \pkg{cubature} and \pkg{R2Cuba}, which are more
+#' appropriate for cubature over simple hypercubes.
+NULL
+
+
+.Options <- new.env()
+
+.onLoad <- function (libname, pkgname)
+{
+ .Options$gpclib <- FALSE
+}
+
+gpclibCheck <- function (fatal = TRUE)
+{
+ gpclibOK <- .Options$gpclib
+ if (!gpclibOK && fatal) {
+ message("Note: The gpclib license is accepted by ",
+ sQuote("gpclibPermit()"), ".")
+ stop("acceptance of the gpclib license is required")
+ }
+ gpclibOK
+}
+
+##' \pkg{gpclib} Licence Acceptance
+##'
+##' Similar to the handling in package \pkg{maptools}, these functions
+##' explicitly accept the restricted \pkg{gpclib} licence (commercial use
+##' prohibited) and report its acceptance status, respectively.
+##' \pkg{gpclib} functionality is only required for
+##' \code{\link{polyCub.exact.Gauss}}.
+##' @export
+gpclibPermit <- function ()
+{
+ if (requireNamespace("gpclib")) .Options$gpclib <- TRUE
+ gpclibPermitStatus()
+}
+##' @rdname gpclibPermit
+##' @export
+gpclibPermitStatus <- function () gpclibCheck(fatal=FALSE)
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..c1f2732
--- /dev/null
+++ b/README.md
@@ -0,0 +1,36 @@
+polyCub [`@CRAN`](https://CRAN.R-project.org/package=polyCub)
+============================================================
+
+An `R` package providing methods for **cubature** (numerical integration) **over
+polygonal domains**. Note that for cubature over simple hypercubes, the packages
+[`cubature`](https://CRAN.R-project.org/package=cubature)
+and [`R2Cuba`](https://CRAN.R-project.org/package=R2Cuba)
+might be more appropriate (cf.
+[`CRAN Task View: Numerical Mathematics`](https://CRAN.R-project.org/view=NumericalMathematics)).
+
+The function `polyCub()` is the main entry point of the package. It is a
+wrapper around the following specific cubature methods.
+
+#### General-purpose cubature rules:
+* `polyCub.midpoint()`: Two-dimensional midpoint rule (a simple wrapper around
+ `spatstat::as.im.function`)
+* `polyCub.SV()`: Product Gauss cubature as proposed by
+ [Sommariva and Vianello (2007, *Bit Numerical Mathematics*,
+ **47** (2), 441-453)](http://dx.doi.org/10.1007/s10543-007-0131-2)
+
+#### Cubature rules for specific types of functions:
+* `polyCub.iso()`: Efficient adaptive cubature for *isotropic* functions via
+ line `integrate()` along the polygon boundary
+* `polyCub.exact.Gauss()` and `circleCub.Gauss()`:
+ Quasi-exact methods specific to the integration of the
+ *bivariate Gaussian density* over polygonal and circular domains, respectively
+
+
+A Short History of the Package
+------------------------------
+For any spatio-temporal point process model, the likelihood contains an integral of the conditional intensity function over the observation region. If this is a polygon, analytical solutions are only available for trivial cases of the intensity function -- thus the need of a cubature method over polygonal domains.
+
+My Master's Thesis (2010) on ["Spatio-Temporal Infectious Disease Epidemiology based on Point Processes"](http://epub.ub.uni-muenchen.de/11703/) is the origin of this package. Section 3.2 therein offers a more detailed description and benchmark experiment of some of the above cubature methods (and others).
+
+The implementation then went into the [`surveillance`](https://CRAN.R-project.org/package=surveillance) package, where it is used to fit `twinstim()`, self-exciting two-component spatio-temporal point process models, described in [Meyer et al (2012, *Biometrics*, **68** (2), 607-616)](http://dx.doi.org/10.1111/j.1541-0420.2011.01684.x).
+In May 2013, I decided to move the cubature functionality into a stand-alone package, since it might be useful for other projects as well. Subsequently, the isotropic cubature method `polyCub.iso()` was developed to efficiently estimate point process models with a power-law distance decay of interaction ([Meyer and Held, 2014, *The Annals of Applied Statistics*, **8** (3), 1612-1639](http://dx.doi.org/10.1214/14-AOAS743)).
diff --git a/debian/README.source b/debian/README.source
deleted file mode 100644
index fe3c6e1..0000000
--- a/debian/README.source
+++ /dev/null
@@ -1,8 +0,0 @@
-Explanation for binary files inside source package according to
- http://lists.debian.org/debian-devel/2013/09/msg00332.html
-
-Files: R/sysdata.rda
-Documented: R/polyCub.SV.R
- Cached data
-
- -- Andreas Tille <tille at debian.org> Thu, 07 Sep 2017 21:55:13 +0200
diff --git a/debian/README.test b/debian/README.test
deleted file mode 100644
index f2fd841..0000000
--- a/debian/README.test
+++ /dev/null
@@ -1,21 +0,0 @@
-Notes on how this package can be tested.
-────────────────────────────────────────
-
-To run the unit tests provided by the package you can do
-
- sh run-unit-test
-
-in this directory.
-
-One test of the test suite is requiring the gpclib GNU R
-package with non-free license. It is available from
- http://cran.r-project.org/web/packages/gpclib/
-and the Debian packaging available from
- svn://svn.debian.org/debian-science/packages/R/r-cran-gpclib/trunk/
-could be used to create a Debian package from this code.
-
-Since we can not rely from non-free software the test in question is
-excluded by a quilt patch in the packaging source.
-
- -- Andreas Tille <tille at debian.org> Mon, 23 Jun 2014 15:18:25 +0200
-
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index 12104ed..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,36 +0,0 @@
-r-cran-polycub (0.6.0-1) unstable; urgency=medium
-
- * New upstream version
- * debhelper 10
- * Standards-Version: 4.1.0 (no changes needed)
-
- -- Andreas Tille <tille at debian.org> Thu, 07 Sep 2017 21:55:13 +0200
-
-r-cran-polycub (0.5-2-1) unstable; urgency=medium
-
- * New upstream version
- * Convert to dh-r
- * Canonical homepage for CRAN
- * d/watch: version=4
-
- -- Andreas Tille <tille at debian.org> Mon, 14 Nov 2016 12:07:41 +0100
-
-r-cran-polycub (0.5-1-1) unstable; urgency=medium
-
- * New upstream version
- * cme fix dpkg-control
-
- -- Andreas Tille <tille at debian.org> Sat, 25 Oct 2014 10:34:31 +0200
-
-r-cran-polycub (0.5-0-1) unstable; urgency=medium
-
- * New upstream version
- * Add autopkgtest (skipped gpclib test since gpclib is non-free)
-
- -- Andreas Tille <tille at debian.org> Mon, 23 Jun 2014 15:18:25 +0200
-
-r-cran-polycub (0.4-1-1) unstable; urgency=low
-
- * Initial release (Closes: #736509)
-
- -- Andreas Tille <tille at debian.org> Fri, 24 Jan 2014 11:41:19 +0100
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index f599e28..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-10
diff --git a/debian/control b/debian/control
deleted file mode 100644
index 0693fc6..0000000
--- a/debian/control
+++ /dev/null
@@ -1,33 +0,0 @@
-Source: r-cran-polycub
-Maintainer: Debian Science Team <debian-science-maintainers at lists.alioth.debian.org>
-Uploaders: Andreas Tille <tille at debian.org>
-Section: gnu-r
-Priority: optional
-Build-Depends: debhelper (>= 10),
- dh-r,
- r-base-dev,
- r-cran-sp,
- r-cran-spatstat,
- r-cran-statmod
-Standards-Version: 4.1.0
-Vcs-Browser: https://anonscm.debian.org/viewvc/debian-science/packages/R/trunk/packages/r-cran-polycub/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-science/packages/R/trunk/packages/r-cran-polycub/trunk/
-Homepage: https://cran.r-project.org/package=polyCub
-
-Package: r-cran-polycub
-Architecture: any
-Depends: ${misc:Depends},
- ${shlibs:Depends},
- ${R:Depends}
-Recommends: ${R:Recommends}
-Suggests: ${R:Suggests}
-Description: GNU R Cubature over Polygonal Domains
- PolyCub is a GNU R package providing methods for cubature (numerical
- integration) over polygonal domains. Currently, four cubature methods
- are implemented: the two-dimensional midpoint rule (a simple wrapper
- around spatstat::as.im.function), the product Gauss cubature proposed by
- Sommariva and Vianello (2007), an adaptive cubature for isotropic
- functions via line integrate() along the boundary, and quasi-exact
- methods specific to the integration of the bivariate Gaussian density
- over polygonal and circular domains. For cubature over simple
- hypercubes, the packages cubature and R2Cuba are more appropriate.
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index d4aa8b2..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,29 +0,0 @@
-Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: polyCub
-Upstream-Contact: Sebastian Meyer <Sebastian.Meyer at ifspm.uzh.ch>
-Source: https://cran.r-project.org/package=polyCub
-
-Files: *
-Copyright: (C) 2005-2016 Sebastian Meyer, Michael Hoehle
-License: GPL-2+
-
-Files: debian/*
-Copyright: 2014-2016 Andreas Tille <tille at debian.org>
-License: GPL-2+
-
-License: GPL-2+
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- .
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- .
- On Debian systems, the complete text of the GNU General Public
- License can be found in `/usr/share/common-licenses/GPL-2'.
diff --git a/debian/docs b/debian/docs
deleted file mode 100644
index 3adf0d6..0000000
--- a/debian/docs
+++ /dev/null
@@ -1,3 +0,0 @@
-debian/README.test
-debian/tests/run-unit-test
-tests
diff --git a/debian/patches/series b/debian/patches/series
deleted file mode 100644
index 2270961..0000000
--- a/debian/patches/series
+++ /dev/null
@@ -1 +0,0 @@
-skip_gpclib_test.patch
diff --git a/debian/patches/skip_gpclib_test.patch b/debian/patches/skip_gpclib_test.patch
deleted file mode 100644
index 2131e68..0000000
--- a/debian/patches/skip_gpclib_test.patch
+++ /dev/null
@@ -1,25 +0,0 @@
-Author: Andreas Tille <tille at debian.org>
-Last-Update: Mon, 23 Jun 2014 15:18:25 +0200
-Description: One test of the test suite is requiring the gpclib GNU R
- package with non-free license. It is available from
- http://cran.r-project.org/web/packages/gpclib/
- and the Debian packaging available from
- svn://svn.debian.org/debian-science/packages/R/r-cran-gpclib/trunk/
- could be used to create a Debian package from this code.
- .
- Since we can not rely from non-free software the test in question is
- excluded by this patch.
-
---- a/tests/testthat/test-polyCub.R
-+++ b/tests/testthat/test-polyCub.R
-@@ -18,7 +18,9 @@ sd <- 3
-
- ## exact value of the integral over the _polygonal_ circle
- intExact <- 0.65844436
--if (requireNamespace("mvtnorm") && gpclibPermit()) {
-+# if (requireNamespace("mvtnorm") && gpclibPermit()) {
-+# Skip this test since gpclib is non-free
-+if ( 0 == 1 ) {
- ## run this conditionally since gpclib might not be available on all
- ## platforms (as pointed out by Uwe Ligges, 2014-04-20)
- test_that("polyCub.exact.Gauss returns validated result", {
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 68d9a36..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,4 +0,0 @@
-#!/usr/bin/make -f
-
-%:
- dh $@ --buildsystem R
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/tests/control b/debian/tests/control
deleted file mode 100644
index b044b0c..0000000
--- a/debian/tests/control
+++ /dev/null
@@ -1,3 +0,0 @@
-Tests: run-unit-test
-Depends: @, r-cran-testthat
-Restrictions: allow-stderr
diff --git a/debian/tests/run-unit-test b/debian/tests/run-unit-test
deleted file mode 100644
index db09f19..0000000
--- a/debian/tests/run-unit-test
+++ /dev/null
@@ -1,12 +0,0 @@
-#!/bin/sh -e
-
-pkg=r-cran-polycub
-
-if [ "$ADTTMP" = "" ] ; then
- ADTTMP=`mktemp -d /tmp/${pkg}-test.XXXXXX`
-fi
-cd $ADTTMP
-cp -a /usr/share/doc/${pkg}/tests/* $ADTTMP
-# Make sure we are using C locale to pass all tests
-LC_ALL=C R --no-save < test-all.R
-rm -rf $ADTTMP/*
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index cd18c59..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,2 +0,0 @@
-version=4
-http://cran.r-project.org/src/contrib/polyCub_([-\d.]*)\.tar\.gz
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..1e2d565
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,22 @@
+
+### Outer header
+
+citHeader("To cite", sQuote("polyCub"),
+ "in publications refer to Supplement B (Section 2.4) of:")
+
+
+### power-law paper
+
+bibentry(
+ bibtype = "Article", key = "Meyer+Held_2014",
+ author = "Sebastian Meyer and Leonhard Held",
+ title = "Power-law models for infectious disease spread",
+ journal = "The Annals of Applied Statistics",
+ issn = "1932-6157",
+ year = "2014",
+ volume = "8",
+ number = "3",
+ pages = "1612--1639",
+ doi = "10.1214/14-AOAS743",
+ eprint = "http://arxiv.org/abs/1308.5115"
+)
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
new file mode 100644
index 0000000..4157db5
--- /dev/null
+++ b/inst/NEWS.Rd
@@ -0,0 +1,248 @@
+\newcommand{\CRANpkg}{\href{https://CRAN.R-project.org/package=#1}{\pkg{#1}}}
+%% some pre-defined commands: \R, \code, \acronym, \url, \file, \pkg
+%% Since R 3.2.0, there are some additional system Rd macros available,
+%% e.g., \CRANpkg and \doi. See the definitions in the file
+%% file.path(R.home("share"), "Rd", "macros", "system.Rd")
+
+\name{NEWS}
+\title{News for Package 'polyCub'}
+
+
+\section{Changes in polyCub version 0.6.0 (2017-05-24)}{
+
+ \itemize{
+ \item Added full C-implementation of \code{polyCub.iso()},
+ which is exposed as \code{"polyCub_iso"} for use by other
+ \R packages (notably future versions of \CRANpkg{surveillance})
+ via \samp{LinkingTo: polyCub} and \samp{#include <polyCubAPI.h>}.
+
+ \item Accommodate CRAN checks:
+ add missing import from \pkg{graphics},
+ register native routines and disable symbol search
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.5-2 (2015-02-25)}{
+
+ \itemize{
+ \item \code{polyCub.midpoint()} works directly with input
+ polygons of classes \code{"gpc.poly"} and \code{"SpatialPolygons"},
+ since package \pkg{polyCub} now registers corresponding
+ \code{as.owin}-methods.
+
+ \item \code{polyCub.exact.Gauss()} did not work if the
+ \code{tristrip} of the transformed input polygon contained
+ degenerate triangles (spotted by Ignacio Quintero).
+
+ \item Line integration in \code{polyCub.iso()} could break due to
+ division by zero if the \code{center} point was part of the polygon
+ boundary.
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.5-1 (2014-10-24)}{
+
+ \itemize{
+ \item Nodes and weights for \code{polyCub.SV()} were only cached
+ up to \code{nGQ=59}, not 60 as announced in version 0.5-0. Fixed
+ that which also makes examples truly run without \pkg{statmod}.
+
+ \item In \code{polyCub.SV()}, the new special setting
+ \code{f=NULL} means to only compute nodes and weights.
+
+ \item Internal changes to the \code{"gpc.poly"} converters
+ to accommodate \CRANpkg{spatstat} 1.39-0.
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.5-0 (2014-05-07)}{
+
+ \itemize{
+ \item \code{polyCub.SV()} gained an argument \code{engine} to choose
+ among available implementations. The new and faster C-implementation
+ is the default. There should not be any numerical differences in the
+ result of the cubature.
+
+ \item Package \CRANpkg{statmod} is no longer strictly required
+ (imported). Nodes and weights for Gauss-Legendre quadrature in
+ \code{polyCub.SV()} are now cached in the \pkg{polyCub} package
+ up to \code{nGQ=60}. \pkg{statmod}\code{::gauss.quad} is only
+ queried for a higher number of nodes.
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.4-3 (2014-03-14)}{
+
+ \itemize{
+ \item \code{polyCub.iso()} ...
+ \itemize{
+ \item could not handle additional arguments for
+ \code{integrate()} given in the \code{control} list.
+
+ \item applies the \code{control} arguments also
+ to the numerical approximation of \code{intrfr}.
+ }
+
+ \item The \code{checkintrfr()} function is exported and documented.
+
+ \item Added a \file{CITATION} file.
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.4-2 (2014-02-12)}{
+
+ \itemize{
+ \item \code{plotpolyf()} ...
+
+ \itemize{
+ \item gained an additional argument
+ \code{print.args}, an optional list of arguments passed to
+ \code{print.trellis()} if \code{use.lattice=TRUE}.
+
+ \item passed a \emph{data frame} of coordinates
+ to \code{f} instead of a matrix as documented.
+ }
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.4-1 (2013-12-05)}{
+
+ \itemize{
+ \item This version solely fixes a missing \file{NAMESPACE} import to
+ make package \pkg{polyCub} again compatible with older versions of
+ \CRANpkg{spatstat} (< 1.33-0).
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.4-0 (2013-11-19)}{
+
+ \subsection{INFRASTRUCTURE}{
+ \itemize{
+ \item \CRANpkg{rgeos} (and therefore the GEOS library) is no longer
+ strictly required (moved from Imports to Suggests).
+
+ \item Added \code{coerce}-methods from \code{"Polygons"} (or
+ \code{"SpatialPolygons"} or \code{"Polygon"}) to \code{"owin"}
+ (\code{as(..., "owin")}).
+
+ \item \acronym{S4}-style \code{coerce}-methods between
+ \code{"gpc.poly"} and \code{"Polygons"}/\code{"owin"} have been
+ removed from the package (since we no longer import the formal class
+ \code{"gpc.poly"} from \pkg{gpclib} or \pkg{rgeos}).
+ However, there are two new functions \code{gpc2owin} and
+ \code{owin2gpc} similar to those dropped from \CRANpkg{spatstat}
+ since version 1.34-0.
+
+ \item Moved \code{discpoly()} back to \CRANpkg{surveillance} since
+ it is only used there.
+
+ \item The latter two changes cause \CRANpkg{surveillance} version
+ 1.6-0 to be incompatible with this new version of \pkg{polyCub}.
+ Appropriate modifications have been made in the new version
+ 1.7-0 of \pkg{surveillance}.
+ }
+ }
+
+ \subsection{SPEED-UP \code{polyCub.SV()}}{
+ \itemize{
+ \item thorough optimization of \code{polyCub.SV()}-related code
+ resulted in about 27\% speed-up:
+ \itemize{
+ \item use \code{mapply()} instead of a \code{for}-loop
+ \item avoid \code{cbind()}
+ \item use \code{tcrossprod()}
+ \item less object copying
+ }
+ }
+ }
+
+ \subsection{MINOR CHANGES}{
+ \itemize{
+ \item \code{xylist()} is now exported. It simply extracts
+ polygon coordinates from various spatial classes (with same unifying
+ intention as \code{xy.coords()}).
+
+ \item A \code{polyregion} of class \code{"SpatialPolygons"}
+ of length more than 1 now works in \code{polyCub}-methods.
+
+ \item Use aspect ratio of 1 in \code{plotpolyf()}.
+ }
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.3-1 (2013-08-22)}{
+
+ \itemize{
+ \item This version solely fixes a few typos and a technical note
+ from \command{R CMD check} in the current R development version
+ (also import packages into the \file{NAMESPACE} which are listed
+ in the \dQuote{Depends:} field).
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.3-0 (2013-07-06)}{
+
+ \itemize{
+ \item New cubature method \code{polyCub.iso()} specific to isotropic
+ functions (thanks to Emil Hedevang for the basic idea).
+ \item New function \code{plotpolyf()} to plot a polygonal domain on
+ top of an image of a bivariate function.
+ \item The package now depends on \R >= 2.15.0 (for \code{.rowSums()}).
+ \item The package no longer registers \code{"owin"} as an \acronym{S4}-class
+ since we depend on the \pkg{sp} package which does the job. This
+ avoids a spurious warning (in \code{.simpleDuplicateClass()}) upon
+ package installation.
+ \item In \code{discpoly()}, the argument \code{r} has been renamed
+ to \code{radius}. This is backward compatible by partial argument
+ matching in old code.
+ }
+
+}
+
+
+
+\section{Changes in polyCub version 0.2-0 (2013-05-09)}{
+
+ \itemize{
+ \item This is the initial version of the \pkg{polyCub} package
+ mainly built on functions previously maintained within the
+ \CRANpkg{surveillance} package. These methods for cubature of
+ polygonal domains have been outsourced into this separate
+ \pkg{polyCub} package since they are of general use for other
+ packages as well.
+
+ \item The \pkg{polyCub} package has more documentation and tests,
+ avoids the use of \CRANpkg{gpclib} as far as possible (using
+ \CRANpkg{rgeos} instead), and solves a compatibility issue with
+ package \CRANpkg{maptools} (use \code{setClass("owin")} instead
+ of \code{setOldClass("owin")}).
+ }
+
+}
diff --git a/inst/examples/plotpolyf.R b/inst/examples/plotpolyf.R
new file mode 100644
index 0000000..0946438
--- /dev/null
+++ b/inst/examples/plotpolyf.R
@@ -0,0 +1,11 @@
+### a polygonal domain
+data("letterR", package="spatstat")
+
+### f: isotropic exponential decay
+fr <- function(r, rate=1) dexp(r, rate=rate)
+fcenter <- c(2,3)
+f <- function (s, rate=1) fr(sqrt(rowSums(t(t(s)-fcenter)^2)), rate=rate)
+
+### plot
+plotpolyf(letterR, f, use.lattice=FALSE)
+plotpolyf(letterR, f, use.lattice=TRUE)
diff --git a/inst/examples/polyCub.R b/inst/examples/polyCub.R
new file mode 100644
index 0000000..2a9b564
--- /dev/null
+++ b/inst/examples/polyCub.R
@@ -0,0 +1,55 @@
+### Short comparison of the various cubature methods
+
+## 2D-function to integrate (here: isotropic zero-mean Gaussian density)
+f <- function (s, sigma = 5) exp(-rowSums(s^2)/2/sigma^2) / (2*pi*sigma^2)
+
+## simple polygonal integration domain
+disc.owin <- spatstat::disc(radius=5, centre=c(3,2), npoly=8)
+
+## plot image of the function and integration domain
+plotpolyf(disc.owin, f, xlim=c(-8,8), ylim=c(-8,8))
+
+
+### Two-dimensional midpoint rule
+
+testmidpoint <- function (eps, main=paste("2D midpoint rule with eps =",eps))
+{
+ plotpolyf(disc.owin, f, xlim=c(-8,8), ylim=c(-8,8), use.lattice=FALSE)
+ ## add evaluation points to plot
+ with(spatstat::as.mask(disc.owin, eps=eps),
+ points(expand.grid(xcol, yrow), col=m, pch=20))
+ polyCub.midpoint(disc.owin, f, eps=eps)
+}
+testmidpoint(5)
+testmidpoint(3)
+testmidpoint(0.5)
+testmidpoint(0.2)
+
+
+### Product Gauss cubature using an increasing number of nodes
+
+for (nGQ in c(1:5,10,20,60)) {
+ cat("nGQ =", sprintf("%2i",nGQ), ": ",
+ format(polyCub.SV(disc.owin, f, nGQ=nGQ), digits=16),
+ "\n")
+}
+
+## 'rotation' affects location of nodes
+opar <- par(mfrow=c(1,2))
+polyCub.SV(disc.owin, f, nGQ=2, rotation=FALSE, plot=TRUE)
+polyCub.SV(disc.owin, f, nGQ=2, rotation=TRUE, plot=TRUE)
+par(opar)
+
+
+### Line integration along the boundary for isotropic functions
+
+polyCub.iso(disc.owin, f, center=c(0,0)) # see ?polyCub.iso
+
+
+### Quasi-exact cubature of the bivariate Gaussian density
+### using gpclib::tristrip and mvtnorm::pmvnorm()
+
+if (requireNamespace("mvtnorm") && gpclibPermit()) {
+ print(polyCub.exact.Gauss(disc.owin, mean=c(0,0), Sigma=5^2*diag(2),
+ plot=TRUE), digits=16)
+}
diff --git a/inst/examples/polyCub.iso.R b/inst/examples/polyCub.iso.R
new file mode 100644
index 0000000..506bb2a
--- /dev/null
+++ b/inst/examples/polyCub.iso.R
@@ -0,0 +1,22 @@
+## we use the example polygon and f (exponential decay) from
+example(plotpolyf)
+
+## numerical approximation of 'intrfr'
+(intISOnum <- polyCub.iso(letterR, f, center=fcenter))
+
+## analytical 'intrfr' (recall: f_r(r)=dexp(r), we need int_0^R r*f(r) dr)
+intrfr <- function (R, rate=1) pgamma(R, 2, rate) / rate
+(intISOana <- polyCub.iso(letterR, intrfr=intrfr, center=fcenter))
+
+stopifnot(all.equal(intISOana, intISOnum, check.attributes=FALSE))
+
+
+### polygon area: f(r) = 1, f(x,y) = 1, center does not really matter
+
+intrfr.const <- function (R) R^2/2
+(area.ISO <- polyCub.iso(letterR, intrfr=intrfr.const, center=c(0,0)))
+
+stopifnot(all.equal(spatstat::area.owin(letterR),
+ area.ISO,
+ check.attributes=FALSE))
+## the hole is subtracted correctly
diff --git a/inst/include/polyCubAPI.h b/inst/include/polyCubAPI.h
new file mode 100644
index 0000000..fa4061e
--- /dev/null
+++ b/inst/include/polyCubAPI.h
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * Header file with wrapper functions for the C-routines provided by polyCub
+ *
+ * Copyright (C) 2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+#include <stdlib.h> // NULL
+#include <Rinternals.h> // SEXP
+#include <R_ext/Rdynload.h> // R_GetCCallable
+
+typedef double (*intrfr_fn) (double, double*);
+
+void polyCub_iso(
+ double *x, double *y, // vertex coordinates (open)
+ int *L, // number of vertices
+ intrfr_fn intrfr, // F(R)
+ double *pars, // parameters for F(R)
+ double *center_x, double *center_y, // center of isotropy
+ int *subdivisions, double *epsabs, double *epsrel, // Rdqags options
+ int *stop_on_error, // !=0 means to stop at first ier > 0
+ double *value, double *abserr, int *neval) // results
+{
+ static void(*fun)(double*,double*,int*,intrfr_fn,double*,double*,double*,
+ int*,double*,double*,int*,double*,double*,int*) = NULL;
+ if (fun == NULL)
+ fun = (void(*)(double*,double*,int*,intrfr_fn,double*,double*,double*,
+ int*,double*,double*,int*,double*,double*,int*))
+ R_GetCCallable("polyCub", "polyiso");
+ fun(x, y, L, intrfr, pars, center_x, center_y,
+ subdivisions, epsabs, epsrel, stop_on_error,
+ value, abserr, neval);
+ return;
+}
diff --git a/man/checkintrfr.Rd b/man/checkintrfr.Rd
new file mode 100644
index 0000000..4648441
--- /dev/null
+++ b/man/checkintrfr.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.iso.R
+\name{checkintrfr}
+\alias{checkintrfr}
+\title{Check the Integral of \eqn{r f_r(r)}}
+\usage{
+checkintrfr(intrfr, f, ..., center, control = list(), rs = numeric(0L),
+ tolerance = control$rel.tol)
+}
+\arguments{
+\item{intrfr}{analytical antiderivative of \eqn{r f_r(r)} from 0 to \code{R}
+(first argument, not necessarily named \code{"R"}, must be vectorized).
+If missing, \code{intrfr} is approximated numerically using
+\code{\link{integrate}} configured with \code{control}.}
+
+\item{f}{a two-dimensional real function.
+As its first argument it must take a coordinate matrix, i.e., a
+numeric matrix with two columns, and it must return a numeric vector of
+length the number of coordinates.}
+
+\item{...}{further arguments for \code{f} or \code{intrfr}.}
+
+\item{center}{numeric vector of length 2, the center of isotropy.}
+
+\item{control}{list of arguments passed to \code{\link{integrate}}, the
+quadrature rule used for the line integral along the polygon boundary.}
+
+\item{rs}{numeric vector of upper bounds for which to check the validity of
+\code{intrfr}. If it has length 0, no checks are performed.}
+
+\item{tolerance}{of \code{\link{all.equal.numeric}} when comparing
+\code{intrfr} results with numerical integration. Defaults to the
+relative tolerance used for \code{integrate}.}
+}
+\value{
+The \code{intrfr} function. If it was not supplied, its quadrature
+version using \code{integrate} is returned.
+}
+\description{
+This function is auxiliary to \code{\link{polyCub.iso}}.
+The (analytical) integral of \eqn{r f_r(r)} from 0 to \eqn{R} is checked
+against a numeric approximation using \code{\link{integrate}} for various
+values of the upper bound \eqn{R}. A warning is issued if inconsistencies
+are found.
+}
diff --git a/man/circleCub.Gauss.Rd b/man/circleCub.Gauss.Rd
new file mode 100644
index 0000000..152fb89
--- /dev/null
+++ b/man/circleCub.Gauss.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/circleCub.R
+\name{circleCub.Gauss}
+\alias{circleCub.Gauss}
+\title{Integration of the Isotropic Gaussian Density over Circular Domains}
+\usage{
+circleCub.Gauss(center, r, mean, sd)
+}
+\arguments{
+\item{center}{numeric vector of length 2 (center of the circle).}
+
+\item{r}{numeric (radius of the circle). Several radii may be supplied.}
+
+\item{mean}{numeric vector of length 2
+(mean of the bivariate Gaussian density).}
+
+\item{sd}{numeric (common standard deviation of the isotropic
+Gaussian density in both dimensions).}
+}
+\value{
+The integral value (one for each supplied radius).
+}
+\description{
+This function calculates the integral of the bivariate, isotropic Gaussian
+density (i.e. \eqn{\Sigma} = \code{sd^2*diag(2)}) over circular domains via
+the cumulative distribution function of the (non-central) Chi-Squared
+distribution (\code{pchisq}), cp. Formula 26.3.24 in Abramowitz and Stegun
+(1972).
+}
+\note{
+The non-centrality parameter of the evaluated chi-squared distribution
+equals the squared distance between the \code{mean} and the
+\code{center}. If this becomes too large, the result becomes inaccurate, see
+\code{\link{pchisq}}.
+}
+\examples{
+circleCub.Gauss(center=c(1,2), r=3, mean=c(4,5), sd=6)
+
+if (requireNamespace("mvtnorm") && gpclibPermit()) {
+ ## compare with cubature over a polygonal approximation of a circle
+ disc.poly <- spatstat::disc(radius=3, centre=c(1,2), npoly=32)
+ polyCub.exact.Gauss(disc.poly, mean=c(4,5), Sigma=6^2*diag(2))
+}
+}
+\references{
+Abramowitz, M. and Stegun, I. A. (1972).
+Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+Tables. New York: Dover Publications.
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/coerce-gpc-methods.Rd b/man/coerce-gpc-methods.Rd
new file mode 100644
index 0000000..74ae597
--- /dev/null
+++ b/man/coerce-gpc-methods.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/coerce-gpc-methods.R
+\name{coerce-gpc-methods}
+\alias{coerce-gpc-methods}
+\alias{owin2gpc}
+\alias{gpc2owin}
+\alias{as.owin.gpc.poly}
+\title{Conversion between polygonal \code{"owin"} and \code{"gpc.poly"}}
+\usage{
+owin2gpc(object)
+
+gpc2owin(object, ...)
+
+\method{as.owin}{gpc.poly}(W, ...)
+}
+\arguments{
+\item{object}{an object of class \code{"gpc.poly"} or \code{"owin"},
+respectively.}
+
+\item{...}{further arguments passed to \code{\link{owin}}.}
+
+\item{W}{an object of class \code{"gpc.poly"}.}
+}
+\value{
+The converted polygon of class \code{"gpc.poly"} or \code{"owin"},
+respectively. If neither package \pkg{rgeos} nor \pkg{gpclib} are available,
+\code{owin2gpc} will just return the \code{pts} slot of the
+\code{"gpc.poly"} (no formal class) with a warning.
+}
+\description{
+Package \pkg{polyCub} implements converters between the classes
+\code{"\link[=owin.object]{owin}"} of package \pkg{spatstat} and
+\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} of package \pkg{rgeos}
+(originally from \pkg{gpclib}).
+Support for the \code{"gpc.poly"} class was dropped from
+\pkg{spatstat} as of version 1.34-0.
+}
+\note{
+The converter \code{owin2gpc} requires the package \pkg{rgeos} (or
+\pkg{gpclib}) for the formal class definition of a \code{"gpc.poly"}.
+It will produce vertices ordered according to the \pkg{sp} convention,
+i.e. clockwise for normal boundaries and anticlockwise for holes, where,
+however, the first vertex is \emph{not} repeated!
+}
+\seealso{
+\code{\link{xylist}}, and the package \pkg{rgeos} for
+conversions of \code{"gpc.poly"} objects from and to \pkg{sp}'s
+\code{"\linkS4class{SpatialPolygons}"} class.
+}
+\author{
+Sebastian Meyer
+}
+\keyword{methods}
+\keyword{spatial}
diff --git a/man/coerce-sp-methods.Rd b/man/coerce-sp-methods.Rd
new file mode 100644
index 0000000..353b55f
--- /dev/null
+++ b/man/coerce-sp-methods.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/coerce-sp-methods.R
+\name{coerce-sp-methods}
+\alias{coerce-sp-methods}
+\alias{as.owin.SpatialPolygons}
+\alias{as.owin.Polygons}
+\alias{as.owin.Polygon}
+\alias{coerce,SpatialPolygons,owin-method}
+\alias{coerce,Polygons,owin-method}
+\alias{coerce,Polygon,owin-method}
+\alias{coerce,Polygon,Polygons-method}
+\title{Coerce \code{"SpatialPolygons"} to \code{"owin"}}
+\usage{
+\method{as.owin}{SpatialPolygons}(W, ...)
+
+\method{as.owin}{Polygons}(W, ...)
+
+\method{as.owin}{Polygon}(W, ...)
+}
+\arguments{
+\item{W}{an object of class \code{"SpatialPolygons"},
+\code{"Polygons"}, or \code{"Polygon"}.}
+
+\item{...}{further arguments passed to \code{\link{owin}}.}
+}
+\description{
+Package \pkg{polyCub} implements \code{coerce}-methods
+(\code{as(object, Class)}) to convert \code{"\linkS4class{SpatialPolygons}"}
+(or \code{"\linkS4class{Polygons}"} or \code{"\linkS4class{Polygon}"})
+to \code{"\link[=owin.object]{owin}"}.
+They are also registered as \code{\link{as.owin}}-methods to support
+\code{\link{polyCub.midpoint}}.
+Note that the \pkg{maptools} package contains an alternative implementation
+of coercion from \code{"SpatialPolygons"} to \code{"owin"} (and reverse),
+and \R will use the S4 \code{coerce}-method that was loaded last,
+and prefer the \code{as.owin.SpatialPolygons} S3-method exported from
+\pkg{maptools} if attached.
+}
+\author{
+Sebastian Meyer
+}
+\keyword{methods}
+\keyword{spatial}
diff --git a/man/dotprod.Rd b/man/dotprod.Rd
new file mode 100644
index 0000000..7f66717
--- /dev/null
+++ b/man/dotprod.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{dotprod}
+\alias{dotprod}
+\title{Dot/Scalar Product of Two Vectors}
+\usage{
+dotprod(x, y)
+}
+\arguments{
+\item{x, y}{numeric vectors (of compatible lengths).}
+}
+\value{
+\code{sum(x*y)}
+}
+\description{
+This is nothing else than \code{sum(x*y)}.
+}
+\keyword{internal}
+\keyword{math}
diff --git a/man/gpclibPermit.Rd b/man/gpclibPermit.Rd
new file mode 100644
index 0000000..aae0231
--- /dev/null
+++ b/man/gpclibPermit.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/zzz.R
+\name{gpclibPermit}
+\alias{gpclibPermit}
+\alias{gpclibPermitStatus}
+\title{\pkg{gpclib} Licence Acceptance}
+\usage{
+gpclibPermit()
+
+gpclibPermitStatus()
+}
+\description{
+Similar to the handling in package \pkg{maptools}, these functions
+explicitly accept the restricted \pkg{gpclib} licence (commercial use
+prohibited) and report its acceptance status, respectively.
+\pkg{gpclib} functionality is only required for
+\code{\link{polyCub.exact.Gauss}}.
+}
diff --git a/man/isClosed.Rd b/man/isClosed.Rd
new file mode 100644
index 0000000..568acaa
--- /dev/null
+++ b/man/isClosed.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{isClosed}
+\alias{isClosed}
+\title{Check if Polygon is Closed}
+\usage{
+isClosed(coords)
+}
+\arguments{
+\item{coords}{numeric coordinate matrix. It is interpreted by
+\code{\link{xy.coords}}.}
+}
+\value{
+logical
+}
+\description{
+Check if the first and last coordinates of a coordinate matrix are
+identical.
+}
+\keyword{internal}
+\keyword{spatial}
diff --git a/man/isScalar.Rd b/man/isScalar.Rd
new file mode 100644
index 0000000..c23a9c7
--- /dev/null
+++ b/man/isScalar.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{isScalar}
+\alias{isScalar}
+\title{Checks if Argument is Scalar}
+\usage{
+isScalar(x)
+}
+\arguments{
+\item{x}{any object}
+}
+\value{
+logical
+}
+\description{
+Check if the argument is scalar, i.e. a numeric vector of length 1.
+}
+\keyword{internal}
diff --git a/man/makegrid.Rd b/man/makegrid.Rd
new file mode 100644
index 0000000..d885659
--- /dev/null
+++ b/man/makegrid.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{makegrid}
+\alias{makegrid}
+\title{Constructs Equally-Spaced Grid}
+\usage{
+makegrid(range, n)
+}
+\arguments{
+\item{range}{numeric vector of length 2.}
+
+\item{n}{length of the desired grid, i.e. number of bins + 1.}
+}
+\value{
+the desired grid, a numeric vector of length \code{n} covering
+\code{range}.
+}
+\description{
+Construct an equally-spaced grid given a range and the number of cut points
+(one more than the number of resulting bins).
+This is nothing else than \code{seq(range[1], range[2], length.out=n)}.
+}
+\keyword{internal}
diff --git a/man/plot_polyregion.Rd b/man/plot_polyregion.Rd
new file mode 100644
index 0000000..2992a49
--- /dev/null
+++ b/man/plot_polyregion.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{plot_polyregion}
+\alias{plot_polyregion}
+\title{Plots a Polygonal Domain (of Various Classes)}
+\usage{
+plot_polyregion(polyregion, lwd = 2, add = FALSE)
+}
+\arguments{
+\item{polyregion}{a polygonal domain.
+The following classes are supported: \code{"\link[spatstat]{owin}"},
+\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"},
+\code{"\linkS4class{SpatialPolygons}"}, \code{"\linkS4class{Polygons}"},
+and \code{"\linkS4class{Polygon}"}
+(for these we have an internal \code{\link{xylist}} method).}
+
+\item{lwd}{line width of the polygon edges.}
+
+\item{add}{logical. Add to existing plot?}
+}
+\description{
+Plots a Polygonal Domain (of Various Classes)
+}
diff --git a/man/plotpolyf.Rd b/man/plotpolyf.Rd
new file mode 100644
index 0000000..6c708f2
--- /dev/null
+++ b/man/plotpolyf.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plotpolyf.R
+\name{plotpolyf}
+\alias{plotpolyf}
+\title{Plot Polygonal Domain on Image of Bivariate Function}
+\usage{
+plotpolyf(polyregion, f, ..., npixel = 100, cuts = 15,
+ col = rev(heat.colors(cuts + 1)), lwd = 3, xlim = NULL, ylim = NULL,
+ use.lattice = TRUE, print.args = list())
+}
+\arguments{
+\item{polyregion}{a polygonal domain.
+The following classes are supported: \code{"\link[spatstat]{owin}"},
+\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"},
+\code{"\linkS4class{SpatialPolygons}"}, \code{"\linkS4class{Polygons}"},
+and \code{"\linkS4class{Polygon}"}
+(for these we have an internal \code{\link{xylist}} method).}
+
+\item{f}{a two-dimensional real function.
+As its first argument it must take a coordinate matrix, i.e., a
+numeric matrix with two columns, and it must return a numeric vector of
+length the number of coordinates.}
+
+\item{...}{further arguments for \code{f}.}
+
+\item{npixel}{numeric vector of length 1 or 2 setting the number of pixels
+in each dimension.}
+
+\item{cuts}{number of cut points in the \eqn{z} dimension.
+The range of function values will be divided into \code{cuts+1} levels.}
+
+\item{col}{colour vector used for the function levels.}
+
+\item{lwd}{line width of the polygon edges.}
+
+\item{xlim, ylim}{numeric vectors of length 2 setting the axis limits.
+\code{NULL} means using the bounding box of \code{polyregion}.}
+
+\item{use.lattice}{logical indicating if \pkg{lattice} graphics
+(\code{\link[lattice]{levelplot}}) should be used.}
+
+\item{print.args}{a list of arguments passed to \code{\link{print.trellis}}
+for plotting the produced \code{\link[=trellis.object]{"trellis"}} object
+(given \code{use.lattice = TRUE}). The latter will be returned without
+explicit \code{print}ing if \code{print.args} is not a list.}
+}
+\description{
+Produces a combined plot of a polygonal domain and an image of a bivariate
+function, using either \code{\link[lattice:levelplot]{lattice::levelplot}}
+or \code{\link{image}}.
+}
+\examples{
+### a polygonal domain
+data("letterR", package="spatstat")
+
+### f: isotropic exponential decay
+fr <- function(r, rate=1) dexp(r, rate=rate)
+fcenter <- c(2,3)
+f <- function (s, rate=1) fr(sqrt(rowSums(t(t(s)-fcenter)^2)), rate=rate)
+
+### plot
+plotpolyf(letterR, f, use.lattice=FALSE)
+plotpolyf(letterR, f, use.lattice=TRUE)
+}
+\author{
+Sebastian Meyer
+}
+\keyword{hplot}
diff --git a/man/polyCub-package.Rd b/man/polyCub-package.Rd
new file mode 100644
index 0000000..6d3ec7f
--- /dev/null
+++ b/man/polyCub-package.Rd
@@ -0,0 +1,73 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/zzz.R
+\docType{package}
+\name{polyCub-package}
+\alias{polyCub-package}
+\title{Cubature over Polygonal Domains}
+\description{
+The \R package \pkg{polyCub} provides methods for \strong{cubature}
+(numerical integration) \strong{over polygonal domains}.
+The function \code{\link{polyCub}()} is the main entry point of the package.
+It is a wrapper around the specific cubature methods listed below.
+}
+\details{
+\describe{
+\item{\code{\link{polyCub.midpoint}}:}{
+Two-dimensional midpoint rule.
+Polygons are converted to binary pixel images
+using the \code{\link[spatstat]{as.im.function}} method from package
+\pkg{spatstat} (Baddeley and Turner, 2005).
+The integral is then obtained as the sum over
+(pixel area * f(pixel midpoint)).
+}
+\item{\code{\link{polyCub.SV}}:}{
+Product Gauss cubature as proposed by Sommariva and Vianello (2007).
+}
+\item{\code{\link{polyCub.iso}}:}{
+Efficient adaptive cubature for \emph{isotropic} functions via line
+\code{\link{integrate}()} along the polygon boundary, see Meyer and Held
+(2014, Supplement B, Section 2.4).
+}
+\item{\code{\link{polyCub.exact.Gauss}}:}{
+Quasi-exact method specific to the integration of the \emph{bivariate Gaussian
+density} over polygonal domains. It is based on formulae from Chapter 26 of
+the Abramowitz and Stegun (1972) handbook, i.e. triangulation of the
+polygonal domain (using \code{\link[gpclib]{tristrip}} of package
+\pkg{gpclib}) and appropriate evaluations of
+\code{\link[mvtnorm]{pmvnorm}} from package \pkg{mvtnorm}.
+Note that there is also a function \code{\link{circleCub.Gauss}}
+to perform integration of the \emph{isotropic} Gaussian density over
+\emph{circular domains}.
+}
+}
+See Section 3.2 of Meyer (2010) for a more detailed description and benchmark
+experiment of some of the above cubature methods (and others).
+}
+\references{
+Abramowitz, M. and Stegun, I. A. (1972).
+Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+Tables. New York: Dover Publications.
+
+Baddeley, A. and Turner, R. (2005).
+\pkg{spatstat}: an \R package for analyzing spatial point patterns.
+\emph{Journal of Statistical Software}, \bold{12} (6), 1-42.
+
+Meyer, S. (2010).
+Spatio-Temporal Infectious Disease Epidemiology based on Point Processes.
+Master's Thesis, LMU Munich.
+Available as \url{http://epub.ub.uni-muenchen.de/11703/}.
+
+Meyer, S. and Held, L. (2014).
+Power-law models for infectious disease spread.
+\emph{The Annals of Applied Statistics}, \bold{8} (3), 1612-1639.\cr
+DOI-Link: \url{http://dx.doi.org/10.1214/14-AOAS743},
+\href{http://arxiv.org/abs/1308.5115}{arXiv:1308.5115}
+
+Sommariva, A. and Vianello, M. (2007).
+Product Gauss cubature over polygons based on Green's integration formula.
+\emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+}
+\seealso{
+The packages \pkg{cubature} and \pkg{R2Cuba}, which are more
+appropriate for cubature over simple hypercubes.
+}
diff --git a/man/polyCub.Rd b/man/polyCub.Rd
new file mode 100644
index 0000000..7dcb6f5
--- /dev/null
+++ b/man/polyCub.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.R
+\name{polyCub}
+\alias{polyCub}
+\title{Wrapper Function for the Various Cubature Methods}
+\usage{
+polyCub(polyregion, f, method = c("SV", "midpoint", "iso", "exact.Gauss"),
+ ..., plot = FALSE)
+}
+\arguments{
+\item{polyregion}{a polygonal integration domain.
+The supported classes depend on the specific method, however, the
+\code{"\link[spatstat]{owin}"} class from package \pkg{spatstat} works for
+all methods, as well should a \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"}
+polygon (but see the comments in \code{help("\link{coerce-methods}")}).}
+
+\item{f}{two-dimensional function to be integrated.
+As its first argument the function must take a coordinate matrix, i.e. a
+numeric matrix with two columns. For the \code{"exact.Gauss"} \code{method},
+\code{f} is ignored since it is specific to the bivariate normal density.}
+
+\item{method}{choose one of the implemented cubature methods (partial
+argument matching is applied), see \code{help("\link{polyCub-package}")}
+for an overview. Defaults to using the product Gauss cubature
+implemented in \code{\link{polyCub.SV}}.}
+
+\item{...}{arguments of \code{f} or of the specific \code{method}.}
+
+\item{plot}{logical indicating if an illustrative plot of the numerical
+integration should be produced.}
+}
+\value{
+The approximated integral of \code{f} over \code{polyregion}.
+}
+\description{
+Instead of calling one of the specific cubature methods of this package, the
+wrapper function \code{polyCub} may be used together with the \code{method}
+argument.
+}
+\examples{
+### Short comparison of the various cubature methods
+
+## 2D-function to integrate (here: isotropic zero-mean Gaussian density)
+f <- function (s, sigma = 5) exp(-rowSums(s^2)/2/sigma^2) / (2*pi*sigma^2)
+
+## simple polygonal integration domain
+disc.owin <- spatstat::disc(radius=5, centre=c(3,2), npoly=8)
+
+## plot image of the function and integration domain
+plotpolyf(disc.owin, f, xlim=c(-8,8), ylim=c(-8,8))
+
+
+### Two-dimensional midpoint rule
+
+testmidpoint <- function (eps, main=paste("2D midpoint rule with eps =",eps))
+{
+ plotpolyf(disc.owin, f, xlim=c(-8,8), ylim=c(-8,8), use.lattice=FALSE)
+ ## add evaluation points to plot
+ with(spatstat::as.mask(disc.owin, eps=eps),
+ points(expand.grid(xcol, yrow), col=m, pch=20))
+ polyCub.midpoint(disc.owin, f, eps=eps)
+}
+testmidpoint(5)
+testmidpoint(3)
+testmidpoint(0.5)
+testmidpoint(0.2)
+
+
+### Product Gauss cubature using an increasing number of nodes
+
+for (nGQ in c(1:5,10,20,60)) {
+ cat("nGQ =", sprintf("\%2i",nGQ), ": ",
+ format(polyCub.SV(disc.owin, f, nGQ=nGQ), digits=16),
+ "\\n")
+}
+
+## 'rotation' affects location of nodes
+opar <- par(mfrow=c(1,2))
+polyCub.SV(disc.owin, f, nGQ=2, rotation=FALSE, plot=TRUE)
+polyCub.SV(disc.owin, f, nGQ=2, rotation=TRUE, plot=TRUE)
+par(opar)
+
+
+### Line integration along the boundary for isotropic functions
+
+polyCub.iso(disc.owin, f, center=c(0,0)) # see ?polyCub.iso
+
+
+### Quasi-exact cubature of the bivariate Gaussian density
+### using gpclib::tristrip and mvtnorm::pmvnorm()
+
+if (requireNamespace("mvtnorm") && gpclibPermit()) {
+ print(polyCub.exact.Gauss(disc.owin, mean=c(0,0), Sigma=5^2*diag(2),
+ plot=TRUE), digits=16)
+}
+}
+\seealso{
+Other polyCub-methods: \code{\link{polyCub.SV}},
+ \code{\link{polyCub.exact.Gauss}},
+ \code{\link{polyCub.iso}}, \code{\link{polyCub.midpoint}}
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/polyCub.SV.Rd b/man/polyCub.SV.Rd
new file mode 100644
index 0000000..eae708b
--- /dev/null
+++ b/man/polyCub.SV.Rd
@@ -0,0 +1,97 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.SV.R
+\name{polyCub.SV}
+\alias{polyCub.SV}
+\title{Product Gauss Cubature over Polygonal Domains}
+\usage{
+polyCub.SV(polyregion, f, ..., nGQ = 20, alpha = NULL, rotation = FALSE,
+ engine = "C", plot = FALSE)
+}
+\arguments{
+\item{polyregion}{a polygonal domain.
+The following classes are supported: \code{"\link[spatstat]{owin}"},
+\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"},
+\code{"\linkS4class{SpatialPolygons}"}, \code{"\linkS4class{Polygons}"},
+and \code{"\linkS4class{Polygon}"}
+(for these we have an internal \code{\link{xylist}} method).}
+
+\item{f}{a two-dimensional real function (or \code{NULL} to only compute
+nodes and weights).
+As its first argument it must take a coordinate matrix, i.e., a
+numeric matrix with two columns, and it must return a numeric vector of
+length the number of coordinates.}
+
+\item{...}{further arguments for \code{f}.}
+
+\item{nGQ}{degree of the one-dimensional Gauss-Legendre quadrature rule
+(default: 20) as implemented in function \code{\link[statmod]{gauss.quad}}
+of package \pkg{statmod}. Nodes and weights up to \code{nGQ=60} are cached
+in \pkg{polyCub}, for larger degrees \pkg{statmod} is required.}
+
+\item{alpha}{base-line of the (rotated) polygon at \eqn{x = \alpha} (see
+Sommariva and Vianello (2007) for an explication). If \code{NULL} (default),
+the midpoint of the x-range of each polygon is chosen if no \code{rotation}
+is performed, and otherwise the \eqn{x}-coordinate of the rotated point
+\code{"P"} (see \code{rotation}). If \code{f} has its maximum value at the
+origin \eqn{(0,0)}, e.g., the bivariate Gaussian density with zero mean,
+\code{alpha = 0} is a reasonable choice.}
+
+\item{rotation}{logical (default: \code{FALSE}) or a list of points
+\code{"P"} and \code{"Q"} describing the preferred direction. If
+\code{TRUE}, the polygon is rotated according to the vertices \code{"P"} and
+\code{"Q"}, which are farthest apart (see Sommariva and Vianello, 2007). For
+convex polygons, this rotation guarantees that all nodes fall inside the
+polygon.}
+
+\item{engine}{character string specifying the implementation to use.
+Up to \pkg{polyCub} version 0.4-3, the two-dimensional nodes and weights
+were computed by \R functions and these are still available by setting
+\code{engine = "R"}.
+The new C-implementation is now the default (\code{engine = "C"}) and
+requires approximately 30\% less computation time.\cr
+The special setting \code{engine = "C+reduce"} will discard redundant nodes
+at (0,0) with zero weight resulting from edges on the base-line
+\eqn{x = \alpha} or orthogonal to it.
+This extra cleaning is only worth its cost for computationally intensive
+functions \code{f} over polygons which really have some edges on the
+baseline or parallel to the x-axis. Note that the old \R
+implementation does not have such unset zero nodes and weights.}
+
+\item{plot}{logical indicating if an illustrative plot of the numerical
+integration should be produced.}
+}
+\value{
+The approximated value of the integral of \code{f} over
+\code{polyregion}.\cr
+In the case \code{f = NULL}, only the computed nodes and weights are
+returned in a list of length the number of polygons of \code{polyregion},
+where each component is a list with \code{nodes} (a numeric matrix with
+two columns), \code{weights} (a numeric vector of length
+\code{nrow(nodes)}), the rotation \code{angle}, and \code{alpha}.
+}
+\description{
+Product Gauss cubature over polygons as proposed by
+Sommariva and Vianello (2007).
+}
+\examples{
+# see example(polyCub)
+}
+\references{
+Sommariva, A. and Vianello, M. (2007).
+Product Gauss cubature over polygons based on Green's integration formula.
+\emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+}
+\seealso{
+Other polyCub-methods: \code{\link{polyCub.exact.Gauss}},
+ \code{\link{polyCub.iso}},
+ \code{\link{polyCub.midpoint}}, \code{\link{polyCub}}
+}
+\author{
+Sebastian Meyer\cr
+The product Gauss cubature is based on the
+original \acronym{MATLAB} implementation \code{polygauss} by Sommariva and
+Vianello (2007), which is available under the GNU GPL (>=2) license from
+\url{http://www.math.unipd.it/~alvise/software.html}.
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/polyCub.exact.Gauss.Rd b/man/polyCub.exact.Gauss.Rd
new file mode 100644
index 0000000..7d2242b
--- /dev/null
+++ b/man/polyCub.exact.Gauss.Rd
@@ -0,0 +1,76 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.exact.Gauss.R
+\name{polyCub.exact.Gauss}
+\alias{polyCub.exact.Gauss}
+\title{Quasi-Exact Cubature of the Bivariate Normal Density}
+\usage{
+polyCub.exact.Gauss(polyregion, mean = c(0, 0), Sigma = diag(2),
+ plot = FALSE)
+}
+\arguments{
+\item{polyregion}{a \code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} polygon or
+something that can be coerced to this class, e.g., an \code{"owin"} polygon
+(converted via \code{\link{owin2gpc}} and -- given \pkg{rgeos} is available
+-- \code{"SpatialPolygons"} also work.}
+
+\item{mean, Sigma}{mean and covariance matrix of the bivariate normal density
+to be integrated.}
+
+\item{plot}{logical indicating if an illustrative plot of the numerical
+integration should be produced. Note that the \code{polyregion} will be
+transformed (shifted and scaled).}
+}
+\value{
+The integral of the bivariate normal density over \code{polyregion}.
+Two attributes are appended to the integral value:
+\item{nEval}{
+number of triangles over which the standard bivariate normal density had to
+be integrated, i.e. number of calls to \code{\link[mvtnorm]{pmvnorm}} and
+\code{\link[stats]{pnorm}}, the former of which being the most time-consuming
+operation.
+}
+\item{error}{
+Approximate absolute integration error stemming from the error introduced by
+the \code{nEval} \code{\link[mvtnorm]{pmvnorm}} evaluations.
+For this reason, the cubature method is in fact only
+quasi-exact (as is the \code{pmvnorm} function).
+}
+}
+\description{
+Integration is based on triangulation of the (transformed) polygonal domain
+and formulae from the
+Abramowitz and Stegun (1972) handbook (Section 26.9, Example 9, pp. 956f.).
+This method is quite cumbersome because the A&S formula is only for triangles
+where one vertex is the origin (0,0). For each triangle of the
+\code{\link[gpclib]{tristrip}} we have to check in which of the 6 outer
+regions of the triangle the origin (0,0) lies and adapt the signs in the
+formula appropriately: \eqn{(AOB+BOC-AOC)} or \eqn{(AOB-AOC-BOC)} or
+\eqn{(AOB+AOC-BOC)} or \eqn{(AOC+BOC-AOB)} or \ldots.
+However, the most time consuming step is the
+evaluation of \code{\link[mvtnorm]{pmvnorm}}.
+}
+\note{
+The package \pkg{gpclib} (which is required to produce the
+\code{tristrip}, since this is not yet implemented in \pkg{rgeos})
+has a restricted license (commercial use prohibited).
+It has to be accepted explicitly via
+\code{\link{gpclibPermit}()} prior to using \code{polyCub.exact.Gauss}.
+}
+\examples{
+# see example(polyCub)
+}
+\references{
+Abramowitz, M. and Stegun, I. A. (1972).
+Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
+Tables. New York: Dover Publications.
+}
+\seealso{
+\code{\link{circleCub.Gauss}} for quasi-exact cubature of the
+isotropic Gaussian density over a circular domain.
+
+Other polyCub-methods: \code{\link{polyCub.SV}},
+ \code{\link{polyCub.iso}},
+ \code{\link{polyCub.midpoint}}, \code{\link{polyCub}}
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/polyCub.iso.Rd b/man/polyCub.iso.Rd
new file mode 100644
index 0000000..88d9577
--- /dev/null
+++ b/man/polyCub.iso.Rd
@@ -0,0 +1,134 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.iso.R
+\name{polyCub.iso}
+\alias{polyCub.iso}
+\alias{.polyCub.iso}
+\title{Cubature of Isotropic Functions over Polygonal Domains}
+\usage{
+polyCub.iso(polyregion, f, intrfr, ..., center, control = list(),
+ check.intrfr = FALSE, plot = FALSE)
+
+.polyCub.iso(polys, intrfr, ..., center, control = list(),
+ .witherror = FALSE)
+}
+\arguments{
+\item{polyregion}{a polygonal domain.
+The following classes are supported: \code{"\link[spatstat]{owin}"},
+\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"},
+\code{"\linkS4class{SpatialPolygons}"}, \code{"\linkS4class{Polygons}"},
+and \code{"\linkS4class{Polygon}"}
+(for these we have an internal \code{\link{xylist}} method).}
+
+\item{f}{a two-dimensional real function.
+As its first argument it must take a coordinate matrix, i.e., a
+numeric matrix with two columns, and it must return a numeric vector of
+length the number of coordinates.}
+
+\item{intrfr}{analytical antiderivative of \eqn{r f_r(r)} from 0 to \code{R}
+(first argument, not necessarily named \code{"R"}, must be vectorized).
+If missing, \code{intrfr} is approximated numerically using
+\code{\link{integrate}} configured with \code{control}.}
+
+\item{...}{further arguments for \code{f} or \code{intrfr}.}
+
+\item{center}{numeric vector of length 2, the center of isotropy.}
+
+\item{control}{list of arguments passed to \code{\link{integrate}}, the
+quadrature rule used for the line integral along the polygon boundary.}
+
+\item{check.intrfr}{logical (or numeric vector) indicating if
+(for which \code{r}'s) the supplied \code{intrfr} function should be
+checked against a numeric approximation. This check requires \code{f}
+to be specified. If \code{TRUE}, the set of test
+\code{r}'s defaults to a \code{\link{seq}} of length 20 from 1 to
+the maximum absolute x or y coordinate of any edge of the \code{polyregion}.}
+
+\item{plot}{logical indicating if an image of the function should be plotted
+together with the polygonal domain, i.e.,
+\code{\link{plotpolyf}(polyregion, f, \dots)}.}
+
+\item{polys}{something like \code{owin$bdry}, but see \code{\link{xylist}}.}
+
+\item{.witherror}{logical indicating if an upper bound for the absolute
+integration error should be attached as an attribute to the result?}
+}
+\value{
+The approximate integral of the isotropic function
+\code{f} over \code{polyregion}.\cr
+If the \code{intrfr} function is provided (which is assumed to be exact), an
+upper bound for the absolute integration error is appended to the result as
+attribute \code{"abs.error"}. It equals the sum of the absolute errors
+reported by all \code{\link{integrate}} calls
+(there is one for each edge of \code{polyregion}).
+}
+\description{
+Conducts numerical integration of a two-dimensional isotropic function
+\eqn{f(x,y) = f_r(||(x,y)-\boldsymbol{\mu}||)}{f(x,y) = f_r(||(x,y)-\mu||)},
+with \eqn{\mu} being the center of isotropy, over a polygonal domain.
+It internally solves a line integral along the polygon boundary using
+\code{\link{integrate}} where the integrand requires the antiderivative of
+\eqn{r f_r(r)}), which ideally is analytically available and supplied to the
+function as argument \code{intrfr}.
+The two-dimensional integration problem thereby reduces to an efficient
+adaptive quadrature in one dimension.
+See Meyer and Held (2014, Supplement B, Section 2.4) for mathematical
+details.
+
+\code{.polyCub.iso} is a \dQuote{bare-bone} version of \code{polyCub.iso}.
+}
+\examples{
+## we use the example polygon and f (exponential decay) from
+example(plotpolyf)
+
+## numerical approximation of 'intrfr'
+(intISOnum <- polyCub.iso(letterR, f, center=fcenter))
+
+## analytical 'intrfr' (recall: f_r(r)=dexp(r), we need int_0^R r*f(r) dr)
+intrfr <- function (R, rate=1) pgamma(R, 2, rate) / rate
+(intISOana <- polyCub.iso(letterR, intrfr=intrfr, center=fcenter))
+
+stopifnot(all.equal(intISOana, intISOnum, check.attributes=FALSE))
+
+
+### polygon area: f(r) = 1, f(x,y) = 1, center does not really matter
+
+intrfr.const <- function (R) R^2/2
+(area.ISO <- polyCub.iso(letterR, intrfr=intrfr.const, center=c(0,0)))
+
+stopifnot(all.equal(spatstat::area.owin(letterR),
+ area.ISO,
+ check.attributes=FALSE))
+## the hole is subtracted correctly
+}
+\references{
+Hedevang, E. (2013). Personal communication at the Summer School on Topics in
+Space-Time Modeling and Inference (May 2013, Aalborg, Denmark).
+
+Meyer, S. and Held, L. (2014).
+Power-law models for infectious disease spread.
+\emph{The Annals of Applied Statistics}, \bold{8} (3), 1612-1639.\cr
+DOI-Link: \url{http://dx.doi.org/10.1214/14-AOAS743},
+\href{http://arxiv.org/abs/1308.5115}{arXiv:1308.5115}
+}
+\seealso{
+\code{system.file("include", "polyCubAPI.h", package = "polyCub")}
+for a full C-implementation of this cubature method (for a \emph{single}
+polygon). The corresponding C-routine \code{polyCub_iso} can be used by
+other \R packages, notably \pkg{surveillance}, via \samp{LinkingTo: polyCub}
+(in the \file{DESCRIPTION}) and \samp{#include <polyCubAPI.h>} (in suitable
+\file{/src} files). Note that the \code{intrfr} function must then also be
+supplied as a C-routine. An example can be found in the package tests.
+
+Other polyCub-methods: \code{\link{polyCub.SV}},
+ \code{\link{polyCub.exact.Gauss}},
+ \code{\link{polyCub.midpoint}}, \code{\link{polyCub}}
+}
+\author{
+Sebastian Meyer
+
+The basic mathematical formulation of this efficient integration for radially
+symmetric functions was ascertained with great support by
+Emil Hedevang (2013), Dept. of Mathematics, Aarhus University, Denmark.
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/polyCub.midpoint.Rd b/man/polyCub.midpoint.Rd
new file mode 100644
index 0000000..5ba617a
--- /dev/null
+++ b/man/polyCub.midpoint.Rd
@@ -0,0 +1,61 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.midpoint.R
+\name{polyCub.midpoint}
+\alias{polyCub.midpoint}
+\title{Two-Dimensional Midpoint Rule}
+\usage{
+polyCub.midpoint(polyregion, f, ..., eps = NULL, dimyx = NULL,
+ plot = FALSE)
+}
+\arguments{
+\item{polyregion}{a polygonal integration domain.
+It can be any object coercible to the \pkg{spatstat} class
+\code{"\link[spatstat]{owin}"} via a corresponding
+\code{\link[spatstat]{as.owin}}-method.
+Note that this includes polygons of the classes \code{"gpc.poly"} and
+\code{"\linkS4class{SpatialPolygons}"}, because \pkg{polyCub} defines
+methods \code{\link{as.owin.gpc.poly}} and
+\code{\link{as.owin.SpatialPolygons}}, respectively.}
+
+\item{f}{a two-dimensional real function.
+As its first argument it must take a coordinate matrix, i.e., a
+numeric matrix with two columns, and it must return a numeric vector of
+length the number of coordinates.}
+
+\item{...}{further arguments for \code{f}.}
+
+\item{eps}{width and height of the pixels (squares),
+see \code{\link[spatstat]{as.mask}}.}
+
+\item{dimyx}{number of subdivisions in each dimension,
+see \code{\link[spatstat]{as.mask}}.}
+
+\item{plot}{logical indicating if an illustrative plot of the numerical
+integration should be produced.}
+}
+\value{
+The approximated value of the integral of \code{f} over
+\code{polyregion}.
+}
+\description{
+The surface is converted to a binary pixel image
+using the \code{\link[spatstat]{as.im.function}} method from package
+\pkg{spatstat} (Baddeley and Turner, 2005).
+The integral under the surface is then approximated as the
+sum over (pixel area * f(pixel midpoint)).
+}
+\examples{
+# see example(polyCub)
+}
+\references{
+Baddeley, A. and Turner, R. (2005).
+\pkg{spatstat}: an \R package for analyzing spatial point patterns.
+\emph{Journal of Statistical Software}, \bold{12} (6), 1-42.
+}
+\seealso{
+Other polyCub-methods: \code{\link{polyCub.SV}},
+ \code{\link{polyCub.exact.Gauss}},
+ \code{\link{polyCub.iso}}, \code{\link{polyCub}}
+}
+\keyword{math}
+\keyword{spatial}
diff --git a/man/polygauss.Rd b/man/polygauss.Rd
new file mode 100644
index 0000000..f772be9
--- /dev/null
+++ b/man/polygauss.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/polyCub.SV.R
+\name{polygauss}
+\alias{polygauss}
+\title{Calculate 2D Nodes and Weights of the Product Gauss Cubature}
+\usage{
+polygauss(xy, nw_MN, alpha = NULL, rotation = FALSE, engine = "C")
+}
+\arguments{
+\item{xy}{list with elements \code{"x"} and \code{"y"} containing the
+polygon vertices in \emph{anticlockwise} order (otherwise the result of the
+cubature will have a negative sign) with first vertex not repeated at the
+end (like \code{owin.object$bdry}).}
+
+\item{nw_MN}{unnamed list of nodes and weights of one-dimensional Gauss
+quadrature rules of degrees \eqn{N} and \eqn{M=N+1} (as returned by
+\code{\link[statmod]{gauss.quad}}): \code{list(s_M, w_M, s_N, w_N)}.}
+
+\item{alpha}{base-line of the (rotated) polygon at \eqn{x = \alpha} (see
+Sommariva and Vianello (2007) for an explication). If \code{NULL} (default),
+the midpoint of the x-range of each polygon is chosen if no \code{rotation}
+is performed, and otherwise the \eqn{x}-coordinate of the rotated point
+\code{"P"} (see \code{rotation}). If \code{f} has its maximum value at the
+origin \eqn{(0,0)}, e.g., the bivariate Gaussian density with zero mean,
+\code{alpha = 0} is a reasonable choice.}
+
+\item{rotation}{logical (default: \code{FALSE}) or a list of points
+\code{"P"} and \code{"Q"} describing the preferred direction. If
+\code{TRUE}, the polygon is rotated according to the vertices \code{"P"} and
+\code{"Q"}, which are farthest apart (see Sommariva and Vianello, 2007). For
+convex polygons, this rotation guarantees that all nodes fall inside the
+polygon.}
+
+\item{engine}{character string specifying the implementation to use.
+Up to \pkg{polyCub} version 0.4-3, the two-dimensional nodes and weights
+were computed by \R functions and these are still available by setting
+\code{engine = "R"}.
+The new C-implementation is now the default (\code{engine = "C"}) and
+requires approximately 30\% less computation time.\cr
+The special setting \code{engine = "C+reduce"} will discard redundant nodes
+at (0,0) with zero weight resulting from edges on the base-line
+\eqn{x = \alpha} or orthogonal to it.
+This extra cleaning is only worth its cost for computationally intensive
+functions \code{f} over polygons which really have some edges on the
+baseline or parallel to the x-axis. Note that the old \R
+implementation does not have such unset zero nodes and weights.}
+}
+\description{
+Calculate 2D Nodes and Weights of the Product Gauss Cubature
+}
+\references{
+Sommariva, A. and Vianello, M. (2007):
+Product Gauss cubature over polygons based on Green's integration formula.
+\emph{Bit Numerical Mathematics}, \bold{47} (2), 441-453.
+}
+\keyword{internal}
diff --git a/man/vecnorm.Rd b/man/vecnorm.Rd
new file mode 100644
index 0000000..8e232cc
--- /dev/null
+++ b/man/vecnorm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tools.R
+\name{vecnorm}
+\alias{vecnorm}
+\title{Euclidean Vector Norm (Length)}
+\usage{
+vecnorm(x)
+}
+\arguments{
+\item{x}{numeric vector.}
+}
+\value{
+\code{sqrt(sum(x^2))}
+}
+\description{
+This is nothing else than \code{sqrt(sum(x^2))}.
+}
+\keyword{internal}
+\keyword{math}
diff --git a/man/xylist.Rd b/man/xylist.Rd
new file mode 100644
index 0000000..113eef9
--- /dev/null
+++ b/man/xylist.Rd
@@ -0,0 +1,87 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/xylist.R
+\name{xylist}
+\alias{xylist}
+\alias{xylist.owin}
+\alias{xylist.gpc.poly}
+\alias{xylist.SpatialPolygons}
+\alias{xylist.Polygons}
+\alias{xylist.Polygon}
+\alias{xylist.default}
+\title{Convert Various Polygon Classes to a Simple List of Vertices}
+\usage{
+xylist(object, ...)
+
+\method{xylist}{owin}(object, ...)
+
+\method{xylist}{gpc.poly}(object, ...)
+
+\method{xylist}{SpatialPolygons}(object, reverse = TRUE, ...)
+
+\method{xylist}{Polygons}(object, reverse = TRUE, ...)
+
+\method{xylist}{Polygon}(object, reverse = TRUE, ...)
+
+\method{xylist}{default}(object, ...)
+}
+\arguments{
+\item{object}{an object of one of the supported spatial classes.}
+
+\item{...}{(unused) argument of the generic.}
+
+\item{reverse}{logical (\code{TRUE}) indicating if the vertex order of the
+\pkg{sp} classes should be reversed to get the \code{xylist}/\code{owin}
+convention.}
+}
+\value{
+Applying \code{xylist} to a polygon object, one gets a simple list,
+where each component (polygon) is a list of \code{"x"} and \code{"y"}
+coordinates. These represent vertex coordinates following \pkg{spatstat}'s
+\code{"owin"} convention (anticlockwise order without repeating any vertex).
+The opposite vertex order can be retained for the \pkg{sp}-classes
+by the non-default use with \code{reverse=FALSE}.
+}
+\description{
+Different packages concerned with spatial data use different polygon
+specifications, which sometimes becomes very confusing (see Details below).
+To be compatible with the various polygon classes, package \pkg{polyCub}
+uses an S3 class \code{"xylist"}, which represents
+polygons by their core feature only, a list of lists of vertex coordinates
+(see the "Value" section below).
+The generic function \code{xylist} can deal with the
+following polygon classes:
+\itemize{
+\item{\code{"\link[=owin.object]{owin}"} from package \pkg{spatstat}}
+\item{\code{"\link[rgeos:gpc.poly-class]{gpc.poly}"} from package
+\pkg{rgeos} (or \pkg{gpclib})}
+\item{\code{"\linkS4class{Polygons}"} from package \pkg{sp}
+(as well as \code{"\linkS4class{Polygon}"} and
+\code{"\linkS4class{SpatialPolygons}"})}
+}
+The (somehow useless) default \code{xylist}-method
+does not perform any transformation but only ensures that the polygons are
+not closed (first vertex not repeated).
+}
+\details{
+Different packages concerned with spatial data use different polygon
+specifications with respect to:
+\itemize{
+\item{do we repeat the first vertex?}
+\item{which direction represents holes?}
+}
+Package overview:
+\describe{
+\item{\pkg{sp}:}{\emph{Repeat} first vertex at the end (closed),
+anticlockwise = hole, clockwise = normal boundary}
+\item{\pkg{spatstat}:}{do \emph{not repeat} first vertex,
+anticlockwise = normal boundary, clockwise = hole. This convention is also
+used in \code{xylist}.}
+\item{\pkg{gpclib}:}{Unfortunately, there seems to be no convention
+for the specification of polygons of class \code{"gpc.poly"}.}
+}
+}
+\author{
+Sebastian Meyer
+}
+\keyword{methods}
+\keyword{spatial}
diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..f07dbb1
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,29 @@
+/*******************************************************************************
+ * Registering native routines (entry points in compiled code)
+ *
+ * Copyright (C) 2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+#include "polyCub.SV.h"
+#include "polyCub.iso.h"
+
+static const R_CMethodDef CEntries[] = {
+ {"C_polygauss", (DL_FUNC) &C_polygauss, 13},
+ {NULL, NULL, 0}
+};
+
+void R_init_polyCub(DllInfo *dll)
+{
+ R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
+ R_useDynamicSymbols(dll, FALSE);
+ //R_forceSymbols(dll, TRUE); // would require R >= 3.0.0
+
+ R_RegisterCCallable("polyCub", "polyiso", (DL_FUNC) &polyiso);
+}
diff --git a/src/polyCub.SV.c b/src/polyCub.SV.c
new file mode 100644
index 0000000..7171460
--- /dev/null
+++ b/src/polyCub.SV.c
@@ -0,0 +1,85 @@
+/*******************************************************************************
+ * C-version of .polygauss.side()
+ *
+ * Copyright (C) 2014,2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+#include "polyCub.SV.h"
+
+static void C_polygauss_side(
+ double *x1, double *y1, double *x2, double *y2,
+ double *s_loc, double *w_loc, double *s_N, double *w_N,
+ double *alpha,
+ int *loc, int *N, // lengths (loc is M=N+1 or N)
+ // *loc * *N nodes and weights will be computed
+ double *nodes_x, double *nodes_y, double *weights)
+{
+ double half_pt_x = (*x1 + *x2) / 2.0;
+ double half_length_x = (*x2 - *x1) / 2.0;
+ double half_pt_y = (*y1 + *y2) / 2.0;
+ double half_length_y = (*y2 - *y1) / 2.0;
+
+ double x_gauss_side, y_gauss_side, scaling_fact_minus;
+ int idx;
+ for (int i = 0; i < *loc; i++) {
+ // GAUSSIAN POINTS ON THE SIDE
+ x_gauss_side = half_pt_x + half_length_x * s_loc[i];
+ y_gauss_side = half_pt_y + half_length_y * s_loc[i];
+ scaling_fact_minus = (x_gauss_side - *alpha) / 2.0;
+ // COMPUTE NODES AND WEIGHTS
+ for (int j = 0; j < *N; j++) {
+ idx = j * *loc + i; // use same order as in R implementation
+ nodes_x[idx] = *alpha + scaling_fact_minus * (s_N[j] + 1.0);
+ nodes_y[idx] = y_gauss_side;
+ weights[idx] = half_length_y*scaling_fact_minus * w_loc[i] * w_N[j];
+ }
+ }
+}
+
+
+/***
+ * Function to be called from R to loop over all polygon edges,
+ * calling the above C_polygauss_side() for each
+ ***/
+
+void C_polygauss(
+ double *x, double *y, // vertex coordinates (open) of a polygon
+ double *s_M, double *w_M, // nodes & weights of Gauss-Legendre quadrature
+ double *s_N, double *w_N, // of degree M=N+1 and N, respectively
+ double *alpha, // base-line
+ int *L, int *M, int *N, // L: number of edges/vertices
+ // result: nodes and weights of length (<=) M*N per edge
+ double *nodes_x, double *nodes_y, double *weights)
+{
+ int idxTo, idxBlock;
+ double x1, y1, x2, y2;
+ for (int i = 0; i < *L; i++) {
+
+ x1 = x[i]; y1 = y[i];
+ if (i == *L-1) idxTo = 0; else idxTo = i+1;
+ x2 = x[idxTo]; y2 = y[idxTo];
+
+ // if edge is on base-line or is orthogonal to it -> skip
+ if ((x1 == *alpha && x2 == *alpha) || (y2 == y1))
+ continue;
+
+ idxBlock = i * *M * *N; // start index of nodes of edge i
+ if (x2 == x1)
+ // side is parallel to base-line -> use degree N in both dimensions
+ C_polygauss_side(&x1, &y1, &x2, &y2,
+ s_N, w_N, s_N, w_N, alpha,
+ N, N,
+ nodes_x + idxBlock, nodes_y + idxBlock, weights + idxBlock);
+ else
+ // use degrees M and N, respectively
+ C_polygauss_side(&x1, &y1, &x2, &y2,
+ s_M, w_M, s_N, w_N, alpha,
+ M, N,
+ nodes_x + idxBlock, nodes_y + idxBlock, weights + idxBlock);
+
+ }
+}
diff --git a/src/polyCub.SV.h b/src/polyCub.SV.h
new file mode 100644
index 0000000..b708843
--- /dev/null
+++ b/src/polyCub.SV.h
@@ -0,0 +1,18 @@
+/*******************************************************************************
+ * Header file of polyCub.SV.c
+ *
+ * Copyright (C) 2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+void C_polygauss(
+ double *x, double *y, // vertex coordinates (open) of a polygon
+ double *s_M, double *w_M, // nodes & weights of Gauss-Legendre quadrature
+ double *s_N, double *w_N, // of degree M=N+1 and N, respectively
+ double *alpha, // base-line
+ int *L, int *M, int *N, // L: number of edges/vertices
+ // result: nodes and weights of length (<=) M*N per edge
+ double *nodes_x, double *nodes_y, double *weights);
diff --git a/src/polyCub.iso.c b/src/polyCub.iso.c
new file mode 100644
index 0000000..e71bc28
--- /dev/null
+++ b/src/polyCub.iso.c
@@ -0,0 +1,140 @@
+/*******************************************************************************
+ * C-version of polyCub1.iso()
+ *
+ * Copyright (C) 2015,2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+/* The corresponding math is derived in Supplement B (Section 2.4) of
+ * Meyer and Held (2014): "Power-law models for infectious disease spread."
+ * The Annals of Applied Statistics, 8(3), 1612-1639.
+ * https://doi.org/10.1214/14-AOAS743SUPPB
+ */
+
+#include <math.h>
+#include <R_ext/Print.h> // Rprintf
+#include <R_ext/Memory.h> // R_alloc
+#include <R_ext/Applic.h> // Rdqags
+
+// header file defines the intrfr_fn type
+#include "polyCub.iso.h"
+
+// integrand for the edge (x0,y0) -> (x1,y1), see Equation 7
+static double lineIntegrand(
+ double t,
+ double x0, double y0, double x1, double y1,
+ intrfr_fn intrfr, double *pars)
+{
+ double num = y1*x0 - x1*y0; // numerator term
+ // point on the edge corresponding to t
+ double px = x0 + t*(x1-x0);
+ double py = y0 + t*(y1-y0);
+ double norm2 = px*px + py*py;
+ // evaluate F(R) = int_0^R r*f(r) dr at R=||(px,py)||
+ double inti = intrfr(sqrt(norm2), pars);
+ return num*inti/norm2;
+}
+
+// set of parameters for line integration (passed via the *ex argument)
+typedef struct {
+ double x0, y0, x1, y1;
+ intrfr_fn intrfr;
+ double *pars;
+} Params;
+
+// vectorized lineIntegrand for use with Rdqags
+static void myintegr_fn(double *x, int n, void *ex)
+{
+ Params *param = (Params *) ex;
+ for(int i = 0; i < n; i++) {
+ x[i] = lineIntegrand(x[i],
+ param->x0, param->y0, param->x1, param->y1,
+ param->intrfr, param->pars);
+ }
+ return;
+}
+
+// calculate line integral for one edge (x0,y0) -> (x1,y1)
+// using Gauss-Kronrod quadrature via Rdqags as declared in <R_ext/Applic.h>,
+// implemented in R/src/appl/integrate.c,
+// and used in R/src/library/stats/src/integrate.c
+static void polyiso_side(
+ double x0, double y0, double x1, double y1, // 2 vertices
+ intrfr_fn intrfr, double *pars, // F(R)
+ int subdivisions, double *epsabs, double *epsrel, // control
+ double *result, double *abserr, int *neval, int *ier) // results
+{
+ double num = y1*x0 - x1*y0; // numerator in lineIntegrand
+ // for any point p on the edge
+ if (num == 0.0) { // 'center' is part of this polygon edge
+ *result = 0.0;
+ *abserr = 0.0;
+ //*last = 0;
+ *neval = 0;
+ *ier = 0;
+ return;
+ }
+ // set of parameters for lineIntegrand
+ Params param = {x0, y0, x1, y1, intrfr, pars};
+ // prepare for Rdqags
+ double lower = 0.0;
+ double upper = 1.0;
+ int lenw = 4 * subdivisions;
+ int last; // unused
+ int *iwork = (int *) R_alloc((size_t) subdivisions, sizeof(int));
+ double *work = (double *) R_alloc((size_t) lenw, sizeof(double));
+
+ Rdqags(myintegr_fn, ¶m, &lower, &upper,
+ epsabs, epsrel,
+ result, abserr, neval, ier, // results
+ &subdivisions, &lenw, &last, iwork, work);
+
+ return;
+}
+
+// line integration along the edges of a polygon
+void polyiso(
+ double *x, double *y, // vertex coordinates (open)
+ int *L, // number of vertices
+ intrfr_fn intrfr, // F(R)
+ double *pars, // parameters for F(R)
+ double *center_x, double *center_y, // center of isotropy
+ int *subdivisions, double *epsabs, double *epsrel, // Rdqags options
+ int *stop_on_error, // !=0 means to stop at first ier > 0
+ double *value, double *abserr, int *neval) // results
+{
+ // auxiliary variables
+ double resulti, abserri;
+ int nevali, ieri;
+ double x0, y0, x1, y1;
+ int idxTo;
+ // initialize result at 0 (do += for each polygon edge);
+ *value = 0.0;
+ *abserr = 0.0;
+ *neval = 0;
+ for (int i = 0; i < *L; i++) {
+ x0 = x[i] - *center_x;
+ y0 = y[i] - *center_y;
+ idxTo = (i == *L-1) ? 0 : i+1;
+ x1 = x[idxTo] - *center_x;
+ y1 = y[idxTo] - *center_y;
+ polyiso_side(x0, y0, x1, y1,
+ intrfr, pars,
+ *subdivisions, epsabs, epsrel,
+ &resulti, &abserri, &nevali, &ieri);
+ if (ieri > 0) {
+ if (stop_on_error == 0) {
+ Rprintf("abnormal termination of integration routine (%i)\n", ieri);
+ } else {
+ error("abnormal termination of integration routine (%i)\n", ieri);
+ }
+ }
+ *value += resulti;
+ *abserr += abserri;
+ *neval += nevali;
+ }
+ return;
+}
diff --git a/src/polyCub.iso.h b/src/polyCub.iso.h
new file mode 100644
index 0000000..5e83953
--- /dev/null
+++ b/src/polyCub.iso.h
@@ -0,0 +1,21 @@
+/*******************************************************************************
+ * Header file of polyCub.iso.c
+ *
+ * Copyright (C) 2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+typedef double (*intrfr_fn) (double, double*);
+
+void polyiso(
+ double *x, double *y, // vertex coordinates (open)
+ int *L, // number of vertices
+ intrfr_fn intrfr, // F(R)
+ double *pars, // parameters for F(R)
+ double *center_x, double *center_y, // center of isotropy
+ int *subdivisions, double *epsabs, double *epsrel, // Rdqags options
+ int *stop_on_error, // !=0 means to stop at first ier > 0
+ double *value, double *abserr, int *neval); // results
diff --git a/tests/test-all.R b/tests/test-all.R
new file mode 100644
index 0000000..47bdd3a
--- /dev/null
+++ b/tests/test-all.R
@@ -0,0 +1,3 @@
+if (require("testthat") && packageVersion("testthat") >= "0.9") {
+ test_check("polyCub")
+}
diff --git a/tests/testthat/polyiso_powerlaw.c b/tests/testthat/polyiso_powerlaw.c
new file mode 100644
index 0000000..f1cf8b3
--- /dev/null
+++ b/tests/testthat/polyiso_powerlaw.c
@@ -0,0 +1,46 @@
+/*******************************************************************************
+ * Example of using the C-routine "polyCub_iso", see also test-polyiso.R
+ *
+ * Copyright (C) 2015,2017 Sebastian Meyer
+ *
+ * This file is part of the R package "polyCub",
+ * free software under the terms of the GNU General Public License, version 2,
+ * a copy of which is available at http://www.R-project.org/Licenses/.
+ ******************************************************************************/
+
+#include <math.h>
+#include <polyCubAPI.h>
+
+// F(R) example
+static double intrfr_powerlaw(double R, double *logpars)
+{
+ double sigma = exp(logpars[0]);
+ double d = exp(logpars[1]);
+ if (d == 1.0) {
+ return R - sigma * log(R/sigma + 1);
+ } else if (d == 2.0) {
+ return log(R/sigma + 1) - R/(R+sigma);
+ } else {
+ return (R*pow(R+sigma,1-d) - (pow(R+sigma,2-d) - pow(sigma,2-d))/(2-d)) / (1-d);
+ }
+}
+
+// function to be called from R
+void C_polyiso_powerlaw(
+ double *x, double *y, // vertex coordinates (open)
+ int *L, // number of vertices
+ //intrfr_fn intrfr, // F(R)
+ double *pars, // parameters for F(R)
+ double *center_x, double *center_y, // center of isotropy
+ int *subdivisions, double *epsabs, double *epsrel, // Rdqags options
+ int *stop_on_error, // !=0 means to stop at first ier > 0
+ double *value, double *abserr, int *neval) // results
+{
+ polyCub_iso(x, y, L,
+ intrfr_powerlaw,
+ pars, center_x, center_y,
+ subdivisions, epsabs, epsrel,
+ stop_on_error,
+ value, abserr, neval);
+ return;
+}
diff --git a/tests/testthat/test-NWGL.R b/tests/testthat/test-NWGL.R
new file mode 100644
index 0000000..70d0ab7
--- /dev/null
+++ b/tests/testthat/test-NWGL.R
@@ -0,0 +1,9 @@
+context("Validation of cached Gauss-Legendre nodes/weights")
+
+if (requireNamespace("statmod")) {
+ test_that("statmod::gauss.quad() still gives the same result", {
+ new.NWGL <- lapply(seq_len(61L), function (n)
+ unname(statmod::gauss.quad(n = n, kind = "legendre")))
+ expect_equal(new.NWGL, .NWGL, check.attributes = FALSE)
+ })
+}
diff --git a/tests/testthat/test-polyCub.R b/tests/testthat/test-polyCub.R
new file mode 100644
index 0000000..c25c1df
--- /dev/null
+++ b/tests/testthat/test-polyCub.R
@@ -0,0 +1,59 @@
+context("Correctness of cubature methods")
+
+### set up test case
+
+## bivariate, isotropic Gaussian density
+f <- function (s, mean, sd)
+ dnorm(s[,1], mean=mean[1], sd=sd) * dnorm(s[,2], mean=mean[2], sd=sd)
+
+## circular domain represented by a polygon
+r <- 5
+center <- c(3,2)
+npoly <- 128
+disc.owin <- spatstat::disc(radius=r, centre=center, npoly=npoly)
+
+## parameters for f
+m <- c(1,1)
+sd <- 3
+
+## exact value of the integral over the _polygonal_ circle
+intExact <- 0.65844436
+if (requireNamespace("mvtnorm") && gpclibPermit()) {
+ ## run this conditionally since gpclib might not be available on all
+ ## platforms (as pointed out by Uwe Ligges, 2014-04-20)
+ test_that("polyCub.exact.Gauss returns validated result", {
+ int <- polyCub.exact.Gauss(disc.owin, mean=m, Sigma=sd^2*diag(2))
+ expect_equal(int, intExact, tolerance=1e-8, check.attributes=FALSE)
+ })
+}
+
+
+### perform the tests (check against each other)
+
+test_that("polyCub.exact.Gauss and circleCub.Gauss give similar results", {
+ ## exact value of the integral over the _real_ circle
+ intExact_circle <- circleCub.Gauss(center=center, r=r, mean=m, sd=sd)
+
+ ## how well this fits with the exact integral over a polyonal approximation
+ ## of the circle depends of course on 'npoly'
+ expect_equal(intExact, intExact_circle,
+ tolerance=0.001, check.attributes=FALSE)
+})
+
+test_that("midpoint-cubature is correct", {
+ int <- polyCub.midpoint(disc.owin, f, mean=m, sd=sd, dimyx=500)
+ expect_equal(int, intExact, tolerance=0.001, check.attributes=FALSE)
+})
+
+test_that("SV-cubature is correct", {
+ intC <- polyCub.SV(disc.owin, f, mean=m, sd=sd, nGQ=3, engine="C")
+ intR <- polyCub.SV(disc.owin, f, mean=m, sd=sd, nGQ=3, engine="R")
+ expect_equal(intC, intR)
+ expect_equal(intC, intExact, tolerance=0.0001, check.attributes=FALSE)
+})
+
+test_that("isotropic cubature is correct", {
+ ## using a numerical approximation of intrfr
+ int0 <- polyCub.iso(disc.owin, f, mean=m, sd=sd, center=m)
+ expect_equal(int0, intExact, check.attributes=FALSE)
+})
diff --git a/tests/testthat/test-polyiso.R b/tests/testthat/test-polyiso.R
new file mode 100644
index 0000000..5695df0
--- /dev/null
+++ b/tests/testthat/test-polyiso.R
@@ -0,0 +1,103 @@
+context("polyCub_iso C-routine (API)")
+
+## CAVE (as of R-3.4.0 with testthat 1.0.2):
+## During R CMD check, tools:::.runPackageTests() sets R_TESTS=startup.Rs,
+## a file which is created in the parent directory "tests", see
+## file.path(R.home("share"), "R", "tests-startup.R")
+## for its contents. However, testthat tests are run with the working directory
+## set to here, so auxiliary R sessions initiated here would fail when trying
+## to source() the R_TESTS file on startup, see the system Rprofile file
+## file.path(R.home("library"), "base", "R", "Rprofile")
+## for what happens. Solution: unset R_TESTS (or set to "") for sub-R processes.
+
+## function to call an R CMD with environment variables
+## 'env' specified as a named character vector
+Rcmd <- function (args, env = character(), ...) {
+ stopifnot(is.vector(env, mode = "character"),
+ !is.null(names(env)))
+ if (.Platform$OS.type == "windows") {
+ if (length(env)) {
+ ## the 'env' argument of system2() is not supported on Windows
+ setenv <- function (envs) {
+ old <- Sys.getenv(names(envs), unset = NA, names = TRUE)
+ set <- !is.na(envs)
+ if (any(set)) do.call(Sys.setenv, as.list(envs[set]))
+ if (any(!set)) Sys.unsetenv(names(envs)[!set])
+ invisible(old)
+ }
+ oldenv <- setenv(env)
+ on.exit(setenv(oldenv))
+ }
+ system2(command = file.path(R.home("bin"), "Rcmd.exe"),
+ args = args, ...)
+ } else {
+ system2(command = file.path(R.home("bin"), "R"),
+ args = c("CMD", args),
+ env = paste(names(env), env, sep = "="), ...)
+ }
+}
+
+message("compiling polyiso_powerlaw.c using R CMD SHLIB")
+shlib_error <- Rcmd(args = c("SHLIB", "--clean", "polyiso_powerlaw.c"),
+ env = c("PKG_CPPFLAGS" = paste0("-I", system.file("include", package="polyCub")),
+ "R_TESTS" = ""))
+if (shlib_error)
+ skip("failed to build the shared object/DLL for the polyCub_iso example")
+
+## load shared object/DLL
+myDLL <- paste0("polyiso_powerlaw", .Platform$dynlib.ext)
+loadNamespace("polyCub")
+dyn.load(myDLL)
+
+## R function calling C_polyiso_powerlaw
+polyiso_powerlaw <- function (xypoly, logpars, center,
+ subdivisions = 100L, rel.tol = .Machine$double.eps^0.25,
+ abs.tol = rel.tol, stop.on.error = TRUE)
+{
+ .C("C_polyiso_powerlaw",
+ as.double(xypoly$x), as.double(xypoly$y), as.integer(length(xypoly$x)),
+ as.double(logpars),
+ as.double(center[1L]), as.double(center[2L]),
+ as.integer(subdivisions), as.double(abs.tol), as.double(rel.tol),
+ as.integer(stop.on.error),
+ value = double(1L), abserr = double(1L), neval = integer(1L))[c("value", "abserr", "neval")]
+}
+
+
+## example polygon and function parameters
+set.seed(1)
+xy <- list(x = stats::rnorm(10), y = stats::rnorm(10))
+hidx <- grDevices::chull(xy)
+xypoly <- lapply(xy, "[", rev(hidx)) # anticlockwise vertex order
+logpars <- log(c(0.5, 1))
+
+(res <- polyiso_powerlaw(xypoly, logpars, center = c(0,0)))
+
+
+## compare with R implementation
+intrfr.powerlaw <- function (R, logpars)
+{
+ sigma <- exp(logpars[[1L]])
+ d <- exp(logpars[[2L]])
+ if (d == 1) {
+ R - sigma * log(R/sigma + 1)
+ } else if (d == 2) {
+ log(R/sigma + 1) - R/(R+sigma)
+ } else {
+ (R*(R+sigma)^(1-d) - ((R+sigma)^(2-d) - sigma^(2-d))/(2-d)) / (1-d)
+ }
+}
+(orig <- polyCub:::polyCub1.iso(xypoly, intrfr.powerlaw, logpars, center = c(0,0)))
+
+
+test_that("C and R implementations give equal results", {
+ expect_equal(res$value, orig[1])
+ expect_equal(res$abserr, orig[2])
+})
+
+## microbenchmark::microbenchmark(
+## polyCub:::polyCub1.iso(xypoly, intrfr.powerlaw, logpars, center = c(0,0)), # 250 mus
+## polyiso_powerlaw(xypoly, logpars, center = c(0,0)), times = 1000) # 50 mus
+
+dyn.unload(myDLL)
+file.remove(myDLL)
diff --git a/tests/testthat/test-regression.R b/tests/testthat/test-regression.R
new file mode 100644
index 0000000..70f68a2
--- /dev/null
+++ b/tests/testthat/test-regression.R
@@ -0,0 +1,12 @@
+context("Regression tests")
+
+test_that("isotropic cubature can handle control list for integrate()", {
+ data("letterR", package="spatstat", envir=environment())
+ f <- function (s) (rowSums(s^2)+1)^-2
+ ## previosly, passing control arguments did not work
+ int1 <- polyCub.iso(letterR, f, center=c(0,0), control=list(rel.tol=1e-3))
+ int2 <- polyCub.iso(letterR, f, center=c(0,0), control=list(rel.tol=1e-8))
+ ## results are almost identical
+ expect_equal(int1, int2, tolerance=1e-3)
+ expect_false(identical(int1, int2))
+})
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-polycub.git
More information about the debian-science-commits
mailing list