MWFunctions

Everything that is discussed in the following chapter is available to the application program by including:

#include "MRCPP/MWFunctions"

Multiwavelet (MW) representations of real-valued scalar functions are in MRCPP called FunctionTrees. These are in principle available in any dimension using the template parameter D (in practice D=1,2,3). There are several different ways of constructing MW functions (computing the expansion coefficients in the MW basis):

  • Projection of analytic function

  • Arithmetic operations

  • Application of MW operator

The first two will be discribed in this chapter, while the last one regarding operators will be the topic of the next chapter.

The interface for constructing MW representations in MRCPP has a dual focus: on the one hand we want a simple, intuitive way of producing adaptive numerical approximations with guaranteed precision that does not require detailed knowledge of the internals of the MW code and with a minimal number of parameters that have to be set. On the other hand we want the possibility for more detailed control of the construction and refinement of the numerical grid where such control is possible and even necessary. In the latter case it is important to be able to reuse the existing grids in e.g. iterative algorithms without excessive allocation/deallocation of memory.

MultiResolution Analysis (MRA)

template<int D>
class MultiResolutionAnalysis

Class collecting computational domain and MW basis.

In order to combine different functions and operators in mathematical operations, they need to be compatible. That is, they must be defined on the same computational domain and constructed using the same polynomial basis (order and type). This information constitutes an MRA, which needs to be defined and passed as argument to all function and operator constructors, and only functions and operators with compatible MRAs can be combined in subsequent calculations.

Public Functions

MultiResolutionAnalysis(std::array<int, 2> bb, int order, int depth = MaxDepth)

Constructs a MultiResolutionAnalysis object composed of computational domain (world) and a polynomial basis (Multiwavelets)

Constructor of the MultiResolutionAnalysis class from scratch, without requiring any pre-existing complex structure. The constructor calls the InterpolatingBasis basis constructor to generate the MultiWavelets basis of functions, then the BoundingBox constructor to create the computational domain. The constructor then checks if the generated node depth, or node refinement is beyond the root scale or the maximum depth allowed, in which case it will abort the process. Otherwise, the process goes on to setup the filters with the class’ setupFilter method.

Parameters:

  • bb – [in] 2-element integer array [Lower, Upper] defining the bounds for a BoundingBox object representing the computational domain

  • order – [in] Maximum polynomial order of the multiwavelet basis, immediately used in the constructor of an InterPolatingBasis object which becomes an attribute of the MRA

  • maxDepth – [in] Exponent of the node refinement in base 2, relative to root scale. In other words, it is the maximum amount of refinement that we allow in a node, in other to avoid overflow of values.

Returns:

New MultiResolutionAnalysis (MRA) object

MultiResolutionAnalysis(const BoundingBox<D> &bb, int order, int depth = MaxDepth)

Constructs a MultiResolutionAnalysis object composed of computational domain (world) and a polynomial basis (Multiwavelets) from a pre-existing BoundingBox object.

Constructor of the MultiResolutionAnalysis class from a BoundingBox object. For more details see the first constructor.

Parameters:

  • bb – [in] BoundingBox object representing the computational domain

  • order – [in] (integer) Maximum polynomial order of the multiwavelet basis, immediately used in the constructor of an InterPolatingBasis object which becomes an attribute of the MRA

  • maxDepth – [in] (integer) Exponent of the node refinement in base 2, relative to root scale. In other words, it is the maximum amount of refinement that we allow in a node, in other to avoid overflow of values.

Returns:

New MultiResolutionAnalysis (MRA) object

MultiResolutionAnalysis(const BoundingBox<D> &bb, const ScalingBasis &sb, int depth = MaxDepth)

Constructor for a MultiResolutionAnalysis object from a pre-existing BoundingBox (computational domain) and a ScalingBasis (Multiwavelet basis) objects.

Creates a MRA object from pre-existing BoundingBox and ScalingBasis objects. These objects are taken as reference. For more details about the constructor itself, see the first constructor.

Parameters:

  • bb – [in] Computational domain as a BoundingBox object, taken by constant reference

  • sb – [in] Polynomial basis (MW) as a ScalingBasis object

  • depth – [in] Maximum allowed resolution depth, relative to root scale

Returns:

New MultiResolutionAnalysis object

MultiResolutionAnalysis(const MultiResolutionAnalysis<D> &mra)

Copy constructor for a MultiResolutionAnalysis object composed of computational domain (world) and a polynomial basis (Multiwavelets)

Copy a MultiResolutionAnalysis object without modifying the original. For more details see the first constructor.

Parameters:

mra – [in] Pre-existing MRA object

Returns:

New MultiResolutionAnalysis (MRA) object

double calcMaxDistance() const

Computes the difference between the lower and upper bounds of the computational domain.

Returns:

Maximum possible distance between two points in the MRA domain

bool operator==(const MultiResolutionAnalysis<D> &mra) const

Equality operator for the MultiResolutionAnalysis class, returns true if both MRAs have the same polynomial basis, computational domain and maximum depth, and false otherwise.

Equality operator for the MultiResolutionAnalysis class, returns true if both MRAs have the same polynomial basis represented by a BoundingBox object, computational domain (ScalingBasis object) and maximum depth (integer), and false otherwise. Computations on different MRA cannot be combined, this operator can be used to make sure that the multiple MRAs are compatible. For more information about the meaning of equality for BoundingBox and ScalingBasis objets, see their respective classes.

