cfunc.new: Define and evaluate a contour function
In gensphere: Generalized Spherical Distributions

Description Usage Arguments Details Value Examples

View source: R/gensphere.R

The directional part of a generalized spherical distribution is defined by a contour function, cfunc for short. These functions are used to define a contour function and then evaluate it.

cfunc.new(d)
cfunc.add.term(cfunc, type, k)
cfunc.finish(cfunc, nsubdiv = 2, maxEvals=100000, norm.const.method="integrate",...)
cfunc.eval(cfunc, x)

`d`	dimension of the space
`cfunc`	an object of class "gensphere.contour"
`type`	string describing what type of term to add to the contour
`k`	a vector of constants used to specify a term
`x`	a (d x n) matrix, with columns x[,i] being points in R^d
`nsubdiv`	number of dyadic subdivisions of the sphere, controls the refinement of the tessellation. Typical values are 2 or 3.
`maxEvals`	maximum number of evaluations of the integrand function, see details section below
`norm.const.method`	method used to compute the norming constant. Can be "integrate" (the default, in which case the contour function is numerically integrated to get the norming constant), or "simplex.area" (in which case no integration is done; the surface area of the contour is approximated by the surface area of the tesselation). This later choice is for complex surface or higher dimensional surfaces where the numerical integration may fail or take a very long time. This allows simulation in cases where the numerical integration is not feasible.
`...`	optional arguments to pass to integration routine, see details section below

A contour function describes the directional behavior of a generalized spherical distribution. In this package, a contour function is built by calling three functions: cfunc.new to start the definition of a d-dimensional contour, cfunc.add.term to add a new term to a contour (may be called more than once), and cfunc.finish to complete the defintion of a contour.

When adding a term, type is one of the literal strings "constant", "elliptical", "proj.normal", "lp.norm", "gen.lp.norm", or "cone". The vector k contains the constants necessary to specify the desired shape. k[1]=the first element of k is always a scale, it allows one to expand or contract a shape. The remaining elements of k depend on the value of 'type'.:

"constant": k is a single number; k[1]=radius of the sphere
"elliptical": k is a vector of length d^2+1; k[1]=scale of ellipse, k[2:(d^2+1)] specify the symmetric positive definite matrix B with is used to compute the elliptical contour
"proj.normal": k is a vector of length d+2; k[1]=scale of the bump, mu=k[2:(d+1)] is the vector pointing in the direction where the normal bump is centered, k[d+2]=standard deviation of the isotropic normal bump
"lp.norm": k is a vector of length 2; k[1]=scale and k[2]=p, the power used in the l_p norm
"gen.lp.norm": k is vector of length 2+m*d for some positive integer m; k[1]=scale, k[2]=p, the power used in the l-p norm, k[3:(2+m*d)] spcifies a matrix A that is used to compute || A x ||_p
"cone": k is a vector of length d+2, k[1]=scale, mu=k[2:(d+1)]= the center of the cone, k[d+2]=base of the cone

Note that cfunc.finish does a lot of calculation, and may take a while, especially in dimension d > 2. The most time consuming part is numerically integrating over the contour, a (d-1) dimensional surface in d-dimensional space and in tesselating the contour in a way that focuses on the bulges in the contour from cones and normal bumps. The integration is required to calculate the norming constant needed to compute the density. This integration is performed by using function adaptIntegrateSphereTri in SphericalCubature and is numerically challenging. In dimension d > 3 or if nsubdiv > 4, users may need to adjust the arguments maxEvals and ... The default value maxEvals=100000 workw in most 3 dim. problems, and it takes a few seconds to calculate. (For an idea of the size and time required, a d=4 dim. case used maxEvals=1e+7 and took around 5 minutes. A d=5 dim. case used maxEvals=1e+8, used 160167 simplices and took over 2 days.) Note that this calculation is only done once; calculating densities and simulating is still fast in higher dimensions. It may be useful to save a complicated/large contour object so that it can be reused across R sessions via save(cfunc) and load(cfunc).

Note: the first time cfunc.finish is called, a warning message about "no degenerate regions are returned" is printed by the package geometry. I do not know how to turn that off; so just ignore it.

cfunc.eval is used to evaluate a completed contour function.

cfunc.new and cfunc.add.term return a list that is an incomplete definition of a contour function. cfunc.finish completes the definition and returns an S3 object of class "gensphere.contour" with fields:

`d`	dimension
`m`	number of terms in the contour function, i.e. the number of times `cfunc.add.term` was called
`term`	a vector length m of type list, with each list describing a term in the contour function
`norm.const`	norming constant
`functionEvaluations`	number of times the integrand (contour) function is evaluated by `adaptIntegrateSphereTri` when computing norm.const
`tessellation`	an object of type "mvmesh" that gives a geometrical description of the contour. It is used to plot the contour and to simulate from the contour
`tessellation.weights`	weights used in simulation; surface area of the simplices in the tessellation
`simplex.count`	vector of length 3 giving the number of simplices used at the end of three internal stages of construction: after initial subdivision, after refining the sphere based on cones and bumps, and final count have adaptive integration routine partitions the sphere
`norm.const.method`	value of input argument norm.const.method

cfunc.eval returns a vector of length n=nrow(x); y[i] = cfunc(x[,i]) = value of the contour function at point x[,i].

The plots below show the three contours defined in the examples below.

# 2-dim diamond
cfunc1 <- cfunc.new(d=2)
cfunc1 <- cfunc.add.term( cfunc1,"lp.norm",k=c(1,1))
cfunc1 <- cfunc.finish( cfunc1 )
cfunc1
cfunc.eval( cfunc1, c(sqrt(2)/2, sqrt(2)/2) )
if( interactive()) { plot( cfunc1, col='red', lwd=3, main="diamond contour") }

# 2-dim blob
cfunc2 <- cfunc.new(d=2)
cfunc2 <- cfunc.add.term( cfunc2,"constant",k=1)
cfunc2 <- cfunc.add.term( cfunc2,"proj.normal",k=c( 1, sqrt(2)/2, sqrt(2)/2, 0.1) )
cfunc2 <- cfunc.add.term( cfunc2,"proj.normal",k=c( 1, -1,0, 0.1) )
cfunc2 <- cfunc.finish( cfunc2, nsubdiv=4 )
if(interactive()) {
  plot( cfunc2, col='green', lwd=3, main="contour with two bumps")

  # 3-dim ball with one spike
  cfunc3 <- cfunc.new( d=3 )
  cfunc3 <- cfunc.add.term( cfunc3, "elliptical",k=c( 1,  1,0,0, 0,1,0, 0,0,1 ) )
  cfunc3 <- cfunc.add.term( cfunc3, "proj.normal",k=c( 1, 1,0,0, .25 ) ) 
  cfunc3 <- cfunc.finish( cfunc3, nsubdiv=3 ) # takes ~20 seconds, get warnings
  plot( cfunc3, show.faces=TRUE, col='blue')
  nS <- dim(cfunc3$tessellation$S)[3]
  title3d( paste("ball with bump with",nS,"simplices"))
}