Class inte_hcubature (o2scl)

O2scl : Class List

template<class func_t>
class o2scl::inte_hcubature : public o2scl::inte_cubature_base

Adaptive multidimensional integration on hyper-rectangles using cubature rules from the Cubature library.

This class is experimental.

Documentation adapted from Cubature

A cubature rule takes a function and a hypercube and evaluates the function at a small number of points, returning an estimate of the integral as well as an estimate of the error, and also a suggested dimension of the hypercube to subdivide.

Given such a rule, the adaptive integration is simple:

1) Evaluate the cubature rule on the hypercube(s). Stop if converged.

2) Pick the hypercube with the largest estimated error, and divide it in two along the suggested dimension.

3) Goto (1).

The basic algorithm is based on the adaptive cubature described in [Genz80] and subsequently extended to integrating a vector of integrands in [Berntsen91].

Binary heap implementation

Based on [Cormen09] and used as a priority queue of regions to integrate.

typedef region heap_item

Desc.

int use_parallel

Desc.

void heap_resize(heap &h, size_t nalloc)

Desc.

heap heap_alloc(size_t nalloc, size_t fdim)

Desc.

void heap_free(heap &h)

Note that heap_free does not deallocate anything referenced by the items.

int heap_push(heap &h, heap_item &hi)

Desc.

int heap_push_many(heap &h, size_t ni, heap_item *hi)

Desc.

heap_item heap_pop(heap &h)

Desc.

int converged(size_t fdim, const std::vector<esterr> &ee, double reqAbsError, double reqRelError, error_norm norm)

Desc.

template<class vec_t>
int rulecubature(rule &r, size_t fdim, func_t &f, const hypercube &h, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)

Desc.

template<class vec_t>
int cubature(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err, int parallel)

Desc.

inte_hcubature()
template<class vec_t>
int integ(size_t fdim, func_t &f, size_t dim, const vec_t &xmin, const vec_t &xmax, size_t maxEval, double reqAbsError, double reqRelError, error_norm norm, vec_t &val, vec_t &err)

Desc.

Public Types

typedef boost::numeric::ublas::vector<double> ubvector
typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range

Protected Functions

double errMax(const std::vector<esterr> &ee)

Return the maximum error from ee.

double compute_vol(const hypercube &h)

Desc.

template<class vec_t>
void make_hypercube(size_t dim, const vec_t &center, const vec_t &halfwidth, hypercube &h)

Desc.

void make_hypercube2(size_t dim, const std::vector<double> &dat, hypercube &h)

Desc.

template<class vec_t>
void make_hypercube_range(size_t dim, const vec_t &xmin, const vec_t &xmax, hypercube &h)

Desc.

void destroy_hypercube(hypercube &h)

Desc.

void make_region(const hypercube &h, size_t fdim, region &R)

Desc.

int cut_region(region &R, region &R2)

Desc.

void alloc_rule_pts(rule &r, size_t num_regions)

Desc.

void make_rule(size_t dim, size_t fdim, size_t num_points, rule &r)

Desc.

int eval_regions(size_t nR, std::vector<region> &R, func_t &f, rule &r)

Desc.

Note

All regions must have same fdim

size_t ls0(size_t n)

Functions to loop over points in a hypercube.

Based on orbitrule.cpp in HIntLib-0.0.10

ls0 returns the least-significant 0 bit of n (e.g. it returns 0 if the LSB is 0, it returns 1 if the 2 LSBs are 01, etcetera).

void evalR_Rfs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector<double> &p, size_t p_ix, const std::vector<double> &c, size_t c_ix, const std::vector<double> &r, size_t r_ix)

Evaluate the integration points for all \( 2^n \) points (+/-r,…+/-r)

A Gray-code ordering is used to minimize the number of coordinate updates in p, although this doesn’t matter as much now that we are saving all pts.

void evalRR0_0fs2(ubvector &pts, size_t pts_ix, size_t dim, std::vector<double> &p, size_t p_ix, const std::vector<double> &c, size_t c_ix, const std::vector<double> &r, size_t r_ix)

Desc.

void evalR0_0fs4d2(ubvector &pts, size_t pts_ix, size_t dim, std::vector<double> &p, size_t p_ix, const std::vector<double> &c, size_t c_ix, const std::vector<double> &r1, size_t r1_ix, const std::vector<double> &r2, size_t r2_ix)

Desc.

size_t num0_0(size_t dim)

Desc.

size_t numR0_0fs(size_t dim)

Desc.

size_t numRR0_0fs(size_t dim)

Desc.

size_t numR_Rfs(size_t dim)

Desc.

int rule75genzmalik_evalError(rule &runder, size_t fdim, func_t &f, size_t nR, std::vector<region> &R)

Desc.

void make_rule75genzmalik(size_t dim, size_t fdim, rule75genzmalik &r)

Desc.

int rule15gauss_evalError(rule &r, size_t fdim, func_t &f, size_t nR, std::vector<region> &R)

1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in GNU GSL (which in turn is based on QUADPACK).

void make_rule15gauss(size_t dim, size_t fdim, rule &r)

Desc.

class esterr

A value and error.

Public Functions

esterr()
esterr(const esterr &e)
esterr &operator=(const esterr &e)

Public Members

double val

Desc.

double err

Desc.

class heap

Desc.

Public Functions

heap()
heap(const heap &e)
heap &operator=(const heap &e)

Public Members

size_t n

Desc.

size_t nalloc

Desc.

std::vector<heap_item> items

Desc.

size_t fdim

Desc.

std::vector<esterr> ee

array of length fdim of the total integrand & error

class hypercube

Desc.

Public Functions

hypercube()
hypercube(const hypercube &e)
hypercube &operator=(const hypercube &e)

Public Members

size_t dim

Desc.

std::vector<double> data

length 2*dim = center followed by half-widths

double vol

cache volume = product of widths

class region

Desc.

Public Functions

region()
region(const region &e)
region &operator=(const region &e)

Public Members

hypercube h

Desc.

size_t splitDim

Desc.

size_t fdim

dimensionality of vector integrand

std::vector<esterr> ee

array of length fdim

double errmax

max ee[k].err

class rule

Desc.

Subclassed by o2scl::inte_hcubature< func_t >::rule75genzmalik

Public Functions

rule()
rule(const rule &e)
rule &operator=(const rule &e)

Public Members

size_t dim

The dimensionality.

size_t fdim

The number of functions.

size_t num_points

The number of evaluation points.

size_t num_regions

The max number of regions evaluated at once.

ubvector pts

points to eval: num_regions * num_points * dim

size_t vals_ix

num_regions * num_points * fdim

class rule75genzmalik : public o2scl::inte_hcubature<func_t>::rule

Desc.

Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded cubature rule of degree 7 (embedded rule degree 5) from [Genz83].

Public Functions

rule75genzmalik()
rule75genzmalik(const rule75genzmalik &e)
rule75genzmalik &operator=(const rule75genzmalik &e)

Public Members

std::vector<double> p

Desc.

double weight1

dimension-dependent constants

double weight3

Desc.

double weight5

Desc.

double weightE1

Desc.

double weightE3

Desc.