Parameters:

mra – [in] MRA object, taken by constant reference

Returns:

Whether the two MRA objects are equal.

bool operator!=(const MultiResolutionAnalysis<D> &mra) const

Inequality operator for the MultiResolutionAnalysis class, returns false if both MRAs have the same polynomial basis, computational domain and maximum depth, and true otherwise.

Inequality operator for the MultiResolutionAnalysis class, returns true if both MRAs have the same polynomial basis represented by a BoundingBox object, computational domain (ScalingBasis object) and maximum depth (integer), and false otherwise. Opposite of the == operator. For more information about the meaning of equality for BoundingBox and ScalingBasis objets, see their respective classes.

Parameters:

mra – [in] MRA object, taken by constant reference

Returns:

Whether the two MRA objects are not equal.

void print() const

Displays the MRA’s attributes in the outstream defined in the Printer class.

This function displays the attributes of the MRA in the using the Printer class. By default, the Printer class writes all information in the output file, not the terminal.

template<int D>
class BoundingBox

Class defining the computational domain.

The computational domain is made up of a collection of D-dimensional boxes on a particular length scale \( n \). The size of each box is then \( [2^{-n}]^D \), i.e. higher scale means smaller boxes, and the scale may be negative. The number of boxes can be different in each dimension \( [n_x, n_y, \dots] \), but they must all be on the same scale (size). Box translations relative to the world origin must be an integer multiple of the given scale size \( 2^{-n} \).

Subclassed by mrcpp::NodeBox< D, T >

Public Functions

explicit BoundingBox(std::array<int, 2> box)

Constructor for BoundingBox object.

Creates a box with appropriate root scale and scaling factor to fit the given bounds, which applies to all dimensions. Root scale is chosen such that the scaling factor becomes 1 < sfac < 2.

Limitations: Box must be either [0,L] or [-L,L], with L a positive integer. This is the most general constructor, which will create a world with no periodic boundary conditions.

Parameters:

box – [in] [lower, upper] bound in all dimensions.

Returns:

New BoundingBox object.

explicit BoundingBox(int n = 0, const std::array<int, D> &l = {}, const std::array<int, D> &nb = {}, const std::array<double, D> &sf = {}, bool pbc = false)

Constructor for BoundingBox object.

Creates a box with given parameters. The parameter n defines the length scale, which, together with sf, determines the unit length of each side of the boxes by \( [2^{-n}]^D \). The parameter l defines the corner translation of the lower corner of the box relative to the world origin. The parameter nb defines the number of boxes in each dimension. The parameter sf defines the scaling factor, which determines the box translations around the origin, i.e. the amount of boxes around origin. The parameter pbc defines whether the world is periodic or not. In this constructor this value makes all dimensions periodic. This constructor is used for work in periodic systems.

Parameters:

  • n – [in] Length scale, default 0.

  • l – [in] Corner translation, default [0, 0, …].

  • nb – [in] Number of boxes, default [1, 1, …].

  • sf – [in] Scaling factor, default [1.0, 1.0, …].

  • pbc – [in] Periodic boundary conditions, default false.

Returns:

New BoundingBox object.

explicit BoundingBox(const NodeIndex<D> &idx, const std::array<int, D> &nb = {}, const std::array<double, D> &sf = {})

Constructor for BoundingBox object.

Creates a box with given parameters. The parameter idx defines the index of the lower corner of the box relative to the world origin. The parameter nb defines the number of boxes in each dimension. The parameter sf defines the scaling factor, which determines the box translations around the origin, i.e. the amount of boxes around origin. This constructor creates a world with no periodic boundary conditions.

Parameters:

  • idx – [in] index of the lower corner of the box.

  • nb – [in] Number of boxes, default [1, 1, …].

  • sf – [in] Scaling factor, default [1.0, 1.0, …].

Returns:

New BoundingBox object.

explicit BoundingBox(const std::array<double, D> &sf, bool pbc = true)

Constructor for BoundingBox object.

Creates a box with given parameters. The parameter sf defines the scaling factor, which determines the box translations around the origin, i.e. the amount of boxes around origin. The parameter pbc defines whether the world is periodic or not. In this constructor this value makes all dimensions periodic. This construtor is used for work in periodic systems.

Parameters:

  • sf – [in] Scaling factor, default [1.0, 1.0, …].

  • pbc – [in] Periodic boundary conditions, default true.

BoundingBox(const std::array<double, D> &sf, std::array<bool, D> pbc)

Constructor for BoundingBox object.

Creates a box with given parameters. The parameter sf defines the scaling factor, which determines the box translations around the origin, i.e. the amount of boxes around origin. The parameter pbc defines whether the world is periodic or not. In this constructor this value makes specific dimensions periodic. This is used for work in periodic systems.

Parameters:

  • sf – [in] Scaling factor, default [1.0, 1.0, …].

  • pbc – [in] Periodic boundary conditions, default true.

Returns:

New BoundingBox object.

BoundingBox(const BoundingBox<D> &box)

Constructor for BoundingBox object.

Creates a box identical to the input box paramter. This constructor uses all parameters from the other BoundingBox to create a new one.

Parameters:

box – [in] Other BoundingBox object.

