Cubature Vectorization Results

Introduction

This R cubature package exposes both the hcubature and pcubature routines of the underlying C cubature library, including the vectorized interfaces.

Per the documentation, use of pcubature is advisable only for smooth integrands in dimensions up to three at most. In fact, the pcubature routines perform significantly worse than the vectorized hcubature in inappropriate cases. So when in doubt, you are better off using hcubature.

Version 2.0 of this package integrates the Cuba library as well, once again providing vectorized interfaces.

The main point of this note is to examine the difference vectorization makes. My recommendations are below in the summary section.

A Timing Harness

Our harness will provide timing results for hcubature, pcubature (where appropriate) and Cuba cuhre calls. We begin by creating a harness for these calls.

library(benchr)
library(cubature)

harness <- function(which = NULL,
                    f, fv, lowerLimit, upperLimit, tol = 1e-3, times = 20, ...) {

    fns <- c(hc = "Non-vectorized Hcubature",
             hc.v = "Vectorized Hcubature",
             pc = "Non-vectorized Pcubature",
             pc.v = "Vectorized Pcubature",
             cc = "Non-vectorized cubature::cuhre",
             cc_v = "Vectorized cubature::cuhre")
    cc <- function() cubature::cuhre(f = f,
                                     lowerLimit = lowerLimit, upperLimit = upperLimit,
                                     relTol = tol,
                                     ...)
    cc_v <- function() cubature::cuhre(f = fv,
                                       lowerLimit = lowerLimit, upperLimit = upperLimit,
                                       relTol = tol,
                                       nVec = 1024L,
                                       ...)

    hc <- function() cubature::hcubature(f = f,
                                         lowerLimit = lowerLimit,
                                         upperLimit = upperLimit,
                                         tol = tol,
                                         ...)

    hc.v <- function() cubature::hcubature(f = fv,
                                           lowerLimit = lowerLimit,
                                           upperLimit = upperLimit,
                                           tol = tol,
                                           vectorInterface = TRUE,
                                           ...)

    pc <- function() cubature::pcubature(f = f,
                                     lowerLimit = lowerLimit,
                                     upperLimit = upperLimit,
                                     tol = tol,
                                     ...)

    pc.v <- function() cubature::pcubature(f = fv,
                                           lowerLimit = lowerLimit,
                                           upperLimit = upperLimit,
                                           tol = tol,
                                           vectorInterface = TRUE,
                                           ...)
    
    ndim = length(lowerLimit)

    if (is.null(which)) {
        fnIndices <- seq_along(fns)
    } else {
        fnIndices <- na.omit(match(which, names(fns)))
    }
    fnList <- lapply(names(fns)[fnIndices], function(x) call(x))

    argList <- c(fnList, times = times, progress = FALSE)
    result <- do.call(benchr::benchmark, args = argList)
    d <- summary(result)[seq_along(fnIndices), ]
    d$expr <- fns[fnIndices]
    d
}

We reel off the timing runs.

Example 1.

func <- function(x) sin(x[1]) * cos(x[2]) * exp(x[3])
func.v <- function(x) {
    matrix(apply(x, 2, function(z) sin(z[1]) * cos(z[2]) * exp(z[3])), ncol = ncol(x))
}

d <- harness(f = func, fv = func.v,
             lowerLimit = rep(0, 3),
             upperLimit = rep(1, 3),
             tol = 1e-5,
             times = 100)
knitr::kable(d, digits = 3, row.names = FALSE)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	100	0.000	0.000	0.000	0.000	0.000	0.003	0.030	1.00
Vectorized Hcubature	100	0.000	0.000	0.000	0.000	0.000	0.003	0.043	1.45
Non-vectorized Pcubature	100	0.001	0.001	0.001	0.001	0.001	0.003	0.091	3.10
Vectorized Pcubature	100	0.001	0.001	0.001	0.001	0.001	0.004	0.125	4.27
Non-vectorized cubature::cuhre	100	0.001	0.001	0.001	0.001	0.001	0.003	0.065	2.19
Vectorized cubature::cuhre	100	0.001	0.001	0.001	0.001	0.001	0.001	0.062	2.21

Multivariate Normal

