Nothing
# SimplicialCubature - an R package to integrate over m-dimensional simplices in n-dimensions.
# Is is required that n >= 2, while the simplices can be of dimension m, 1 <= m <= n.
#
# Written by John Nolan, American University, Aug-Sep 2014. Modified in Spring 2015.
# e-mail jpnolan@american.edu
#
# This research was supported by an agreement with Cornell University, Operations
# Research & Information Engineering, under contract W911NF-12-1-0385 from the Army
# Research Development and Engineering Command.
#
# There are two main functions in this package:
# 1. Function adaptIntegrateSimplex
# This function adaptively integrates a function f(x) over a simplex in
# n-dimensional space, n >= 2.
# (a) If the the simplicies are n-dimensional, then the integration
# is performed by function adsimp(...), an R
# translation of Alan Genz's matlab function adsimp, which
# adaptively integrates an arbitrary function over an n-dim. simplex.
# See comments in function adsimp for references on the method. The only
# change to that method is that a new argument partitionInfo is added;
# if that is TRUE, adsimp will return information about
# the subdivision of the simplices.
# (b) If the simplex is of dimension 1 <= m < n, then the problem is transformed
# into an integral over m-dimensional space, and evaluated using adsimp(...)
#
# In both cases, the adaptIntegrateSimplex function allows the integrand function
# to have optional arguments by using "..." ; this is achieved by defining
# a new integrand function locally that incorporates the extra arguments,
# leaving adsimp unaware of that construct.
#
# 2. Function integrateSimplexPolynomial
# R implementations of the Lasserre-Avrachenkov and Grundmann-Moller
# methods for exactly integrating a polynomial over a simplex.
# The polynomial p(x) must be defined by using function polyDefine.
# It is specified as a list with fields:
# p$coef[1:m] vector of doubles, the coefficients of the terms
# p$k[1:m,1:n] matrix of integers
# m and n are given implicitly in the sizes of these arrays
# Precisely, p(x[1],...,x[n])=sum (coef[i] * x[1]^k[i,1] * ... * x[n]^k[i,n]),
# where the sum is over i=1,...,m.
#
# The functions polyDefine, polyEval, and polyPrint are used to define,
# evaluate, and print a polynomial in coded form.
#
# The functions CanonicalSimplex, UnitSimplexV, SimplexVolume, and
# JacobianS2Canonical are utility routines to define the canonical
# simplex, the unit simples, volume of an n-dim. simplex in R^n,
# and the Jacobian of the transformation from an m-dim. simplex S2
# to the m-dim. canonical simplex.
#
# All other functions are internal and are used at your own risk.
#
################################################################################
adaptIntegrateSimplex <- function( f, S, fDim=1L, maxEvals=10000L, absError=0.0,
tol=1.0e-5, integRule=3L, partitionInfo=FALSE, ... ) {
# Adaptively integrate f(x) over one or more simplices in R^n, where
# S[1:n,1:m,1:nS] may be a single simplex or a list of simplices.
# The region of integration can be an m-dimensional simplex, 1 <= m <= n.
#
# When m=n, this is a wrapper function for Alan Genz's adsimp program to adaptively
# integrate a function f(x) over an n-dimensional simplex S. f(x) can be vector valued
# and S can be an array of simplices. We change variable names to parallel
# the variable names in CRAN packages cubature and SphericalCubature.
# When m < n, the problem is transformed to integration in m-space and then
# adsimp is used to evaluate that m-dim. integral.
#
# Input variables:
# f integrand function, returns a vector of length fDim
# S either a single n x (m+1) matrix, with S[,i] being the i-th vertex
# of the simplex, or an n x (m+1) x nS array of simplices, with S[,,i]
# being the i-th simplex in the first form
# fDim integer that gives the dimension of (the range of) function f
# maxEvals integer that gives an upper limit for how many evaluations of
# of f are allowed. If this value is exceeded (or about to be exceeded),
# the return variable will contain a non-zero value.
# absError requested absolute error in the integral
# tol requested relative error in the integral
# integRule integer in the range 1:4 specifiying the integration rule used
# (see function adsimp for more info)
# partitionInfo TRUE/FALSE value indicating whether partition information
# is be be returned. This takes more memory (and time if m < n), but
# the results may be of use in some applications. (see function adsimp for more info)
# ... optional arguments to f
#
# Returns a list with multiple fields:
# integral estimated values of integrals, vector if fDim > 1
# estAbsError estimated absolute error
# functionEvaluations count of number of function evaluations
# returnCode integer status: returnCode=0 is a successful return
# message a text string describing what returnCode means
# if partitionInfo=TRUE, other information about the final partition/subdivision
# of the region of integration is returned, see function adsimp for more details.
if( !is.function(f) ) { stop("first argument is not a function" ) }
if( is.matrix(S) ) { S <- array( S, dim=c(nrow(S),ncol(S),1)) }
if( !is.array(S) | (length(dim(S)) != 3) ) stop("S must be a single simplex or an array of simplices")
dimS <- dim(S); n <- dimS[1]; m <- dimS[2]-1; nS <- dimS[3]
if (n < 2) stop("dimension of the space is less than 2")
if( m > n ) stop("invalid simplices: m > n")
if( n == m ) { # region is a set of n-dim simplices inside R^n
# define a new function to evaluate f using optional arguments in ...
fNew <- function( x ) { f( x, ... ) }
# call adsimp to do integration
a <- adsimp( ND=n, VRTS=S, NF=fDim, F=fNew, MXFS=maxEvals, EA=absError,
ER=tol, KEY=integRule, partitionInfo )
} else {
# integration region is a set of m-dim. simplices inside R^n with m < n.
# transform the problem to an m-dimensional one: map simplex S[,,i]
# to a translation of the canonical simplex by (2*i,0,0,...,0).
# This translation allows the function transformedF to know
# which region of integration a point is in simply by looking at its first
# coordinate. This allows us to incorporate the change of variable
# coefficient, JacobianCoef[i] for simplex S[,,i], into the transformed
# integrand function.
S0 <- CanonicalSimplex( m )
transformedS <- array( 0.0, dim=c(m,m+1,nS) )
JacobianCoef <- rep(0.0,nS)
for (i in 1:nS) {
JacobianCoef[i] <- JacobianS2Canonical( S[,,i] )
S0[1, ] <- S0[1, ] + 2
transformedS[,,i] <- S0
}
# integrate transformed f over transformed simplices
if (m==1) { # evaluate a line integral
# define the transformed function which evaluates the original function f(x) at
# multiple points in R^n, returning multiple values.
transformedF <- function( u ) {
# in 1-dim case, u is a vector of distinct points in R.
# since the function adaptIntegrateVectorFunc does one interval at a
# time, all points in u must be in the same original parition interval.
nu <- length(u)
y <- matrix(0.0,nrow=fDim,ncol=nu)
for (i in 1:nu) {
b <- original.coordinates( u[i], S )
y[,i] <- JacobianCoef[b$i] * f( b$x, ... )
}
return(y) }
a <- adaptIntegrateVectorFunc( transformedS, fDim, transformedF, maxEvals, absError, tol, partitionInfo, ... )
} else { # evaluate m-dim. integral where 2 <= m < n
# define the transformed function which evaluates the original function f(x) at a
# a single point in R^n, returning a single number
transformedF <- function( u ) {
b <- original.coordinates( u, S )
y <- JacobianCoef[b$i] * f( b$x, ... )
return(y) }
a <- adsimp( ND=m, VRTS=transformedS, NF=fDim, F=transformedF, MXFS=maxEvals, EA=absError,
ER=tol, KEY=integRule, partitionInfo )
}
}
if (m==1) { # in univariate case, just return result of adaptIntegrateVectorFunc
result <- a
} else { # rename results and return to caller
if (partitionInfo) {
result <- list(integral=a$VL, estAbsError=a$AE, functionEvaluations=a$NV,
returnCode=a$FL, subsimplices=a$VRTS, subsimplicesIntegral=a$VLS,
subsimplicesAbsError=a$AES, subsimplicesVolume=a$VOL )
} else {
result <- list(integral=a$VL, estAbsError=a$AE, functionEvaluations=a$NV,
returnCode=a$FL )
}
}
if(result$returnCode == 0) {
if (partitionInfo & ( m < n )) {
# the current partition info applies to transformed simplices,
# translate those simplices back to the original coordinates.
n3 <- dim(result$subsimplices)[3]
newSubsimplices <- array( 0.0, dim=c(n,m+1,n3) )
for (i in 1:n3) {
for (j in 1:(m+1)) {
b <- original.coordinates( result$subsimplices[,j,i], S )
newSubsimplices[,j,i] <- b$x
}
result$subsimplicesVolume[b$i] <- JacobianCoef[b$i]*result$subsimplicesVolume[b$i]
}
result$subsimplices <- newSubsimplices
}
}
# convert returnCode to a message string
result$message <- adsimp.return.message( result$returnCode )
return(result) }
################################################################################
adsimp.return.message <- function( rcode ) {
# return a text string explaining the return code recode from adsimp
if (rcode == 0) msg <- "OK"
if (rcode == 1) msg <- "error: maxEvals exceeded - too many function evaluations"
if (rcode == 2) msg <- "error: integRule < 0 or integRule > 4"
if (rcode == 3) msg <- "error: dimension n of the space < 2"
if (rcode == 4) msg <- "error: fDim = dimension of f < 1"
if (rcode == 5) msg <- "error: absError < 0 and tol < 0"
if ((rcode < 0) | (rcode > 5)) msg <- paste("error: unknown return code = ",rcode )
return(msg) }
################################################################################
adaptIntegrateVectorFunc <- function( intervals, fDim, f, maxEvals, absError, tol, partitionInfo=FALSE,... ){
# Evaluate integral of a vector valued function f(x)=(f[1](x),...,f[fDim](x)) of a real
# variable x (not a vector) over a list of intervals. The built-in R function integrate
# uses quadpack function dqags when lower and upper are both finite, and that appears
# to evaluate f at 21 points each time.
dqags.num.points <- 21L
if( !is.function(f) ) { stop("first argument is not a function" ) }
# define function to return a single coordinate of vector valued f(x)
# This is inefficient, but
newF <- function( x, which.coord ) {
# make sure right number of values requested
if(length(x)!=dqags.num.points) { warning(paste("adaptIntegrateVectorFunc: length(x) != ",dqags.num.points)) }
x <- as.vector(x)
y <- rep(0.0,length(x))
for (j in 1:length(x)) {
y[j] <- f(x[j],...)[which.coord]
}
return(y) }
n.intervals <- dim(intervals)[3]
max.subdiv <- as.integer( maxEvals/dqags.num.points )
abs.tol <- absError/n.intervals
y <- rep(0.0,fDim)
if (partitionInfo) {
y.partition <- matrix( 0.0, nrow=fDim, ncol=n.intervals)
est.abs.error.partition <- y.partition
}
est.abs.error <- rep(0.0,fDim)
actual.subdiv <- rep( 0L, fDim )
# loop over coordinates of vector valued f
for (i in 1:fDim) {
# loop over intervals
for (j in 1:n.intervals) {
tmp <- integrate( newF, lower=intervals[1,1,j], upper=intervals[1,2,j], rel.tol=tol,
abs.tol=abs.tol, subdivisions=max.subdiv, which.coord=i )
if (tmp$message != "OK") { break }
y[i] <- y[i] + tmp$value
est.abs.error[i] <- est.abs.error[i] + tmp$abs.error
actual.subdiv[i] <- actual.subdiv[i] + tmp$subdivisions
if (partitionInfo) {
y.partition[i,j] <- y.partition[i,j] + tmp$value
est.abs.error.partition[i,j] <- est.abs.error.partition[i,j]+tmp$abs.error
}
}
}
result <- list(integral=y, estAbsError=est.abs.error,
functionEvaluations=dqags.num.points*sum(actual.subdiv),
returnCode=ifelse(tmp$message=="OK",0L,-1L),
message=tmp$message )
if(partitionInfo) {
result$subsimplices <- intervals
result$subsimplicesIntegral <- y.partition
result$subsimplicesAbsError <- est.abs.error.partition
result$subsimplicesVolume <- intervals[1,2,]-intervals[1,1,]
}
return(result)}
################################################################################
original.coordinates <- function( u, S ) {
# when m < n, this function is used to transform from the m-dimensional
# new coordinates to original n-dimensional coordinates.
# u should be a single point in R^m, e.g. a vector of length m.
# The basic idea is barycentric coordinates, after correcting for the
# deliberate shift of the i-th simplex by (2*i,0,0,...,0).
# x <- u[1]*S[,1,i] + u[2]*S[,2,i]+...+u[m]*S[,m,i]+ (1-sum(u))*S[,m+1,i]
#
# This function returns a list with two fields:
# $i for the index of the simplex
# $x for the point in R^n corresponding to u in R^m
i <- as.integer(u[1]/2.0) # figures out which simplex u is in
u[1] <- u[1]- 2*i # subtract the shift to get barycentric coordinates
v <- c( u, 1-sum(u) )
return( list( i=i, x= S[,,i] %*% v ) ) }
################################################################################
integrateSimplexPolynomial <- function( p, S, method="GM" ) {
#
# Exactly integrate a polynomial p(x) over a simplex.
#
# Input variables:
# p a polynomial, which should be defined through function polyDefine.
# S is either an n x (m+1) matrix with S[,j] the j-th vertex of the simplex.
# Alternatively, S can be an (n x (m+1) x k) array, with S[,,i] giving the
# i-th simplex as in the matrix case.
# method is a string, one of:
# "GM" for Grundmann-Moller method, works for 1 < m <= n.
# "LA" for the Lasserre-Avrachenkov method, only works for m=n.
#
# The formulas used here are exact for polynomials, not approximate, but
# use floating point arithmetic to evaluate the polynomial. So there will be
# round-off due to this floating point evaluation. If the coefficients of the
# polynomial are all rational numbers, and the vertices are all rational numbers,
# it is possible to evaluate the polynomials exactly as rational numbers, and
# the integral can be evaluated exactly. We do not do this here, though it should
# be possible using CRAN packages gmp or Rmpfr. Exact rational arithmetic is done
# in the C program LattE, which is available from the website:
# https://www.math.ucdavis.edu/~latte
#
# If method="GM", then the Grundmann-Moler algorithm of order ceiling( (degree(p)-1) / 2 )
# is used. The GM method is exact for this order on a polynomial. See
# function grnmol( ) for details of the method and references.
#
# If method="LA", this function splits the polynomial into a list of q-homogeneous
# polynomials and calls function LasserreAvrachenkov( ) to exactly integrate a
# q-homogeneous polynomial over each simplex S[,,i].
# if S is a single simplex, make it an array of simplices of length 1
if (is.matrix(S)) { S <- array( S, dim=c(dim(S),1) ) }
# check input parameters
if( !is.array(S) ) stop("S must be a single simplex (a matrix) or a list of matrices (a 3-dim. array)")
# check input parameters
dimS <- dim(S)
if ( length(dimS) != 3) stop( "S must be a single simplex (a matrix) or a list of matrices (a 3-dim. array)" )
stopifnot( is.vector(p$coef), is.matrix(p$k), dim(p$k)[2]==dimS[1], method %in% c("GM","LA") )
n <- dimS[1]; m <- dimS[2]-1; k <- dimS[3]
if( m > n) stop("m must be <= n")
if( n < 1 ) stop("n=nrow(S) must be bigger >= 1")
if( m < 1 ) stop("m=ncol(S)-1 must be bigger than 1")
# compute the degrees of each of the terms in p(.)
deg <- rowSums(p$k)
# Grundmann-Moller algorithm
if (method=="GM") {
gm.order <- ceiling( (max(deg)-1) / 2 ) # order of Grundmann-Moller method
if (gm.order < 1) { gm.order <- 1 }
integral <- 0; functionEvaluations <- 0
if (m == n) {
# polynomial interface to grnmol
fnew <- function( x ) { evalPoly( x, p ) }
for (i in 1:k) {
a <- grnmol( fnew, S[ , , i], gm.order )
integral <- integral + a$Q[gm.order+1]
functionEvaluations <- functionEvaluations + a$nv
}
} else { # 1 <= m < n, map to an integral over m-dim. canonical simplex
# integration region is a set of m-dim. simplices inside R^n with m < n.
# transform the problem to an m-dimensional one: map simplex S[,,i]
# to a translation of the canonical simplex by (2*i,0,0,...,0).
# This translation allows the function transformedF to know
# which region of integration a point is in simply by looking at its first
# coordinate. This allows us to incorporate the change of variable
# coefficient, JacobianCoef[i] for simplex S[,,i], into the transformed
# integrand function.
S0 <- CanonicalSimplex( m )
transformedS <- array( 0.0, dim=c(m,m+1,k) )
JacobianCoef <- rep(0.0,k)
for (i in 1:k) {
JacobianCoef[i] <- JacobianS2Canonical( S[,,i] )
S0[1, ] <- S0[1, ] + 2
transformedS[,,i] <- S0
}
# define the transformed function which evaluates the original function f(x) at a
# a single point in R^n, returning a single number
transformedF <- function( u ) {
b <- original.coordinates( u, S )
y <- JacobianCoef[b$i] * evalPoly( b$x, p )
return(y) }
for (i in 1:k) {
curS <- transformedS[ , , i]
if (!is.matrix(curS)) { curS <- matrix( curS, nrow=1 ) }
a <- grnmol( transformedF, curS, gm.order )
integral <- integral + a$Q[gm.order+1]
functionEvaluations <- functionEvaluations + a$nv
}
}
}
# Lasserre-Avrachenkov method
if (method=="LA") {
# in principle, L-A method can be made to work when m < n, but it is
# messy. Since it requires that you separate the terms of the polynomial
# into homogeneous terms, you have to express the transformed polynomial
# in m-dimensions in polynomial form. This requires some effort; we do
# not see it worth the effort now.
if (m != n) stop("Lasserre-Avranchenkov method only works with m=n")
if (n == 1) stop("Lasserre-Avranchenkov method only works for n > 1")
unique.deg <- unique(deg)
# Loop through terms that are homog. of each degree q and all simplices.
# Use a trick of setting some coef. to 0 when the corresponding term is
# not of degree q.
integral <- 0.0; functionEvaluations <- 0
for ( q in unique.deg ) {
useTerm <- ( deg == q ) # useTerm[i] = 1 if term i is q-homog., 0 otherwise
for (i in 1:k) {
tmp <- LasserreAvrachenkov( q, p, useTerm, S[ , , i] )
integral <- integral + tmp$integral
functionEvaluations <- functionEvaluations + tmp$functionEvaluations
}
}
}
return( list(integral=integral,functionEvaluations=functionEvaluations )) }
################################################################################
definePoly <- function( coef, k ) {
# Define a polynomial in n variables with coefficients coef[1:m] and powers k[1:m,1:n]
# this is essentially an S3 method
# check sizes
stopifnot( is.vector(coef), is.numeric(coef), is.matrix(k), length(coef) == nrow(k) )
# as.integer drops the matrix dimensions
return( list( coef=coef, k=matrix( as.integer(k), nrow=nrow(k) ) ) ) }
################################################################################
evalPoly <- function( x, p, useTerm=rep(TRUE,length(p$coef)) ) {
# Evaluate the polynomial y[j] = p(x[,j]) at j=1,...,nx
#
# The vector useTerm indicates which terms in polynomial p are
# used in this evaluation: if useTerm[i]=TRUE, include the i-th term
# in the function evaluation. This is mainly for efficiency reasons:
# function integrateSimplexPolynomial splits the polynomial into terms of
# different degrees, and this vector of indicators allows us to use only
# the terms required. This avoids repeated copying of terms, and will
# be simpler if these functions are coded in C.
if (is.vector(x)) { x <- matrix(x,ncol=1) }
stopifnot( nrow(x) == ncol(p$k) )
nx <- ncol(x); m <- length(p$coef)
y <- rep(0.0,nx)
for (i in 1:m) {
if (useTerm[i]) { # skip terms when told to
for (j in 1:nx) {
y[j] <- y[j] + p$coef[i]*prod( x[,j]^p$k[i,] )
}
}
}
return(y) }
################################################################################
printPoly <- function( p, num.digits=5 ) {
# print the polynomial p in human readable form
# round coefficients for readability
coef <- round(p$coef,num.digits)
m <- length(coef)
k <- p$k
dimk <- dim(k)
if (m != dimk[1]) stop( "m != dimk[1]; not valid polynomial" )
n <- dimk[2]
cat("Polynomial has",m,"terms and",n,"variables\n")
deg <- rowSums(k)
for (i in 1:m) {
if (i > 1) { if (coef[i] > 0) cat ("\n + ") else cat("\n ") }
cat( coef[i] )
for (j in 1:n) {
if (k[i,j] > 0 ) cat( " * x[",j,"]^", k[i,j], sep="" )
}
cat(" (degree=", deg[i],")" )
}
cat("\n")}
################################################################################
LasserreAvrachenkov <- function( q, p, useTerm, S ){
# Exactly integrate a homogeneous polynomial p over a single simplex S using the
# algorithm in Lasserre and Avrachenkov (2001), MAA Monthly, v 108, n 2, pg 151-154.
# We use the formula (18) of Baldoni, et. al. (2011), Math. of Comput., v 80, pg 297-325.
# We also incorporate the fact that you can eliminate half the terms in their inner sum and
# multiply the result by 2. This is done below by restricting eps[1]=+1, or equivalently
# j[1]=0.
# The polynomial p(.) should be defined using function definePoly.
# Note that p may not be homogeneous, however argument useTerm is a boolean
# vector telling which terms of p to use. It is assumed that all the terms
# in p with useTerm[i]=TRUE are homogeneous of order q. Since this function
# is only called internally from function integrateSimplexPolynomial, we do
# not check this here.
n <- nrow(S); m <- ncol(S)-1
vol <- SimplexVolume(S)
if (q == 0) { # constant terms are treated separately
const <- vol
total <- sum( p$coef[useTerm] )
count <- 1
} else { #
const <- vol/ ( 2^(q-1) * factorial(q) * choose(q+n,n) )
count <- 0
total <- 0.0
i <- rep( 1L, q )
while ( !is.na( i[1] ) ) {
j <- rep( 0L, q )
while ( j[1]==0L ) {
eps <- 1 - 2*j # convert to +1s and -1s
tmp <- eps[1]*S[,i[1]]
if( q > 1 ) {
for (l in 2:q) {
tmp <- tmp + eps[l]*S[ ,i[l]]
}
}
val <- evalPoly( tmp, p, useTerm )
total <- total + prod(eps) * val
count <- count + 1
j <- nextIntBaseB( j, b=2 )
}
i <- nextIndexLA( i, b=m+1 )
}
}
return( list(integral=const*total, functionEvaluations=count) ) }
################################################################################
SimplexVolume <- function( S ){
# compute the volume of the n-dim simplex S in R^n. S is an n x (n+1) matrix,
# with the columns of S giving the vertices of the simplex.
# For S=CanonicalSimplex(n), the volume is 1/n!
stopifnot( is.matrix(S), nrow(S) >= 1, ncol(S) == nrow(S)+1 )
n <- nrow(S)
V <- S[ ,-1] - matrix( S[,1], nrow=n,ncol=n ) # n x n matrix
return( abs( det(V) )/factorial(n) ) }
################################################################################
SimplexSurfaceArea <- function( S3 ){
# compute the surface area of the (n-1)-dim simplex S3 in R^n. S3 is an
# n x n matrix, with the columns of S3 giving the vertices of the simplex.
# For S3 = UnitSimplexV(n),
# SimplexSurfaceArea(UnitSimplexV(n))/SimplexVolume( CanonicalSimplex(n-1) )
# = JacobianS2Canonical(SimplexSurfaceArea(UnitSimplexV(n)))= sqrt(n)
stopifnot( is.matrix(S3), nrow(S3) >= 1, ncol(S3) == nrow(S3) )
n <- nrow(S3)
return( JacobianS2Canonical(S3)*SimplexVolume(CanonicalSimplex(n-1)) ) }
################################################################################
JacobianS2Canonical <- function( S2 ){
# compute the Jacobian of transformation between the m-dimensional simplex S2
# (sitting inside R^n) and the canonical simplex in R^m.
# S2 is an n x (m+1) matrix, with columns of S2 giving the vertices of the simplex.
# The method is to compute the Gram matrix G=t(V) %*% V, where V is as below.
stopifnot( is.matrix(S2), nrow(S2) > 1, ncol(S2) <= nrow(S2)+1 )
n <- nrow(S2); m <- ncol(S2)-1
V <- S2[ ,-1] - matrix( S2[,1], nrow=n,ncol=m ) # an n x m matrix
return( sqrt( abs( det( t(V) %*% V ) ) ) ) }
################################################################################
CanonicalSimplex <- function( n ) {
# return the vertices of the canonical simplex in R^n, the n-dim. simplex with
# vertices the origin (0,...,0) and the unit vectors e[1],...,e[n].
# It has n-dim. volume 1/n!
return( cbind( rep(0.0,n), diag(rep(1.0,n)) ) ) }
################################################################################
UnitSimplexV <- function( n ) {
# return the vertices (V-representaion) of the unit simplex in R^n, the (n-1) dim. simplex with
# vertices the unit vectors e[1],...,e[n].
return( diag(rep(1.0,n)) ) }
################################################################################
nextIntBaseB <- function( current.n, b ) {
# Compute the next integer in base b representation.
# current.n is a vector of integers containing the base b representation of integer n:
# n=current.n[1]*b^k + current.n[2]*b^(k-1) + ... + current.n[2]*b + current.n[1].
# It is assumed that n >= 0.
# Compute (n+1) in base b, by adding 1 and carrying. If necessary, the
# length of the integer vector is increased by 1 to add a new leading digit.
next.n <- current.n
k <- length(next.n)
repeat {
next.n[k] <- next.n[k] + 1
if (next.n[k] >= b) { next.n[k] <- 0; carry <- TRUE
} else { carry <- FALSE }
k <- k-1
if ( !carry | (k==0 ) ) break
}
if (carry) { next.n <- c(1L,next.n) }
return( next.n ) }
################################################################################
nextIndexLA <- function( current.n, b ) {
# Compute the next multi-index used in the outer sum in function LasserreAvrachenkov
# Do this by incrementing current.n until we find one that meets the requirements.
# The number of values returned is choose( b+d, b ); after that
# a vector of NAs is returned to signal that there are no more indices.
d <- length( current.n )
ones <- rep( 1L, d )
next.n <- current.n - ones # subtract off the initial ones
repeat {
next.n <- nextIntBaseB( next.n, b )
if ( all( diff(next.n) >= 0) | (length(next.n) > d ) ) break
}
if (length(next.n) > d) { next.n <- rep( NA, d ) }
else { next.n <- next.n + ones }
return( next.n ) }
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.