Returns:

New BoundingBox object.

BoundingBox<D> &operator=(const BoundingBox<D> &box)

Assignment operator overload for BoundingBox object.

Allocates all parameters in this BoundingBox to be that of the other BoundingBox.

Parameters:

box – [in] Other BoundingBox object.

Returns:

New BoundingBox object.

NodeIndex<D> getNodeIndex(int bIdx) const

Fetches a NodeIndex object from a given box index.

During the adaptive refinement, each original box will contain an increasing number of smaller cells, each of which will be part of a specific node in the tree. These cells are divided adaptivelly. This function returns the NodeIndex object of the cell at the lower back corner of the box object indexed by bIdx. Specialized for D=1 below

Parameters:

bIdx – [in] Box index, the index of the box we want to fetch the cell index from.

Returns:

The NodeIndex object of the index given as it is in the Multiresolutoin analysis.

int getBoxIndex(Coord<D> r) const

Fetches the index of a box from a given coordinate.

Specialized for D=1 below

Parameters:

r – [in] D-dimensional array representaing a coordinate in the simulation box

Returns:

The index value of the boxes in the position given as it is in the generated world.

int getBoxIndex(NodeIndex<D> nIdx) const

Fetches the index of a box from a given NodeIndex.

During the multiresolution analysis the boxes will be divided into smaller boxes, which means that each individual box will be part of a specific node in the tree. Each node will get its own index value, but will still be part of one of the original boxes of the world. Specialized for D=1 below

Parameters:

nIdx – [in] NodeIndex object, representing the node and its index in the adaptive tree.

Returns:

The index value of the boxes in which the NodeIndex object is mapping to.

NodeIndex<1> getNodeIndex(int bIdx) const

Fetches a NodeIndex object from a given box index, specialiced for 1-D.

During the adaptive refinement, each original box will contain an increasing number of smaller cells, each of which will be part of a specific node in the tree. These cells are divided adaptivelly. This function returns the NodeIndex object of the cell at the lower back corner of the box object indexed by bIdx.

Parameters:

bIdx – [in] Box index, the index of the box we want to fetch the cell index from.

Returns:

The NodeIndex object of the index given as it is in the Multiresolutoin analysis.

int getBoxIndex(Coord<1> r) const

Fetches the index of a box from a given coordinate, specialized for 1D.

Parameters:

r – [in] 1-dimensional array representaing a coordinate in the simulation box

Returns:

The index value of the boxes in the position given as it is in the generated world.

int getBoxIndex(NodeIndex<1> nIdx) const

Fetches the index of a box from a given NodeIndex specialized for 1-D.

During the multiresolution analysis the boxes will be divided into smaller boxes, which means that each individual box will be part of a specific node in the tree. Each node will get its own index value, but will still be part of one of the original boxes of the world.

Parameters:

nIdx – [in] NodeIndex object, representing the node and its index in the adaptive tree.

Returns:

The index value of the boxes in which the NodeIndex object is mapping to.

class LegendreBasis : public mrcpp::ScalingBasis

Legendre scaling functions as defined by Alpert, SIAM J Math Anal 24 (1), 246 (1993).

Public Functions

inline LegendreBasis(int k)
Parameters:

k – [in] Polynomial order of basis, 1 < k < 40

Returns:

New LegendreBasis object

class InterpolatingBasis : public mrcpp::ScalingBasis

Interpolating scaling functions as defined by Alpert etal, J Comp Phys 182, 149-190 (2002).

Public Functions

inline InterpolatingBasis(int k)
Parameters:

k – [in] Polynomial order of basis, 1 < k < 40

Returns:

New InterpolatingBasis object

FunctionTree

template<int D, typename T>
class FunctionTree : public mrcpp::MWTree<D, T>, public mrcpp::RepresentableFunction<D, T>

Function representation in MW basis.

Constructing a full grown FunctionTree involves a number of steps, including setting up a memory allocator, constructing root nodes according to the given MRA, building an adaptive tree structure and computing MW coefficients. The FunctionTree constructor does only half of these steps: It takes an MRA argument, which defines the computational domain and scaling basis (these are fixed parameters that cannot be changed after construction). The tree is initialized with a memory allocator and a set of root nodes, but it does not compute any coefficients and the function is initially undefined. An undefined FunctionTree will have a well defined tree structure (at the very least the root nodes of the given MRA, but possibly with additional refinement) and its MW coefficient will be allocated but uninitialized, and its square norm will be negative (minus one).

Public Functions

inline FunctionTree(const MultiResolutionAnalysis<D> &mra, const std::string &name)

Constructor.

Constructs an uninitialized tree, containing only empty root nodes. Nonshared memory is used for the tree, which means that the tree is allocated locally.

Parameters:

  • mra – [in] Which MRA the function is defined

  • name – [in] Name of the tree

Returns:

New FunctionTree object

FunctionTree(const MultiResolutionAnalysis<D> &mra, SharedMemory<T> *sh_mem = nullptr, const std::string &name = "nn")

Constructor.

Constructs an uninitialized tree, containing only empty root nodes. If a shared memory pointer is provided the tree will be allocated in this shared memory window, otherwise it will be local to each MPI process.

Constructs an uninitialized tree, containing only empty root nodes. If a shared memory pointer is provided the tree will be allocated in this shared memory window, otherwise it will be local to each MPI process.