Using cubature, we evaluate $\int_R\phi(x)dx$ where $\phi(x)$ is the three-dimensional multivariate normal density with mean 0, and variance $\Sigma = \left(\begin{array}{rrr} 1 &\frac{3}{5} &\frac{1}{3}\\ \frac{3}{5} &1 &\frac{11}{15}\\ \frac{1}{3} &\frac{11}{15} & 1 \end{array} \right)$ and $R$ is $[-\frac{1}{2}, 1] \times [-\frac{1}{2}, 4] \times [-\frac{1}{2}, 2].$

We construct a scalar function (my_dmvnorm) and a vector analog (my_dmvnorm_v). First the functions.

m <- 3
sigma <- diag(3)
sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
sigma[3,2] <- sigma[2, 3] <- 11/15
logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
my_dmvnorm <- function (x, mean, sigma, logdet) {
    x <- matrix(x, ncol = length(x))
    distval <- stats::mahalanobis(x, center = mean, cov = sigma)
    exp(-(3 * log(2 * pi) + logdet + distval)/2)
}

my_dmvnorm_v <- function (x, mean, sigma, logdet) {
    distval <- stats::mahalanobis(t(x), center = mean, cov = sigma)
    exp(matrix(-(3 * log(2 * pi) + logdet + distval)/2, ncol = ncol(x)))
}

Now the timing.

d <- harness(f = my_dmvnorm, fv = my_dmvnorm_v,
             lowerLimit = rep(-0.5, 3),
             upperLimit = c(1, 4, 2),
             tol = 1e-5,
             times = 10,
             mean = rep(0, m), sigma = sigma, logdet = logdet)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	10	0.883	0.887	0.893	1.020	1.213	1.263	10.200	491.00
Vectorized Hcubature	10	0.002	0.002	0.002	0.002	0.003	0.003	0.023	1.07
Non-vectorized Pcubature	10	0.378	0.399	0.457	0.473	0.544	0.604	4.733	251.00
Vectorized Pcubature	10	0.001	0.001	0.002	0.002	0.002	0.002	0.018	1.00
Non-vectorized cubature::cuhre	10	0.366	0.367	0.369	0.429	0.519	0.522	4.286	203.00
Vectorized cubature::cuhre	10	0.003	0.004	0.004	0.004	0.005	0.006	0.043	2.00

The effect of vectorization is huge. So it makes sense for users to vectorize the integrands as much as possible for efficiency.

Furthermore, for this particular example, we expect mvtnorm::pmvnorm to do pretty well since it is specialized for the multivariate normal. The good news is that the vectorized versions of hcubature and pcubature are quite competitive if you compare the table above to the one below.

library(mvtnorm)
g1 <- function() pmvnorm(lower = rep(-0.5, m),
                                  upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
                                  alg = Miwa(), abseps = 1e-5, releps = 1e-5)
g2 <- function() pmvnorm(lower = rep(-0.5, m),
                         upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
                         alg = GenzBretz(), abseps = 1e-5, releps = 1e-5)
g3 <- function() pmvnorm(lower = rep(-0.5, m),
                         upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
                         alg = TVPACK(), abseps = 1e-5, releps = 1e-5)

knitr::kable(summary(benchr::benchmark(g1(), g2(), g3(), times = 20, progress = FALSE)),
             digits = 3, row.names = FALSE)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
g1()	20	0.001	0.001	0.001	0.001	0.001	0.003	0.030	1
g2()	20	0.001	0.001	0.001	0.001	0.001	0.004	0.027	1
g3()	20	0.001	0.001	0.001	0.001	0.001	0.003	0.027	1

Product of cosines

testFn0 <- function(x) prod(cos(x))
testFn0_v <- function(x) matrix(apply(x, 2, function(z) prod(cos(z))), ncol = ncol(x))

d <- harness(f = testFn0, fv = testFn0_v,
             lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 1000)
knitr::kable(d, digits = 3)

expr	n.eval	max	total	relative
Non-vectorized Hcubature	1000	0.002	0.047	1.00
Vectorized Hcubature	1000	0.002	0.075	1.67
Non-vectorized Pcubature	1000	0.002	0.057	1.30
Vectorized Pcubature	1000	0.002	0.128	2.91
Non-vectorized cubature::cuhre	1000	0.004	0.379	8.42
Vectorized cubature::cuhre	1000	0.003	0.405	9.15

Gaussian function

testFn1 <- function(x) {
    val <- sum(((1 - x) / x)^2)
    scale <- prod((2 / sqrt(pi)) / x^2)
    exp(-val) * scale
}

