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
>
classmrcpp
::
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
(const BoundingBox<D> &bb, const ScalingBasis &sb, int depth = MaxDepth)¶ - Return
New MultiResolutionAnalysis object
- Parameters
[in] bb
: Computational domain[in] sb
: Polynomial basis[in] depth
: Maximum allowed resolution depth, relative to root scale
-
-
template<int
D
>
classmrcpp
::
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 >
Public Functions
-
BoundingBox
(int n = 0, const std::array<int, D> &l = {}, const std::array<int, D> &nb = {}, const std::array<double, D> &sf = {})¶ - Return
New BoundingBox object
- Parameters
[in] n
: Length scale, default 0[in] l
: Corner translation, default [0, 0, …][in] nb
: Number of boxes, default [1, 1, …][in] sf
: Scaling factor, default [1.0, 1.0, …]
-
-
class
mrcpp
::
LegendreBasis
: public mrcpp::ScalingBasis¶ Legendre scaling functions as defined by Alpert, SIAM J Math Anal 24 (1), 246 (1993).
Public Functions
-
LegendreBasis
(int k)¶ - Return
New LegendreBasis object
- Parameters
[in] k
: Polynomial order of basis,1 < k < 40
-
-
class
mrcpp
::
InterpolatingBasis
: public mrcpp::ScalingBasis¶ Interpolating scaling functions as defined by Alpert etal, J Comp Phys 182, 149-190 (2002).
Public Functions
-
InterpolatingBasis
(int k)¶ - Return
New InterpolatingBasis object
- Parameters
[in] k
: Polynomial order of basis,1 < k < 40
-
FunctionTree¶
-
template<int
D
>
classmrcpp
::
FunctionTree
: public mrcpp::MWTree<D>, public mrcpp::RepresentableFunction<D>¶ 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
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.
- Return
New FunctionTree object
- Parameters
[in] mra
: Which MRA the function is defined[in] sh_mem
: Pointer to MPI shared memory block
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:
Start with an initial guess for the grid
Compute the MW coefficients for the output function on the current grid
Refine the grid where necessary based on the local wavelet norm
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, fixed grid.
Keeps the node structure of the tree, even though the zero function is representable at depth zero. Use cropTree to remove unnecessary nodes.
-
template<int
D
>
voidmrcpp
::
project
(double prec, FunctionTree<D> &out, RepresentableFunction<D> &inp, int maxIter, bool absPrec)¶ Project an analytic function onto the MW basis, adaptive grid.
The output function will be computed using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Input function[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
copy_func
(FunctionTree<D> &out, FunctionTree<D> &inp)¶ Copy function from one tree onto the grid of another tree, fixed grid.
The output function will be computed using the general algorithm:
Loop through current leaf nodes of the output tree
Copy MW coefs from the corresponding input node
- Parameters
[out] out
: Output function[in] inp
: Input function
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called and will overwrite any existing coefs.
-
template<int
D
>
voidmrcpp
::
add
(double prec, FunctionTree<D> &out, FunctionTreeVector<D> &inp, int maxIter, bool absPrec)¶ Addition of several MW function representations, adaptive grid.
The output function will be computed as the sum of all input functions in the vector (including their numerical coefficients), using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Vector of input function[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
add
(double prec, FunctionTree<D> &out, double a, FunctionTree<D> &inp_a, double b, FunctionTree<D> &inp_b, int maxIter, bool absPrec)¶ Addition of two MW function representations, adaptive grid.
The output function will be computed as the sum of the two input functions (including the numerical coefficient), using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] a
: Numerical coefficient of function a[in] inp_a
: Input function a[in] b
: Numerical coefficient of function b[in] inp_b
: Input function b[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
multiply
(double prec, FunctionTree<D> &out, FunctionTreeVector<D> &inp, int maxIter, bool absPrec, bool useMaxNorms)¶ Multiplication of several MW function representations, adaptive grid.
The output function will be computed as the product of all input functions in the vector (including their numerical coefficients), using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Vector of input function[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision[in] useMaxNorms
: Build output tree based on norm estimates from input
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
multiply
(double prec, FunctionTree<D> &out, double c, FunctionTree<D> &inp_a, FunctionTree<D> &inp_b, int maxIter, bool absPrec, bool useMaxNorms)¶ Multiplication of two MW function representations, adaptive grid.
The output function will be computed as the product of the two input functions (including the numerical coefficient), using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] c
: Numerical coefficient[in] inp_a
: Input function a[in] inp_b
: Input function b[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision[in] useMaxNorms
: Build output tree based on norm estimates from input
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
square
(double prec, FunctionTree<D> &out, FunctionTree<D> &inp, int maxIter, bool absPrec)¶ Out-of-place square of MW function representations, adaptive grid.
The output function will be computed as the square of the input function, using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Input function to square[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
power
(double prec, FunctionTree<D> &out, FunctionTree<D> &inp, double p, int maxIter, bool absPrec)¶ Out-of-place power of MW function representations, adaptive grid.
The output function will be computed as the input function raised to the given power, using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Input function to square[in] p
: Numerical power[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
-
template<int
D
>
voidmrcpp
::
dot
(double prec, FunctionTree<D> &out, FunctionTreeVector<D> &inp_a, FunctionTreeVector<D> &inp_b, int maxIter, bool absPrec)¶ Dot product of two MW function vectors, adaptive grid.
The output function will be computed as the dot product of the two input vectors (including their numerical coefficients). The precision parameter is used only in the multiplication part, the final addition will be on the fixed union grid of the components.
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp_a
: Input function vector[in] inp_b
: Input function vector[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
- Note
The length of the input vectors must be the same.
-
template<int
D
>
voidmrcpp
::
map
(double prec, FunctionTree<D> &out, FunctionTree<D> &inp, FMap fmap, int maxIter, bool absPrec)¶ map a MW function onto another representations, adaptive grid
The output function tree will be computed by mapping the input tree values through the fmap function, using the general algorithm:
Compute MW coefs on current grid
Refine grid where necessary based on
prec
Repeat until convergence or
maxIter
is reachedprec < 0
ormaxIter = 0
means NO refinementmaxIter < 0
means no bound
- Parameters
[in] prec
: Build precision of output function[out] out
: Output function to be built[in] inp
: Input function[in] fmap
: mapping function[in] maxIter
: Maximum number of refinement iterations in output tree[in] absPrec
: Build output tree based on absolute precision
No assumption is made for how the mapping function looks. It is left to the end-user to guarantee that the mapping function does not lead to numerically unstable/inaccurate situations (e.g. divide by zero, overflow, etc…)
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called (this grid should however be EMPTY, e.i. no coefs).
Creating undefined FunctionTrees¶
The grid of a FunctionTree
can also be constructed without computing any
MW coefficients:
-
template<int
D
>
voidmrcpp
::
build_grid
(FunctionTree<D> &out, const RepresentableFunction<D> &inp, int maxIter)¶ Build empty grid based on info from analytic function.
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 reachedmaxIter < 0
means no bound
- Parameters
[out] out
: Output tree to be built[in] inp
: Input function[in] maxIter
: Maximum number of refinement iterations in output tree
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called. It requires that the functionsisVisibleAtScale()
andisZeroOnInterval()
is implemented in the particularRepresentableFunction
.
-
template<int
D
>
voidmrcpp
::
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 reachedmaxIter < 0
means no bound
- Parameters
[out] out
: Output tree to be built[in] inp
: Input Gaussian expansion[in] maxIter
: Maximum number of refinement iterations in output tree
- 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.
-
template<int
D
>
voidmrcpp
::
build_grid
(FunctionTree<D> &out, FunctionTree<D> &inp, int maxIter)¶ Build empty grid based on another MW function representation.
The grid of the output function will be EXTENDED with all existing nodes in corresponding input function, using the general algorithm:
Loop through current leaf nodes of the output tree
Refine node if the corresponding node in the input has children
Repeat until all input nodes are covered or
maxIter
is reachedmaxIter < 0
means no bound
- Parameters
[out] out
: Output tree to be built[in] inp
: Input tree[in] maxIter
: Maximum number of refinement iterations in output tree
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called. This means that all nodes on the input tree will also be in the final output tree (unlessmaxIter
is reached, but NOT vice versa.
-
template<int
D
>
voidmrcpp
::
build_grid
(FunctionTree<D> &out, FunctionTreeVector<D> &inp, int maxIter)¶ Build empty grid based on several MW function representation.
The grid of the output function will be EXTENDED with all existing nodes in all corresponding input functions, using the general algorithm:
Loop through current leaf nodes of the output tree
Refine node if the corresponding node in one of the inputs has children
Repeat until all input nodes are covered or
maxIter
is reachedmaxIter < 0
means no bound
- Parameters
[out] out
: Output tree to be built[in] inp
: Input tree vector[in] maxIter
: Maximum number of refinement iterations in output tree
- Note
This algorithm will start at whatever grid is present in the
out
tree when the function is called. This means that the final output grid will contain (at least) the union of the nodes of all input trees (unlessmaxIter
is reached).
-
template<int
D
>
voidmrcpp
::
copy_grid
(FunctionTree<D> &out, FunctionTree<D> &inp)¶ Build empty grid that is identical to another MW grid.
- Note
The difference from the corresponding
build_grid
function is that this will first clear the grid of theout
function, whilebuild_grid
will extend the existing grid.- Parameters
[out] out
: Output tree to be built[in] inp
: Input tree
-
template<int
D
>
voidmrcpp
::
clear_grid
(FunctionTree<D> &out)¶ Clear the MW coefficients of a function representation.
- Note
This will only clear the MW coefs in the existing nodes, it will not change the grid of the function. Use
FunctionTree::clear()
to remove all grid refinement as well.- Parameters
[inout] out
: Output function to be cleared
-
void
mrcpp::FunctionTree
::
clear
()¶ Remove all nodes in the tree.
Leaves the tree inn the same state as after construction, i.e. undefined function containing only root nodes without coefficients. The assigned memory (nodeChunks in SerialTree) 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
(double 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
[in] c
: Scalar coefficient
-
void
mrcpp::FunctionTree
::
normalize
()¶ In-place rescaling by a function norm \( ||f||^{-1} \), fixed grid.
-
void
mrcpp::FunctionTree
::
add
(double c, FunctionTree<D> &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
[in] c
: Numerical coefficient of input function[in] inp
: Input function to add
-
void
mrcpp::FunctionTree
::
multiply
(double c, FunctionTree<D> &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
[in] c
: Numerical coefficient of input function[in] inp
: 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
[in] p
: Numerical power
-
void
mrcpp::FunctionTree
::
map
(FMap 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
[in] fmap
: 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.
- Parameters
prec
: New precision criterionsplitFac
: Splitting factor: 1, 2 or 3absPrec
: Use absolute precision
- 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 tos=1
.
-
template<int
D
>
intmrcpp
::
refine_grid
(FunctionTree<D> &out, int scales)¶ Refine the grid of a MW function representation.
This will split ALL leaf nodes in the tree the given number of times, then it will compute scaling coefs of the new nodes, thus leaving the function representation unchanged, but on a larger grid.
- Return
The number of nodes that were split
- Parameters
[inout] out
: Output tree to be refined[in] scales
: Number of refinement levels
-
template<int
D
>
intmrcpp
::
refine_grid
(FunctionTree<D> &out, double prec, bool absPrec)¶ Refine the grid of a MW function representation.
This will first perform a split check on the existing leaf nodes in the tree based on the provided precision parameter, then it will compute scaling coefs of the new nodes, thus leaving the function representation unchanged, but (possibly) on a larger grid.
- Return
The number of nodes that were split
- Parameters
[inout] out
: Output tree to be refined[in] prec
: Precision for initial split check[in] absPrec
: Build output tree based on absolute precision
-
template<int
D
>
intmrcpp
::
refine_grid
(FunctionTree<D> &out, FunctionTree<D> &inp)¶ Refine the grid of a MW function representation.
This will first perform a split check on the existing leaf nodes in the output tree based on the structure of the input tree (same as build_grid), then it will compute scaling coefs of the new nodes, thus leaving the function representation unchanged, but on a larger grid.
- Return
The number of nodes that were split
- Parameters
[inout] out
: Output tree to be refined[in] inp
: Input tree that defines the new grid
File I/O¶
-
void
mrcpp::FunctionTree
::
saveTree
(const std::string &file) override¶ Write the tree structure to disk, for later use.
- Parameters
[in] file
: File name, will get “.tree” extension
-
void
mrcpp::FunctionTree
::
loadTree
(const std::string &file) override¶ Read a previously stored tree structure from disk.
- Note
This tree must have the exact same MRA the one that was saved
- Parameters
[in] file
: File name, will get “.tree” extension
Extracting data¶
Given a FunctionTree
that is a well defined function representation, the
following data can be extracted:
-
double
mrcpp::FunctionTree
::
integrate
() const¶ - Return
Integral of the function over the entire computational domain
-
double
mrcpp::FunctionTree
::
evalf
(const Coord<D> &r) const¶ - Return
Function value in a point, out of bounds returns zero
- 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. If you want to include also the final wavelet part you’ll have to manually extend the MW grid by one level before evaluating, using
mrcpp::refine_grid(tree, 1)
This is done to allow a fast and const function evaluation that can be done in OMP parallel.- Parameters
[in] r
: Cartesian coordinate
-
double
mrcpp::MWTree
::
getSquareNorm
() const¶ - Return
Squared L2 norm of the function
-
int
mrcpp::MWTree
::
getNNodes
(int depth = -1) const¶ - Return
Total number of nodes in the tree, at given depth
- Parameters
[in] depth
: Tree depth to count, negative means count all nodes
-
int
mrcpp::MWTree
::
getSizeNodes
() const¶ - Return
Size of all MW coefs in the tree, in kB
-
template<int
D
>
doublemrcpp
::
dot
(FunctionTree<D> &bra, FunctionTree<D> &ket)¶ The dot product is computed with the trees in compressed form, i.e. scaling coefs only on root nodes, wavelet coefs on all nodes. Since wavelet functions are orthonormal through ALL scales and the root scaling functions are orthonormal to all finer level wavelet functions, this becomes a rather efficient procedure as you only need to compute the dot product where the grids overlap.
- Return
Dot product <bra|ket> of two MW function representations
- Parameters
[in] bra
: Bra side input function[in] ket
: Ket side input function
FunctionTreeVector¶
The FunctionTreeVector
is simply an alias for a std::vector
of
std::tuple
containing a numerical coefficient and a FunctionTree
pointer.
-
template<int
D
>
voidmrcpp
::
clear
(FunctionTreeVector<D> &fs, bool dealloc = false)¶ Remove all entries in the vector.
- Parameters
[in] fs
: Vector to clear[in] dealloc
: Option to free FunctionTree pointer before clearing
-
template<int
D
>
doublemrcpp
::
get_coef
(const FunctionTreeVector<D> &fs, int i)¶ - Return
Numerical coefficient at given position in vector
- Parameters
[in] fs
: Vector to fetch from[in] i
: Position in vector
-
template<int
D
>
FunctionTree<D> &mrcpp
::
get_func
(FunctionTreeVector<D> &fs, int i)¶ - Return
FunctionTree at given position in vector
- Parameters
[in] fs
: Vector to fetch from[in] i
: Position in vector
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
\([-32,32]^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
}