Parameters:

  • mra – [in] Which MRA the function is defined

  • sh_mem – [in] Pointer to MPI shared memory block

  • name – [in] Name of the tree

  • mra – [in] Which MRA the function is defined

  • sh_mem – [in] Pointer to MPI shared memory block

Returns:

New FunctionTree object

Returns:

New FunctionTree object

template<typename U = T, typename = std::enable_if_t<std::is_same_v<U, ComplexDouble>>>
FunctionTree(FunctionTree<D, double> &realTree, FunctionTree<D, double> &imagTree, SharedMemory<T> *sh_mem = nullptr, const std::string &name = "nn")

Constructor (complex version only)

Constructs a complex tree, given two real trees which define the real and imaginary parts of the complex function. If a shared memory pointer is provided the tree will be allocated in this shared memory window, otherwise it will be local to each MPI process.

Parameters:

  • mra – [in] Which MRA the function is defined

  • sh_mem – [in] Pointer to MPI shared memory block

  • name – [in] Name of the tree

Returns:

New FunctionTree object

Creating defined FunctionTrees

The following functions will define MW coefficients where there are none, and thus require that the output FunctionTree is in an undefined state. All functions marked with ‘adaptive grid’ will use the same building algorithm:

  1. Start with an initial guess for the grid

  2. Compute the MW coefficients for the output function on the current grid

  3. Refine the grid where necessary based on the local wavelet norm

  4. Iterate points 2 and 3 until the grid is converged

With a negative precision argument, the grid will be fixed, e.i. it will not be refined beyond the initial grid. There is also an argument to limit the number of extra refinement levels beyond the initial grid, in which the adaptive refinement will stop, even if the local precision requirement is not met.

void mrcpp::MWTree::setZero()

Set the MW coefficients to zero, keeping the same tree structure.

Keeps the node structure of the tree, even though the zero function is representable at depth zero. One should then use cropTree to remove unnecessary nodes.

Warning

doxygenfunction: Unable to resolve function “mrcpp::project” with arguments “(double, FunctionTree<D>&, RepresentableFunction<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void project (CompFunction< 3 > &out, RepresentableFunction< 3, double > &f, double prec) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::copy_func” with arguments (FunctionTree<D>&, FunctionTree<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void copy_func(FunctionTree<D, T> &out, FunctionTree<D, T> &inp)

Warning

doxygenfunction: Unable to resolve function “mrcpp::add” with arguments “(double, FunctionTree<D>&, FunctionTreeVector<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void add (CompFunction< 3 > &out, ComplexDouble a, CompFunction< 3 > inp_a, ComplexDouble b, CompFunction< 3 > inp_b, double prec, bool conjugate) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::add” with arguments “(double, FunctionTree<D>&, double, FunctionTree<D>&, double, FunctionTree<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void add (CompFunction< 3 > &out, ComplexDouble a, CompFunction< 3 > inp_a, ComplexDouble b, CompFunction< 3 > inp_b, double prec, bool conjugate) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::multiply” with arguments “(double, FunctionTree<D>&, FunctionTreeVector<D>&, int, bool, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void multiply (CompFunction< 3 > &out, CompFunction< 3 > inp_a, CompFunction< 3 > inp_b, double prec, bool absPrec, bool useMaxNorms, bool conjugate) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::multiply” with arguments “(double, FunctionTree<D>&, double, FunctionTree<D>&, FunctionTree<D>&, int, bool, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void multiply (CompFunction< 3 > &out, CompFunction< 3 > inp_a, CompFunction< 3 > inp_b, double prec, bool absPrec, bool useMaxNorms, bool conjugate) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::square” with arguments (double, FunctionTree<D>&, FunctionTree<D>&, int, bool) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void square(double prec, FunctionTree<D, T> &out, FunctionTree<D, T> &inp, int maxIter, bool absPrec, bool conjugate)

Warning

doxygenfunction: Unable to resolve function “mrcpp::power” with arguments (double, FunctionTree<D>&, FunctionTree<D>&, double, int, bool) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void power(double prec, FunctionTree<D, T> &out, FunctionTree<D, T> &inp, double p, int maxIter, bool absPrec)

Warning

doxygenfunction: Unable to resolve function “mrcpp::dot” with arguments “(double, FunctionTree<D>&, FunctionTreeVector<D>&, FunctionTreeVector<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template ComplexDouble dot (CompFunction< 3 > bra, CompFunction< 3 > ket) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::treeMap” with arguments (double, FunctionTree<D>&, FunctionTree<D>&, FMap, int, bool) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void treeMap(double prec, FunctionTree<D, T> &out, FunctionTree<D, T> &inp, FMap<T, T> fmap, int maxIter, bool absPrec)

Creating undefined FunctionTrees

The grid of a FunctionTree can also be constructed without computing any MW coefficients:

Warning

doxygenfunction: Unable to resolve function “mrcpp::build_grid” with arguments (FunctionTree<D>&, const RepresentableFunction<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, int scales)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, std::vector<FunctionTree<D, T>*> &inp, int maxIter)
- template<int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, int maxIter)
template<int D>
void mrcpp::build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, int maxIter)

Build empty grid based on info from Gaussian expansion.

The grid of the output function will be EXTENDED using the general algorithm:

  • Loop through current leaf nodes of the output tree

  • Refine node based on custom split check from the function

  • Repeat until convergence or maxIter is reached

  • maxIter < 0 means no bound

Note

This algorithm will start at whatever grid is present in the out tree when the function is called. It will loop through the Gaussians in the expansion and extend the grid based on the position and exponent of each term. Higher exponent means finer resolution.

Parameters:

  • out – [out] Output tree to be built

  • inp – [in] Input Gaussian expansion

  • maxIter – [in] Maximum number of refinement iterations in output tree

Warning

doxygenfunction: Unable to resolve function “mrcpp::build_grid” with arguments (FunctionTree<D>&, FunctionTree<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, int scales)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, std::vector<FunctionTree<D, T>*> &inp, int maxIter)
- template<int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, int maxIter)

Warning

doxygenfunction: Unable to resolve function “mrcpp::build_grid” with arguments (FunctionTree<D>&, FunctionTreeVector<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp, int maxIter)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, int scales)
- template<int D, typename T> void build_grid(FunctionTree<D, T> &out, std::vector<FunctionTree<D, T>*> &inp, int maxIter)
- template<int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, int maxIter)