testFn1_v <- function(x) {
    val <- matrix(apply(x, 2, function(z) sum(((1 - z) / z)^2)), ncol(x))
    scale <- matrix(apply(x, 2, function(z) prod((2 / sqrt(pi)) / z^2)), ncol(x))
    exp(-val) * scale
}

d <- harness(f = testFn1, fv = testFn1_v,
             lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 10)

knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	10	0.003	0.003	0.003	0.003	0.003	0.003	0.029	32.90
Vectorized Hcubature	10	0.005	0.005	0.005	0.005	0.005	0.007	0.052	55.90
Non-vectorized Pcubature	10	0.000	0.000	0.000	0.000	0.000	0.000	0.001	1.00
Vectorized Pcubature	10	0.000	0.000	0.000	0.000	0.000	0.000	0.002	2.13
Non-vectorized cubature::cuhre	10	0.013	0.014	0.015	0.015	0.016	0.018	0.148	163.00
Vectorized cubature::cuhre	10	0.020	0.020	0.021	0.021	0.022	0.022	0.211	235.00

Discontinuous function

testFn2 <- function(x) {
    radius <- 0.50124145262344534123412
    ifelse(sum(x * x) < radius * radius, 1, 0)
}

testFn2_v <- function(x) {
    radius <- 0.50124145262344534123412
    matrix(apply(x, 2, function(z) ifelse(sum(z * z) < radius * radius, 1, 0)), ncol = ncol(x))
}

d <- harness(which = c("hc", "hc.v", "cc", "cc_v"),
             f = testFn2, fv = testFn2_v,
             lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 10)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	10	0.041	0.042	0.043	0.043	0.044	0.045	0.431	1.00
Vectorized Hcubature	10	0.049	0.049	0.049	0.049	0.049	0.049	0.489	1.15
Non-vectorized cubature::cuhre	10	0.802	0.807	0.811	0.820	0.823	0.862	8.200	19.00
Vectorized cubature::cuhre	10	0.888	0.894	0.899	0.903	0.906	0.951	9.034	21.00

A Simple Polynomial (product of coordinates)

testFn3 <- function(x) prod(2 * x)
testFn3_v <- function(x) matrix(apply(x, 2, function(z) prod(2 * z)), ncol = ncol(x))

d <- harness(f = testFn3, fv = testFn3_v,
             lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 50)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	50	0.000	0.000	0.000	0.000	0.000	0.000	0.003	1.12
Vectorized Hcubature	50	0.000	0.000	0.000	0.000	0.000	0.000	0.005	1.71
Non-vectorized Pcubature	50	0.000	0.000	0.000	0.000	0.000	0.000	0.003	1.00
Vectorized Pcubature	50	0.000	0.000	0.000	0.000	0.000	0.000	0.004	1.53
Non-vectorized cubature::cuhre	50	0.001	0.001	0.001	0.001	0.001	0.003	0.036	11.90
Vectorized cubature::cuhre	50	0.001	0.001	0.001	0.001	0.001	0.001	0.034	11.90

Gaussian centered at $\frac{1}{2}$

testFn4 <- function(x) {
    a <- 0.1
    s <- sum((x - 0.5)^2)
    ((2 / sqrt(pi)) / (2. * a))^length(x) * exp (-s / (a * a))
}

testFn4_v <- function(x) {
    a <- 0.1
    r <- apply(x, 2, function(z) {
        s <- sum((z - 0.5)^2)
        ((2 / sqrt(pi)) / (2. * a))^length(z) * exp (-s / (a * a))
    })
    matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn4, fv = testFn4_v,
             lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 20)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	20	0.001	0.001	0.001	0.001	0.001	0.002	0.028	1.00
Vectorized Hcubature	20	0.002	0.002	0.002	0.002	0.002	0.002	0.036	1.31
Non-vectorized Pcubature	20	0.002	0.002	0.002	0.002	0.002	0.002	0.041	1.50
Vectorized Pcubature	20	0.003	0.003	0.003	0.003	0.003	0.005	0.055	1.89
Non-vectorized cubature::cuhre	20	0.003	0.003	0.003	0.004	0.003	0.006	0.073	2.40
Vectorized cubature::cuhre	20	0.004	0.004	0.004	0.004	0.004	0.006	0.081	2.78

Double Gaussian

