Using H2Pack in Python

We provide a simple python interface for H2Pack called "PyH2Pack" (python3 only).

Installation

First, make sure that H2Pack is installed (See Installing H2Pack) and verify you have libH2Pack.a in the src directory.

In the H2Pack/pyh2pack directory,

• If you are using ICC+MKL: Run LDSHARED="icc -shared" CC=icc python3 setup_icc.py install. If MKL is not installed for python (e.g., via pip or conda), before starting python3, you need to preload the MKL library by export LD_PRELOAD=$MKLROOT/lib/intel64/libmkl_rt.so where $MKLROOT is the path to the installed MKL.

• If you are using GCC+OpenBLAS: you need to modify OPENBLAS_INSTALL_DIR = <path to your OpenBLAS install directory> in setup.py. Run CC=gcc python3 setup.py install to perform the installation.

Usage

Here is a basic test code:

import pyh2pack
import numpy as np
 
#  generate random points in a 3D box
N = 80000
pt_dim = 3
coord = np.random.uniform(0, 1, size=(pt_dim, N))

#  dimension and parameters of the kernel function
krnl_dim = 1
krnl_param = np.array([1., -2.]) # please use double format but not integers
 
#  build the H2 matrix
A = pyh2pack.H2Mat(kernel="Quadratic_3D", krnl_dim=krnl_dim, pt_coord=coord, pt_dim=3, JIT_mode=1, rel_tol=1e-3, krnl_param=krnl_param)

#  H2 matvec
x = np.random.normal(size=(krnl_dim*N))
y = A.matvec(x)

#  direct matvec for a subset of rows
start_pt = 8000
end_pt = 10000
z = A.direct_matvec(x, start_pt, end_pt)

#  error of H2 matvec 
print(np.linalg.norm(y[(start_pt-1)*krnl_dim:end_pt*krnl_dim] - z) / np.linalg.norm(z))

#  print out statistic info of this h2 matrix representation
A.print_statistic()

#. destroy the h2 data structure (optional)
A.clean()

Most of the above functions are straightforward to understand. The detailed input description for the initializer pyh2pack.H2Mat() is as follows:

• kernel : string, kernel name (see below for supported kernel functions);
• krnl_dim : integer, dimension of the kernel function's output;
• pt_coord : 2d numpy array, point coordinates. Each row or column stores the coordinate of one point;
• pt_dim : integer, dimensions of the points (1, 2, or 3);
• rel_tol : float, accuracy threshold for the H2 matrix representation;
• (optional) JIT_mode : 1 or 0, flag for running matvec in JIT mode (JIT mode reduces storage cost but has slower matvec);
• (optional) krnl_param : 1d numpy array, parameters of the kernel function;
• (optional) max_leaf_points: integer, the maximum number of points in each leaf node.;
• (optional) max_leaf_size : double, the maximum edge length of each leaf box.;
• (optional) pp_filename : string, path to a file that either contains precomputed proxy points or will be written to in the setup process;

Supported kernel functions and their description

To list all the available kernel functions provided, run pyh2pack.print_kernels(). The list contains the name, parameter description, default parameter values, and other information about each supported kernel. Below are the currently supported kernels:

==================== Supported kernel functions ====================
 * KERNEL 0: Gaussian_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = exp(-l|x-y|^2) 
    - default proxy point select   : QR
 * KERNEL 1: Gaussian_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = exp(-l|x-y|^2) 
    - default proxy point select   : QR
 * KERNEL 2: Exponential_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = exp(-l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 3: Exponential_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = exp(-l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 4: Matern32_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = (1+sqrt(3)*l) exp(-sqrt(3)*l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 5: Matern32_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = (1+sqrt(3)*l) exp(-sqrt(3)*l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 6: Matern52_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = (1+sqrt(5)*l+5/3*l) exp(-sqrt(5)*l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 7: Matern52_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : one parameter l: K(x,y) = (1+sqrt(5)*l+5/3*l) exp(-sqrt(5)*l|x-y|) 
    - default proxy point select   : QR
 * KERNEL 8: Quadratic_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 2
    - default kernel parameters    : 1.000000 -0.500000 
    - kernel parameter description : two parameters c,a: K(x,y) = (1+c*|x-y|^2)^a
    - default proxy point select   : QR
 * KERNEL 9: Quadratic_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 2
    - default kernel parameters    : 1.000000 -0.500000 
    - kernel parameter description : two parameters c,a: K(x,y) = (1+c*|x-y|^2)^a
    - default proxy point select   : QR
 * KERNEL 10: Coulomb_3D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 3
    - number of kernel parameters  : 0
    - default kernel parameters    : 
    - kernel parameter description : no parameters
    - default proxy point select   : proxy surface
 * KERNEL 11: Laplace_2D
    - dimension of kernel output   : 1 * 1
    - dimension of input point     : 2
    - number of kernel parameters  : 0
    - default kernel parameters    : 
    - kernel parameter description : no parameters
    - default proxy point select   : proxy surface
 * KERNEL 12: Stokes_3D
    - dimension of kernel output   : 3 * 3
    - dimension of input point     : 3
    - number of kernel parameters  : 2
    - default kernel parameters    : 1.000000 0.100000 
    - kernel parameter description : two parameters: (eta,  a)
    - default proxy point select   : proxy surface
 * KERNEL 13: RPY
    - dimension of kernel output   : 3 * 3
    - dimension of input point     : 3
    - number of kernel parameters  : 1
    - default kernel parameters    : 1.000000 
    - kernel parameter description : two parameters: eta
    - default proxy point select   : proxy surface

Examples

• ./pyh2pack/example.py