Warning

doxygenfunction: Unable to resolve function “mrcpp::copy_grid” with arguments “(FunctionTree<D>&, FunctionTree<D>&)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void copy_grid (CompFunction< 1 > &out, CompFunction< 1 > &inp) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::clear_grid” with arguments (FunctionTree<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void clear_grid(FunctionTree<D, T> &out)
void mrcpp::FunctionTree::clear()

Remove all nodes in the tree.

Leaves the tree in the same state as after construction, i.e. undefined tree structure containing only root nodes without coefficients. The assigned memory, including branch and leaf nodes, (nodeChunks in NodeAllocator) is NOT released, but is immediately available to the new function.

Changing FunctionTrees

There are also a number of in-place operations that change the MW coefficients of a given defined FunctionTree. All changing operations require that the FunctionTree is in a defined state.

void mrcpp::FunctionTree::rescale(T c)

In-place multiplication by a scalar, fixed grid.

The leaf node point values of the function will be in-place multiplied by the given coefficient, no grid refinement.

Parameters:

c – [in] Scalar coefficient

void mrcpp::FunctionTree::normalize()

In-place rescaling by a function norm \( ||f||^{-1} \), fixed grid.

void mrcpp::FunctionTree::add(T c, FunctionTree<D, T> &inp)

In-place addition with MW function representations, fixed grid.

The input function will be added in-place on the current grid of the function, i.e. no further grid refinement.

Parameters:

  • c – [in] Numerical coefficient of input function

  • inp – [in] Input function to add

void mrcpp::FunctionTree::multiply(T c, FunctionTree<D, T> &inp)

In-place multiplication with MW function representations, fixed grid.

The input function will be multiplied in-place on the current grid of the function, i.e. no further grid refinement.

Parameters:

  • c – [in] Numerical coefficient of input function

  • inp – [in] Input function to multiply

void mrcpp::FunctionTree::square()

In-place square of MW function representations, fixed grid.

The leaf node point values of the function will be in-place squared, no grid refinement.

void mrcpp::FunctionTree::power(double p)

In-place power of MW function representations, fixed grid.

The leaf node point values of the function will be in-place raised to the given power, no grid refinement.

Parameters:

p – [in] Numerical power

void mrcpp::FunctionTree::map(FMap<T, T> fmap)

In-place mapping with a predefined function f(x), fixed grid.

The input function will be mapped in-place on the current grid of the function, i.e. no further grid refinement.

Parameters:

fmap – [in] mapping function

int mrcpp::FunctionTree::crop(double prec, double splitFac = 1.0, bool absPrec = true)

Reduce the precision of the tree by deleting nodes.

This will run the tree building algorithm in “reverse”, starting from the leaf nodes, and perform split checks on each node based on the given precision and the local wavelet norm.

Note

The splitting factor appears in the threshold for the wavelet norm as \( ||w|| < 2^{-sn/2} ||f|| \epsilon \). In principal, s should be equal to the dimension; in practice, it is set to s=1.

Parameters:

  • prec – New precision criterion

  • splitFac – Splitting factor: 1, 2 or 3

  • absPrec – Use absolute precision

Warning

doxygenfunction: Unable to resolve function “mrcpp::refine_grid” with arguments (FunctionTree<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, double prec, bool absPrec)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, int scales)

Warning

doxygenfunction: Unable to resolve function “mrcpp::refine_grid” with arguments (FunctionTree<D>&, double, bool) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, double prec, bool absPrec)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, int scales)

Warning

doxygenfunction: Unable to resolve function “mrcpp::refine_grid” with arguments (FunctionTree<D>&, FunctionTree<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, FunctionTree<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, const RepresentableFunction<D, T> &inp)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, double prec, bool absPrec)
- template<int D, typename T> int refine_grid(FunctionTree<D, T> &out, int scales)

File I/O

void mrcpp::FunctionTree::saveTree(const std::string &file)

Write the tree structure to disk, for later use.

Parameters:

file – [in] File name, will get “.tree” extension

void mrcpp::FunctionTree::loadTree(const std::string &file)

Read a previously stored tree structure from disk.

Note

