HierMat package¶
block_cluster_tree module¶
block_cluster_tree.py: BlockClusterTree
,
build_block_cluster_tree()
,
recursion_build_block_cluster_tree()
-
class
HierMat.block_cluster_tree.
BlockClusterTree
(left_clustertree, right_clustertree, sons=None, level=0, is_admissible=False, plot_info=None)[source]¶ Bases:
object
Tree structure containing two
ClusterTree
objects.Parameters: - left_clustertree (ClusterTree) – “left side” cluster tree
- right_clustertree (ClusterTree) – “right side” cluster tree
- sons (list(BlockClusterTree)) – sons if the current tree is not a leaf
- level (int) – level of the current tree
- is_admissible (bool) – whether the current tree is admissible or not
- plot_info (tuple(int, int)) – coordinates of the current tree
-
depth
(root_level=None)[source]¶ Depth of the tree.
The root is at level 0 if not specified otherwise.
Parameters: root_level (int) – internal use. Returns: depth Return type: int Hint
recursive function
-
draw
(axes, admissible_color='#1e26bc', inadmissible_color='#bc1d38')[source]¶ Draw a patch into given axes.
Parameters: - axes (matplotlib.pyplot.axes) – axes instance to draw in
- admissible_color (str) – color for admissible patch (see matplotlib for color specs)
- inadmissible_color (str) – color for inadmissible patch
-
shape
()[source]¶ Return length of left and right cluster tree.
Returns: x and y dimension Return type: tuple(int, int)
-
HierMat.block_cluster_tree.
build_block_cluster_tree
(left_cluster_tree, right_cluster_tree=None, start_level=0, admissible_function=<function admissible>)[source]¶ Factory for building a block cluster tree.
Parameters: - left_cluster_tree (ClusterTree) – “left side” cluster tree
- right_cluster_tree (ClusterTree) – “right side” cluster tree
- start_level (int) – counter that keeps track of the level
- admissible_function (function that returns a bool) – function that determines whether two cluster trees are admissible or not. This is
crucial for the structure of the block cluster tree. (See
utils.admissible
for an example)
Returns: block cluster tree
Return type:
-
HierMat.block_cluster_tree.
recursion_build_block_cluster_tree
(current_tree, admissible_function)[source]¶ Recursion to
build_block_cluster_tree()
.
cluster module¶
cluster.py: Cluster
object and iterator
-
class
HierMat.cluster.
Cluster
(grid, indices=None)[source]¶ Bases:
object
Handles operations on a
Grid
object by manipulating an index list.Parameters: - grid (Grid) – grid to build cluster around
- indices (list) – index list (optional)
-
diameter
()[source]¶ Return diameter.
Return the maximal Euclidean distance between two points. For big Clusters this is costly.
Returns: diameter Return type: float
-
distance
(other)[source]¶ Compute distance to other cluster
Return minimal Euclidean distance between points of self and other
Parameters: other – another instance of Cluster Returns: distance Return type: float
-
get_grid_item_support
(item)[source]¶ Return support of item from grid.
Parameters: item (tuple(float)) – point
-
get_grid_item_support_by_index
(item)[source]¶ Return support of i-th item from the current index list.
This enables the user to get the support of the “last” item by calling get_grid_item_support_by_index(-1)
Parameters: item (int) – index
cluster_tree module¶
cluster_tree.py: ClusterTree
, build_cluster_tree()
, recursion_build_cluster_tree()
-
class
HierMat.cluster_tree.
ClusterTree
(content, sons=None, max_leaf_size=1, level=0)[source]¶ Bases:
object
Tree structure built according to the
Splitable
in useSplits are done until max_leaf_size is reached
-
depth
(root_level=None)[source]¶ Get depth of the tree
Parameters: root_level (int) – internal use. Returns: depth Return type: int Hint
recursive function
-
distance
(other)[source]¶ Return distance to other
Parameters: other (ClusterTree) – other cluster tree Returns: distance Return type: float
-
get_grid_item
(item)[source]¶ Get grid item from content
Parameters: item (int) – index of item to get
-
get_grid_item_support
(item)[source]¶ Return supports ofitem from grid
Parameters: item (tuple(float)) – point
-
get_grid_item_support_by_index
(item)[source]¶ Return supports of i-th item from grid
Parameters: item (int) – index
-
get_index
(item)[source]¶ Get index from content
Parameters: item (int) – index to get Returns: index Return type: int
-
get_patch_coordinates
()[source]¶ Return min and max out of indices
Returns: min and max Return type: tuple(int, int)
-
-
HierMat.cluster_tree.
admissible
(left_clustertree, right_clustertree)[source]¶ Default admissible condition for BlockClusterTree
True if the smaller diameter of the input is smaller or equal to the distance between the two ClusterTrees
Parameters: - left_clustertree (ClusterTree) – “Left-side” ClusterTree
- right_clustertree (ClusterTree) – “Right-side” ClusterTree
Returns: admissible
Return type:
-
HierMat.cluster_tree.
build_cluster_tree
(splitable, max_leaf_size=1, start_level=0)[source]¶ Factory for building a cluster tree
Parameters: Returns: cluster tree
Return type:
-
HierMat.cluster_tree.
recursion_build_cluster_tree
(current_tree)[source]¶ Recursion to
build_cluster_tree()
cuboid module¶
cuboid.py: Cuboid
and minimal_cuboid()
-
class
HierMat.cuboid.
Cuboid
(low_corner, high_corner)[source]¶ Bases:
object
Cuboid that is parallel to the axis in Cartesian coordinates.
Characterized by two diagonal corners.
Attributes:
low_corner: numpy.array with minimal values in each dimension
high_corner: numpy.array with maximal values in each dimension
-
diameter
()[source]¶ Return the diameter of the cuboid.
Diameter is the Euclidean distance between low_corner and high_corner.
Returns: diameter Return type: float
grid module¶
grid.py: Grid
object and iterator
-
class
HierMat.grid.
Grid
(points, supports)[source]¶ Bases:
object
Discretized grid characterized by points and supports
Parameters: Raises: ValueError – if points and supports have different length
-
dim
()[source]¶ Dimension of the Grid
Note
this is the dimension of the first point
if you store points of different dimensions, this will be misleadingReturns: dimension Return type: int
-
hmat module¶
hmat.py: HMat
, build_hmatrix()
, recursion_build_hmatrix()
, StructureWarning
-
class
HierMat.hmat.
HMat
(blocks=(), content=None, shape=(), parent_index=())[source]¶ Bases:
object
Implement a hierarchical Matrix
Parameters: - blocks (list(HMat)) – list of HMat instances (children) (optional)
- content (RMat or numpy.matrix) – the content if the matrix has no children (optional)
- shape (tuple(int, int)) – the shape of the matrix (same as for numpy matrices)
- parent_index (tuple(int, int)) – the index of this matrix with respect to its containing parent-matrix, zero based
Supported operations
+
(formatted addition)*
(formatted multiplication)-a
(unary minus)-
(formatted subtraction)==
(equal)!=
(not equal)abs
(frobenius norm)inv
(inversion)
-
block_structure
()[source]¶ return the block structure of self
produce a dict with parent_index: block paris
Return type: dict Returns: dict with index: HMatrix pairs
-
column_sequence
()[source]¶ Return the sequence of column groups as in thm. 1.9.4. in [Eve66]
only for consistent matrices
Returns: list of column-sizes in first row Return type: list(int) Raises: StructureWarning
if the matrix is not consistent
-
inv
(make_copy=True)[source]¶ Invert the matrix according to chapter 7.5.1 in [Hac15]
Note
If make_copy is True, a deepcopy is made at the beginning. This can be very slow for large matrices.
Parameters: make_copy (bool) – make a copy of the matrix to not change the original (Default True) Returns: approximation of inverse Return type: HMat Raises: numpy.LinAlgError
if matrix is singular
-
is_consistent
()[source]¶ check if the blocks are aligned, i.e. we have consistent rows and columns as in def. 1.9.5 in [Eve66]
Returns: True on consistency, False otherwise Return type: bool
-
norm
(order=None)[source]¶ Norm of the matrix
Parameters: order – order of the norm (see in numpy.linalg.norm()
)Returns: norm Return type: float
-
restructure
(structure)[source]¶ Restructure self into blocks to match structure
Parameters: structure (dict(index: size)) – a dict as returned by HMat.block_structure Returns: HMat that has the specified structure Return type: HMat Raises: NotImplementedError
if self has blocks or if self contains unknown content
-
row_sequence
()[source]¶ Return the sequence of row groups as in thm. 1.9.4. in [Eve66]
defined only for consistent matrices
Returns: list of row-sizes in first column Return type: list(int) Raises: StructureWarning
if the matrix is not consistent
-
split_hmat
(start_x, start_y, end_x, end_y)[source]¶ Fetch the block specified by indices from a numpy.matrix
Parameters: Returns: numpy.matrix with corresponding block
Return type:
-
split_rmat
(start_x, start_y, end_x, end_y)[source]¶ Fetch the block specified by indices from a rmat
Example
r = RMat(numpy.matrix([[2], [2], [2]]), numpy.matrix([[3], [3], [3]])) h = HMat(content=r, shape=(3, 3), parent_index=(0, 0)) res = h.split_rmat(0, 2, 1, 3) res <RMat with left_mat: matrix([[2]]), right_mat: matrix([[3]]) and max_rank: None>
Parameters: Returns: rmat with corresponding block
Return type:
-
to_matrix
()[source]¶ Full matrix representation
Returns: full matrix Return type: numpy.matrix
-
zero
(make_copy=True)[source]¶ Return a copy of self with zero entries
i.e. an HMat instance with the same structure, but all content is replaced by numpy.zeros
Note
If make_copy is True, a deepcopy is made at the beginning. This can be very slow for large matrices.
Parameters: make_copy (bool) – make a copy of the matrix to not change the original (Default True) Returns: zero matrix Return type: HMat
-
exception
HierMat.hmat.
StructureWarning
[source]¶ Bases:
exceptions.Warning
Special Warning used in this module to indicate bad block structure
-
HierMat.hmat.
build_hmatrix
(block_cluster_tree=None, generate_rmat_function=None, generate_full_matrix_function=None)[source]¶ Factory to build an hierarchical matrix
Takes a block cluster tree and generating functions for full matrices and low rank matrices respectively
The generating functions take a block cluster tree as input and return a RMat or numpy.matrix respectively
Parameters: - block_cluster_tree (BlockClusterTree) – block cluster tree giving the structure
- generate_rmat_function – function taking an admissible block cluster tree and returning a rank-k matrix
- generate_full_matrix_function – function taking an inadmissible block cluster tree and returning a numpy.matrix
Returns: hmatrix
Return type:
-
HierMat.hmat.
recursion_build_hmatrix
(current_hmat, block_cluster_tree, generate_rmat, generate_full_mat)[source]¶ Recursion to
build_hmatrix()
rmat module¶
rmat.py: RMat
-
class
HierMat.rmat.
RMat
(left_mat, right_mat=None, max_rank=None)[source]¶ Bases:
object
Rank-k matrix
Implementation of the low rank matrix as described in [HB2015]
-
form_add
(other, rank=None)[source]¶ Formatted addition of self and other, i.e. addition and reduction to rank:
(self + other).reduce(rank)
If rank is omitted, reduction to min(rank(self), rank(other))
Parameters: Returns: reduced result
Return type:
-
static
from_matrix
(matrix, max_rank)[source]¶ Convert a full matrix to the rank k format
Parameters: - matrix (numpy.matrix) – Matrix to convert to RMat
- max_rank (int) – maximum rank
Returns: Best approximation of matrix in RMat format
Return type:
-
norm
(order=None)[source]¶ Norm of the matrix
Parameters: order – order of the norm (see in numpy.linalg.norm()
)Returns: norm Return type: float
-
reduce
(new_k)[source]¶ Perform a reduced QR decomposition to rank new_k
Parameters: new_k (int) – rank to reduce to Returns: reduced matrix Return type: RMat
-
split
(start_x, start_y, end_x, end_y)[source]¶ Return the block of self that matches the supplied indices
Parameters: Returns: rmat that represents the corresponding block
Return type:
-
to_matrix
()[source]¶ Full matrix representation:
left_mat * right_mat.T
Returns: full matrix Return type: numpy.matrix
-
splitable module¶
splitable.py: Splitable
(Interface) and iterator and different strategies
-
Starting from a minimal cuboid that get’s split in half on every level, points are distributed according to their geometric location
-
Same as RegularCuboid, but after every split, the cuboids are reduced to minimal size
-
Distributes the points evenly on every level, depends solely on the order of the initial list
-
class
HierMat.splitable.
Balanced
(cluster)[source]¶ Bases:
HierMat.splitable.Splitable
Balanced strategy, where a cluster is always restructure in half only according to its size, without regarding the geometry
Gives a binary tree
-
diameter
()[source]¶ Return the diameter
Diameter of the contained cluster
Returns: diameter Return type: float
-
distance
(other)[source]¶ Distance to other balanced
Parameters: other (Balanced) – other balanced Returns: distance Return type: float
-
get_grid_item
(item)[source]¶ Get grid item from cluster
Parameters: item (int) – index of item to get
-
get_grid_item_support
(item)[source]¶ Return supports of item from grid
Parameters: item (tuple(float)) – point
-
get_grid_item_support_by_index
(item)[source]¶ Return supports of item from grid
Parameters: item (int) – index
-
get_index
(item)[source]¶ Get index from cluster
Parameters: item (int) – index to get Returns: index Return type: int
-
-
class
HierMat.splitable.
MinimalCuboid
(cluster)[source]¶ Bases:
HierMat.splitable.RegularCuboid
Method of minimal cuboids
Split is implemented by splitting the surrounding cuboid in half along longest axis and distributing indices according to the cuboid they belong to. Afterwards all resulting cuboids are shrunk to minimal.
Gives a binary tree
-
split
()[source]¶ Split the cuboid and distribute items in cluster according to the cuboid they belong to. Reduce every restructure to minimal size
Returns: list of minimal cuboids of smaller size Return type: list(MinimalCuboid)
-
-
class
HierMat.splitable.
RegularCuboid
(cluster, cuboid=None)[source]¶ Bases:
HierMat.splitable.Splitable
Method of regular cuboids
If no cuboid is given, a minimal cuboid is built around initial list of indices. Split is then implemented by splitting the surrounding cuboid in half along longest axis and distributing indices according to the cuboid they belong to.
Gives a binary tree
-
diameter
()[source]¶ Return the diameter of the cuboid
Diameter of the surrounding cuboid
Returns: diameter Return type: float
-
distance
(other)[source]¶ Return the distance between own cuboid and other cuboid
Parameters: other (RegularCuboid) – other regular cuboid Returns: distance Return type: float
-
split
()[source]¶ Split the cuboid and distribute items in cluster according to the cuboid they belong to
Returns: list of smaller regular cuboids Return type: list(RegularCuboid)
-
-
class
HierMat.splitable.
Splitable
[source]¶ Bases:
object
Interface to the different strategies that can be used to restructure a cluster in two or more.
Methods that need to be implemented by subclasses
- __eq__: check for equality against other Splitable
- split: restructure the object in two or more, return new instances
- diameter: return the diameter of the object
- distance: return the distance to other object
-
cluster
= None¶
-
get_grid_item
(item)[source]¶ Get grid item from cluster
Parameters: item (int) – index of item to get
-
get_grid_item_support
(item)[source]¶ Return supports of item from grid
Parameters: item (tuple(float)) – point
-
get_grid_item_support_by_index
(item)[source]¶ Return supports of item from grid
Parameters: item (int) – index
-
get_index
(item)[source]¶ Get index from cluster
Parameters: item (int) – index to get Returns: index Return type: int
utils module¶
utils.py: Utilities for the HMatrix
module
-
HierMat.utils.
block_cluster_tree_plot
(obj, filename=None, ticks=False, face_color='#133f52', admissible_color='#76f7a8', inadmissible_color='#ff234b')[source]¶ Plot the block cluster tree
Parameters: - obj (BlockClusterTree) – block cluster tree to plot
- filename (str) – filename to save the plot. if omitted, the plot will be displayed
- ticks (bool) – show ticks in the plot
- face_color – background color (see matplotlib for color specs)
- admissible_color (str) – color for admissible patch
- inadmissible_color (str) – color for inadmissible patch
Note
depends on
matplotlib.pyplot
-
HierMat.utils.
divisor_generator
(n)[source]¶ Return divisors of n
Parameters: n (int) – integer to find divisors of Returns: divisors Return type: list[int] Warning
This is a generator! To get a list with all divisors call:
list(divisor_generator(n))
Note
found at StackOverflow on 2017.03.08
-
HierMat.utils.
export
(obj, form='xml', out_file='out')[source]¶ Export obj in specified format.
Parameters: - obj (BlockClusterTree or ClusterTree or HMat or RMat) – object to export
- form (str) – format specifier
- out_file (str) – path to output file
Raises: NotImplementedError – if form is not supported
Note
implemented so far:
- xml
- dot
- bin
-
HierMat.utils.
load
(filename)[source]¶ Load a
ClusterTree
orBlockClusterTree
from fileParameters: filename (String) – file to import Returns: object Return type: BlockClusterTree or ClusterTree Note
Depends on
pickle
-
HierMat.utils.
plot
(obj, filename=None, **kwargs)[source]¶ plot an object
Parameters: - obj (BlockClusterTree) – object to plot
- filename (str) – filename to save the plot to (if omitted, the plot will be displayed)
- kwargs – optional arguments to specific plot commands see the respective documentations
Examples¶
Examples to the HierMat package
In this module are example scripts to demonstrate the usage of the package
model_1d module¶
This is an implementation of the one-dimensional model example described in [thesis]. It shows a typical use-case of hierarchical matrices:
Integral Equation
rac{1}{2pi}log |x-y| u(y) dy = g
for \(x \in [0,1]\).
After Galerkin discretization, we end up with a linear system
\[\mathbf{A}^{\text{Gal}} \cdot \alpha = \mathbf{g}\]where
\[\mathbf{A}^{\text{Gal}}_{t,\tau}:= \int_{t}\int_{\tau} \log |x-y| dy dx\]To determine admissible blocks we start by building the geometric objects:
midpoints = [((i + 0.5)/n,) for i in xrange(n)] intervals = {p: (p[0] - 0.5/n, p[0] + 0.5/n) for p in midpoints} grid = HierMat.Grid(points=midpoints, supports=intervals) cluster = HierMat.Cluster(grid=grid) unit_cuboid = HierMat.Cuboid([0], [1]) strategy = HierMat.RegularCuboid(cluster=cluster, cuboid=unit_cuboid) cluster_tree = HierMat.build_cluster_tree(splitable=strategy, max_leaf_size=n_min) block_cluster_tree = HierMat.build_block_cluster_tree(left_cluster_tree=cluster_tree, right_cluster_tree=cluster_tree, admissible_function=HierMat.admissible )With the structure established, we can produce the hierarchical matrix:
hmat = HierMat.build_hmatrix(block_cluster_tree=block_cluster_tree, generate_rmat_function=lambda bct: galerkin_1d_rank_k(bct, max_rank), generate_full_matrix_function=galerkin_1d_full )
-
HierMat.examples.model_1d.
galerkin_1d_full
(block_cluster_tree)[source]¶ Exact calculation of the integral
\[A_{i,j}=A_{\tau,t}^{Gal}=\int_t\int_\tau\log\Vert x-y\Vert \;dydx\]Parameters: block_cluster_tree (HierMat.BlockClusterTree) – inadmissible block cluster tree Returns: matrix with same shape as block_cluster_tree.shape() Return type: numpy.matrix
-
HierMat.examples.model_1d.
galerkin_1d_rank_k
(block_cluster_tree, max_rank)[source]¶ Low-rank approximation of the kernel
\[ \begin{align}\begin{aligned}R: &= A \cdot B^T\\A_{\tau, k}: &= \int_\tau \log \Vert x-y_k\Vert dx\\B_{\tau, k}: &= \int_\tau \mathcal{L}_k(y) dy\end{aligned}\end{align} \]Parameters: - block_cluster_tree (HierMat.BlockClusterTree) – admissible block cluster tree
- max_rank (int) – separation rank
Returns:
-
HierMat.examples.model_1d.
get_chebyshev_interpol_points
(points, lower=0, upper=1)[source]¶ Get Chebyshev interpolation points on interval \((a, b)\).
\[\frac{1}{2} (a + b) + \frac{1}{2} (b - a) \cos\left(\frac{2k+1}{2n}\pi\right)\]for \(k=0, ..., \text{points} - 1\)
Parameters: Returns: \(\frac{1}{2} (a + b) + \frac{1}{2} (b - a) \cos\left(\frac{2k+1}{2n}\pi\right)\)
Return type: list(floats)
-
HierMat.examples.model_1d.
ker
(x, y)[source]¶ Kernel to integrate
\[ker(x, y):= \log\left(\Vert x - y \Vert\right)\]Parameters: Returns: \(\log\left(\Vert x - y \Vert\right)\)
Return type:
-
HierMat.examples.model_1d.
kerlog
(x)[source]¶ kerlog function as in [thesis].
\[kerlog(x):= x^2 \left( \log(\Vert x \Vert) - \frac{1}{2} \right)\]Parameters: x (float) – real number Returns: \(x^2 ( \log(\Vert x \Vert) - \frac{1}{2} )\) Return type: float