testFn5 <- function(x) {
    a <- 0.1
    s1 <- sum((x - 1 / 3)^2)
    s2 <- sum((x - 2 / 3)^2)
    0.5 * ((2 / sqrt(pi)) / (2. * a))^length(x) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
}
testFn5_v <- function(x) {
    a <- 0.1
    r <- apply(x, 2, function(z) {
        s1 <- sum((z - 1 / 3)^2)
        s2 <- sum((z - 2 / 3)^2)
        0.5 * ((2 / sqrt(pi)) / (2. * a))^length(z) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
    })
    matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn5, fv = testFn5_v,
             lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 20)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	20	0.004	0.004	0.004	0.004	0.004	0.006	0.084	1.46
Vectorized Hcubature	20	0.005	0.005	0.005	0.005	0.005	0.005	0.095	1.82
Non-vectorized Pcubature	20	0.003	0.003	0.003	0.003	0.003	0.005	0.055	1.00
Vectorized Pcubature	20	0.003	0.003	0.003	0.003	0.003	0.004	0.068	1.28
Non-vectorized cubature::cuhre	20	0.007	0.007	0.007	0.008	0.009	0.009	0.157	2.83
Vectorized cubature::cuhre	20	0.008	0.008	0.009	0.009	0.009	0.011	0.178	3.27

Tsuda’s Example

testFn6 <- function(x) {
    a <- (1 + sqrt(10.0)) / 9.0
    prod( a / (a + 1) * ((a + 1) / (a + x))^2)
}

testFn6_v <- function(x) {
    a <- (1 + sqrt(10.0)) / 9.0
    r <- apply(x, 2, function(z) prod( a / (a + 1) * ((a + 1) / (a + z))^2))
    matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn6, fv = testFn6_v,
             lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 20)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	20	0.002	0.002	0.002	0.002	0.002	0.002	0.033	1.00
Vectorized Hcubature	20	0.002	0.002	0.002	0.002	0.002	0.005	0.049	1.36
Non-vectorized Pcubature	20	0.008	0.008	0.008	0.008	0.008	0.011	0.169	5.06
Vectorized Pcubature	20	0.010	0.010	0.011	0.011	0.011	0.013	0.217	6.44
Non-vectorized cubature::cuhre	20	0.004	0.004	0.004	0.005	0.005	0.007	0.094	2.75
Vectorized cubature::cuhre	20	0.005	0.005	0.005	0.008	0.005	0.067	0.166	3.00

Morokoff & Calflish Example

testFn7 <- function(x) {
    n <- length(x)
    p <- 1/n
    (1 + p)^n * prod(x^p)
}
testFn7_v <- function(x) {
    matrix(apply(x, 2, function(z) {
        n <- length(z)
        p <- 1/n
        (1 + p)^n * prod(z^p)
    }), ncol = ncol(x))
}

d <- harness(f = testFn7, fv = testFn7_v,
             lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 20)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	20	0.003	0.003	0.004	0.004	0.004	0.006	0.076	1.00
Vectorized Hcubature	20	0.004	0.004	0.004	0.005	0.005	0.007	0.096	1.22
Non-vectorized Pcubature	20	0.008	0.008	0.009	0.009	0.009	0.011	0.176	2.43
Vectorized Pcubature	20	0.009	0.010	0.010	0.010	0.010	0.013	0.208	2.83
Non-vectorized cubature::cuhre	20	0.044	0.044	0.044	0.045	0.045	0.052	0.896	12.60
Vectorized cubature::cuhre	20	0.043	0.044	0.044	0.044	0.044	0.047	0.883	12.50

Wang-Landau Sampling 1d, 2d Examples

I.1d <- function(x) {
    sin(4 * x) *
        x * ((x * ( x * (x * x - 4) + 1) - 1))
}
I.1d_v <- function(x) {
    matrix(apply(x, 2, function(z)
        sin(4 * z) *
        z * ((z * ( z * (z * z - 4) + 1) - 1))),
        ncol = ncol(x))
}
d <- harness(f = I.1d, fv = I.1d_v,
             lowerLimit = -2, upperLimit = 2, times = 100)
knitr::kable(d, digits = 3)