This tree must have the exact same MRA the one that was saved

Parameters:

file – [in] File name, will get “.tree” extension

Extracting data

Given a FunctionTree that is a well defined function representation, the following data can be extracted:

T mrcpp::FunctionTree::integrate() const
Returns:

Integral of the function over the entire computational domain

virtual T mrcpp::FunctionTree::evalf(const Coord<D> &r) const override

Note

This will only evaluate the scaling part of the leaf nodes in the tree, which means that the function values will not be fully accurate. This is done to allow a fast and const function evaluation that can be done in OMP parallel. If you want to include also the final wavelet part you can call the corresponding evalf_precise function, or you can manually extend the MW grid by one level before evaluating, using mrcpp::refine_grid(tree, 1)

Parameters:

r – [in] Cartesian coordinate

Returns:

Function value in a point, out of bounds returns zero

inline double mrcpp::MWTree::getSquareNorm() const
Returns:

Squared L2 norm of the function

inline int mrcpp::MWTree::getNNodes() const
Returns:

the total number of nodes in the tree

int mrcpp::MWTree::getSizeNodes() const
Returns:

Size of all MW coefs in the tree, in kB

Warning

doxygenfunction: Unable to resolve function “mrcpp::dot” with arguments “(FunctionTree<D>&, FunctionTree<D>&)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template ComplexDouble dot (CompFunction< 3 > bra, CompFunction< 3 > ket) ———^

FunctionTreeVector

The FunctionTreeVector is simply an alias for a std::vector of std::tuple containing a numerical coefficient and a FunctionTree pointer.

Warning

doxygenfunction: Unable to resolve function “mrcpp::clear” with arguments (FunctionTreeVector<D>&, bool) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> void clear(FunctionTreeVector<D, T> &fs, bool dealloc = false)

Warning