expr	n.eval	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	100	0.000	0.000	0.000	0.000	0.000	0.013	2.29
Vectorized Hcubature	100	0.000	0.000	0.000	0.000	0.000	0.022	3.75
Non-vectorized Pcubature	100	0.000	0.000	0.000	0.000	0.000	0.006	1.00
Vectorized Pcubature	100	0.000	0.000	0.000	0.000	0.000	0.016	2.77
Non-vectorized cubature::cuhre	100	0.000	0.000	0.000	0.000	0.000	0.024	4.05
Vectorized cubature::cuhre	100	0.001	0.001	0.001	0.001	0.004	0.059	8.92

I.2d <- function(x) {
    x1 <- x[1]; x2 <- x[2]
    sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
}
I.2d_v <- function(x) {
    matrix(apply(x, 2,
                 function(z) {
                     x1 <- z[1]; x2 <- z[2]
                     sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
                 }),
           ncol = ncol(x))
}
d <- harness(f = I.2d, fv = I.2d_v,
             lowerLimit = rep(-1, 2), upperLimit = rep(1, 2), times = 100)
knitr::kable(d, digits = 3)

expr	n.eval	min	lw.qu	median	mean	up.qu	max	total	relative
Non-vectorized Hcubature	100	0.004	0.005	0.005	0.005	0.005	0.009	0.494	10.50
Vectorized Hcubature	100	0.005	0.005	0.006	0.006	0.006	0.009	0.580	12.70
Non-vectorized Pcubature	100	0.000	0.000	0.000	0.000	0.000	0.000	0.044	1.00
Vectorized Pcubature	100	0.001	0.001	0.001	0.001	0.001	0.001	0.066	1.49
Non-vectorized cubature::cuhre	100	0.001	0.001	0.001	0.001	0.001	0.004	0.136	2.95
Vectorized cubature::cuhre	100	0.001	0.001	0.001	0.001	0.001	0.004	0.146	3.20

Implementation Notes

About the only real modification we have made to the underlying cubature library is that we use M = 16 rather than the default M = 19 suggested by the original author for pcubature. This allows us to comply with CRAN package size limits and seems to work reasonably well for the above tests. Future versions will allow for such customization on demand.

The changes made to the Cuba library are managed in a Github repo branch: each time a new release is made, we update the main branch, and keep all changes for Unix platforms in a branch named R_pkg against the current main branch. Customization for windows is done in the package itself using the Makevars.win script.

Summary

The recommendations are:

Vectorize your function. The time spent in so doing pays back enormously. This is easy to do and the examples above show how.
Vectorized hcubature seems to be a good starting point.
For smooth integrands in low dimensions ( $\leq 3$ ), pcubature might be worth trying out. Experiment before using in a production package.

Session Info

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] mvtnorm_1.2-5  cubature_2.1.1 benchr_0.2.5  
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.36      desc_1.4.3         R6_2.5.1           fastmap_1.2.0     
##  [5] xfun_0.45          cachem_1.1.0       knitr_1.48         htmltools_0.5.8.1 
##  [9] rmarkdown_2.27     lifecycle_1.0.4    cli_3.6.3          RcppProgress_0.4.2
## [13] sass_0.4.9         pkgdown_2.1.0      textshaping_0.4.0  jquerylib_0.1.4   
## [17] systemfonts_1.1.0  compiler_4.4.1     tools_4.4.1        ragg_1.3.2        
## [21] evaluate_0.24.0    bslib_0.7.0        Rcpp_1.0.12        yaml_2.3.9        
## [25] jsonlite_1.8.8     rlang_1.1.4        fs_1.6.4

Balasubramanian Narasimhan

2024-07-11

Introduction

A Timing Harness

Example 1.

Multivariate Normal

Product of cosines

Gaussian function

Discontinuous function

A Simple Polynomial (product of coordinates)

Gaussian centered at $\frac{1}{2}$

Double Gaussian

Tsuda’s Example

Morokoff & Calflish Example

Wang-Landau Sampling 1d, 2d Examples

Implementation Notes

Summary

Session Info

Cubature Vectorization Results

Balasubramanian Narasimhan

2024-07-11

Introduction

A Timing Harness

Example 1.

Multivariate Normal

Product of cosines

Gaussian function

Discontinuous function

A Simple Polynomial (product of coordinates)

Gaussian centered at 12\frac{1}{2}

Double Gaussian

Tsuda’s Example

Morokoff & Calflish Example

Wang-Landau Sampling 1d, 2d Examples

Implementation Notes

Summary

Session Info

Gaussian centered at $\frac{1}{2}$