doxygenfunction: Unable to resolve function “mrcpp::get_coef” with arguments (const FunctionTreeVector<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> T get_coef(const FunctionTreeVector<D, T> &fs, int i)

Warning

doxygenfunction: Unable to resolve function “mrcpp::get_func” with arguments (FunctionTreeVector<D>&, int) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> FunctionTree<D, T> &get_func(FunctionTreeVector<D, T> &fs, int i)
- template<int D, typename T> const FunctionTree<D, T> &get_func(const FunctionTreeVector<D, T> &fs, int i)

Warning

doxygenfunction: Unable to resolve function “mrcpp::get_n_nodes” with arguments (const FunctionTreeVector<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> int get_n_nodes(const FunctionTreeVector<D, T> &fs)

Warning

doxygenfunction: Unable to resolve function “mrcpp::get_size_nodes” with arguments (const FunctionTreeVector<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches: 

- template<int D, typename T> int get_size_nodes(const FunctionTreeVector<D, T> &fs)

Examples

Constructing an MRA

An MRA is defined in two steps, first the computational domain is given by a BoundingBox (D is the dimension), e.g. for a total domain of \([-16,16]^3\) in three dimensions (eight root boxes of size \([16]^3\) each):

int n = -4;                                        // Root scale defines box size 2^{-n}
std::array<int, 3> l{-1, -1, -1};                  // Translation of first box [l_x,l_y,l_z]
std::array<int, 3> nb{2, 2, 2};                    // Number of boxes [n_x,n_y,n_z]
mrcpp::BoundingBox<3> world(n, l, nb);

which is combined with a ScalingBasis to give an MRA, e.g. interpolating scaling functions of order \(k=9\):

int N = 20;                                        // Maximum refinement 2^{-(n+N)}
int k = 9;                                         // Polynomial order
mrcpp::InterpolatingBasis basis(k);                // Legendre or Interpolating basis
mrcpp::MultiResolutionAnalysis<D> MRA(world, basis, N);

Two types of ScalingBasis are supported (LegendreBasis and InterpolatingBasis), and they are both available at orders \(k=1,2,\dots,40\) (note that some operators are constructed using intermediates of order \(2k\), so in that case the maximum order available is \(k=20\)).

Working withFunctionTreeVectors

Elements can be appended to the vector using the std::make_tuple, elements are obtained with the get_func and get_coef functions:

mrcpp::FunctionTreeVector<D> tree_vec;             // Initialize empty vector
tree_vec.push_back(std::make_tuple(2.0, &tree_a)); // Push back pointer to FunctionTree
auto coef = mrcpp::get_coef(tree_vec, 0);          // Get coefficient of first entry
auto &tree = mrcpp::get_func(tree_vec, 0);         // Get function of first entry
mrcpp::clear(tree_vec, false);                     // Bool argument for tree destruction

Building empty grids

Sometimes it is useful to construct an empty grid based on some available information of the function that is about to be represented. This can be e.g. that you want to copy the grid of an existing FunctionTree or that an analytic function has more or less known grid requirements (like Gaussians). Sometimes it is even necessary to force the grid refinement beyond the coarsest scales in order for the adaptive refining algorithm to detect a wavelet “signal” that allows it to do its job properly (this happens for narrow Gaussians where none of the initial quadrature points hits a function value significantly different from zero).

The simplest way to build an empty grid is to copy the grid from an existing tree (assume that f_tree has been properly built so that it contains more than just root nodes)

mrcpp::FunctionTree<D> f_tree(MRA);   // Input tree
mrcpp::FunctionTree<D> g_tree(MRA);   // Output tree

mrcpp::project(prec, f_tree, f_func); // Build adaptive grid for f_tree
mrcpp::copy_grid(g_tree, f_tree);     // Copy grid from f_tree to g_tree

Passing an analytic function as argument to the generator will build a grid based on some predefined information of the function (if there is any, otherwise it will do nothing)

mrcpp::RepresentableFunction<D> func; // Analytic function
mrcpp::FunctionTree<D> tree(MRA);     // Output tree
mrcpp::build_grid(tree, func);        // Build grid based on f_func

The lambda analytic functions do not provide such information, this must be explicitly implemented as a RepresentableFunction sub-class (see MRCPP programmer’s guide for details).

Actually, the effect of the build_grid is to extend the existing grid with any missing nodes relative to the input. There is also a version of build_grid taking a FunctionTree argument. Its effect is very similar to the copy_grid above, with the only difference that now the output grid is extended with the missing nodes (e.i. the nodes that are already there are not removed first). This means that we can build the union of two grids by successive applications of build_grid

mrcpp::FunctionTree<D> f_tree(MRA);   // Construct empty grid of root nodes
mrcpp::build_grid(f_tree, g_tree);    // Extend f with missing nodes relative to g
mrcpp::build_grid(f_tree, h_tree);    // Extend f with missing nodes relative to h

In contrast, doing the same with copy_grid would clear the f_tree grid in between, and you would only get a (identical) copy of the last h_tree grid, with no memory of the g_tree grid that was once there. One can also make the grids of two functions equal to their union

mrcpp::build_grid(f_tree, g_tree);    // Extend f with missing nodes relative to g
mrcpp::build_grid(g_tree, f_tree);    // Extend g with missing nodes relative to f

The union grid of several trees can be constructed in one go using a FunctionTreeVector

mrcpp::FunctionTreeVector<D> inp_vec;
inp_vec.push_back(std::make_tuple(1.0, tree_1));
inp_vec.push_back(std::make_tuple(1.0, tree_2));
inp_vec.push_back(std::make_tuple(1.0, tree_3));

mrcpp::FunctionTree<D> f_tree(MRA);
mrcpp::build_grid(f_tree, inp_vec);   // Extend f with missing nodes from all trees in inp_vec

Projection

The project function takes an analytic D-dimensional scalar function (which can be defined as a lambda function or one of the explicitly implemented sub-classes of the RepresentableFunction base class in MRCPP) and projects it with the given precision onto the MRA defined by the FunctionTree. E.g. a unit charge Gaussian is projected in the following way (the MRA must be initialized as above)

// Defining an analytic function
double beta = 10.0;
double alpha = std::pow(beta/pi, 3.0/2.0);
auto func = [alpha, beta] (const mrcpp::Coord<3> &r) -> double {
    double R = std::sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
    return alpha*std::exp(-beta*R*R);
};

double prec = 1.0e-5;
mrcpp::FunctionTree<3> tree(MRA);
mrcpp::project(prec, tree, func);

This projection will start at the default initial grid (only the root nodes of the given MRA), and adaptively build the full grid. Alternatively, the grid can be estimated a priori if the analytical function has some known features, such as for Gaussians:

double prec;                                       // Precision of the projection
int max_iter;                                      // Maximum levels of refinement

mrcpp::GaussFunc<D> func;                          // Analytic Gaussian function
mrcpp::FunctionTree<D> tree(MRA);                  // Output tree

mrcpp::build_grid(tree, func);                     // Empty grid from analytic function
mrcpp::project(prec, tree, func, max_iter);        // Starts projecting from given grid

This will first produce an empty grid suited for representing the analytic function func (this is meant as a way to make sure that the projection starts on a grid where the function is actually visible, as for very narrow Gaussians, it’s not meant to be a good approximation of the final grid) and then perform the projection on the given numerical grid. With a negative prec (or max_iter = 0) the projection will be performed strictly on the given initial grid, with no further refinements.

Addition

Arithmetic operations in the MW representation are performed using the FunctionTreeVector, and the general sum \(f = \sum_i c_i f_i(x)\) is done in the following way

double a, b, c;                                    // Addition parameters
mrcpp::FunctionTree<D> a_tree(MRA);                // Input function
mrcpp::FunctionTree<D> b_tree(MRA);                // Input function
mrcpp::FunctionTree<D> c_tree(MRA);                // Input function

mrcpp::FunctionTreeVector<D> inp_vec;              // Vector to hold input functions
inp_vec.push_back(std::make_tuple(a, &a_tree));    // Append to vector
inp_vec.push_back(std::make_tuple(b, &b_tree));    // Append to vector
inp_vec.push_back(std::make_tuple(c, &c_tree));    // Append to vector

mrcpp::FunctionTree<D> f_tree(MRA);                // Output function
mrcpp::add(prec, f_tree, inp_vec);                 // Adaptive addition

The default initial grid is again only the root nodes, and a positive prec is required to build an adaptive tree structure for the result. The special case of adding two functions can be done directly without initializing a FunctionTreeVector

mrcpp::FunctionTree<D> f_tree(MRA);
mrcpp::add(prec, f_tree, a, a_tree, b, b_tree);

Addition of two functions is usually done on their (fixed) union grid

mrcpp::FunctionTree<D> f_tree(MRA);                // Construct empty root grid
mrcpp::build_grid(f_tree, a_tree);                 // Copy grid of g
mrcpp::build_grid(f_tree, b_tree);                 // Copy grid of h
mrcpp::add(-1.0, f_tree, a, a_tree, b, b_tree);    // Add functions on fixed grid

Note that in the case of addition there is no extra information to be gained by going beyond the finest refinement levels of the input functions, so the union grid summation is simply the best you can do, and adding a positive prec will not make a difference. There are situations where you want to use a smaller grid, though, e.g. when performing a unitary transformation among a set of FunctionTrees. In this case you usually don’t want to construct all the output functions on the union grid of all the input functions, and this can be done by adding the functions adaptively starting from root nodes.

If you have a summation over several functions but want to perform the addition on the grid given by the first input function, you first copy the wanted grid and then perform the operation on that grid

mrcpp::FunctionTreeVector<D> inp_vec;
inp_vec.push_back(std::make_tuple(a, a_tree));
inp_vec.push_back(std::make_tuple(b, b_tree));
inp_vec.push_back(std::make_tuple(c, c_tree));

mrcpp::FunctionTree<D> f_tree(MRA);                // Construct empty root grid
mrcpp::copy_grid(f_tree, get_func(inp_vec, 0));    // Copy grid of first input function
mrcpp::add(-1.0, f_tree, inp_vec);                 // Add functions on fixed grid

Here you can of course also add a positive prec to the addition and the resulting function will be built adaptively starting from the given initial grid.

Multiplication

The multiplication follows the exact same syntax as the addition, where the product \(f = \prod_i c_i f_i(x)\) is done in the following way

double a, b, c;                                    // Multiplication parameters
mrcpp::FunctionTree<D> a_tree(MRA);                // Input function
mrcpp::FunctionTree<D> b_tree(MRA);                // Input function
mrcpp::FunctionTree<D> c_tree(MRA);                // Input function

mrcpp::FunctionTreeVector<D> inp_vec;              // Vector to hold input functions
inp_vec.push_back(std::make_tuple(a, &a_tree));    // Append to vector
inp_vec.push_back(std::make_tuple(b, &b_tree));    // Append to vector
inp_vec.push_back(std::make_tuple(c, &c_tree));    // Append to vector

mrcpp::FunctionTree<D> f_tree(MRA);                // Output function
mrcpp::multipy(prec, f_tree, inp_vec);             // Adaptive multiplication

In the special case of multiplying two functions the coefficients are collected into one argument

mrcpp::FunctionTree<D> f_tree(MRA);
mrcpp::multiply(prec, f_tree, a*b, a_tree, b_tree);

For multiplications, there might be a loss of accuracy if the product is restricted to the union grid. The reason for this is that the product will contain signals of higher frequency than each of the input functions, which require a higher grid refinement for accurate representation. By specifying a positive prec you will allow the grid to adapt to the higher frequencies, but it is usually a good idea to restrict to one extra refinement level beyond the union grid (by setting max_iter=1) as the grids are not guaranteed to converge for such local operations (like arithmetics, derivatives and function mappings)

mrcpp::FunctionTree<D> f_tree(MRA);                    // Construct empty root grid
mrcpp::build_grid(f_tree, a_tree);                     // Copy grid of a
mrcpp::build_grid(f_tree, b_tree);                     // Copy grid of b
mrcpp::multiply(prec, f_tree, a*b, a_tree, b_tree, 1); // Allow 1 extra refinement

Re-using grids

Given a FunctionTree that is a valid function representation, we can clear its MW expansion coefficients as well as its grid refinement

mrcpp::FunctionTree<D> tree(MRA);                      // tree is an undefined function
mrcpp::project(prec, tree, f_func);                    // tree represents analytic function f
tree.clear();                                          // tree is an undefined function
mrcpp::project(prec, tree, f_func);                    // tree represents analytic function g

This action will leave the FunctionTree in the same state as after construction (undefined function, only root nodes), and its coefficients can now be re-computed.

In certain situations it might be desireable to separate the actions of computing MW coefficients and refining the grid. For this we can use the refine_grid, which will adaptively refine the grid one level (based on the wavelet norm and the given precision) and project the existing function representation onto the new finer grid

mrcpp::refine_grid(tree, prec);

E.i., this will not change the function that is represented in tree, but it might increase its grid size. The same effect can be made using another FunctionTree argument instead of the precision parameter

mrcpp::refine_grid(tree_out, tree_in);

which will extend the grid of tree_out in the same way as build_grid as shown above, but it will keep the function representation in tree_out.

This functionality can be combined with clear_grid to make a “manual” adaptive building algorithm. One example where this might be useful is in iterative algorithms where you want to fix the grid size for all calculations within one cycle and then relax the grid in the end in preparation for the next iteration. The following is equivalent to the adaptive projection above (refine_grid returns the number of new nodes that were created in the process)

int n_nodes = 1;
while (n_nodes > 0) {
    mrcpp::project(-1.0, tree, func);                  // Project f on fixed grid
    n_nodes = mrcpp::refine_grid(tree, prec);          // Refine grid based on prec
    if (n_nodes > 0) mrcpp::clear_grid(tree);          // Clear grid for next iteration
}