src.functions package

Submodules

src.functions.analytics module

author: Shahin (Amir Hossein) Rabbani
contact: shahin.rab@gmail.com
copyright: See License

src.functions.analytics.compare_jacobi_3methods_2d(opts)

Compare solutions to the 2D Poisson equation for 3 cases: Jacobi in the matrix form vs vector form vs Poisson Filters method (convolution, reduced or full kernel). No boundary treatment is done, which means data domain is treated as infinite domain.

See demos.convergence.demo_3methods_comparison_no_bc_2d() for the example demo.

Parameters

opts (OptionsGeneral) –

Returns

err_abs_jm_vs_jv : absolute error Jacobi matrix form vs Jacobi vector form
err_rel_jm_vs_jv : relative error Jacobi matrix form vs Jacobi vector form in percent
err_abs_jv_vs_pk : absolute error Jacobi vector form vs Poisson kernel
err_rel_jv_vs_pk : relative error Jacobi vector form vs Poisson kernel in percent
err_abs_jm_vs_pk : absolute error Jacobi matrix form vs Poisson kernel
err_rel_jm_vs_pk : relative error Jacobi matrix form vs Poisson kernel in percent
jacobi_solution_matrix_form : Jacobi solution in matrix form
jacobi_solution_vector_form : Jacobi solution in vector form
poisson_solution : Poisson kernel solution
data_padded : generate sample data

src.functions.analytics.compare_jacobi_3methods_neumann_2d(opts)

Compare solutions to the 2D Poisson equation for 3 cases: Jacobi in the matrix form vs vector form vs Poisson Filters method (convolution, reduced or full kernel).

This is with Neumann boundary treatment.

See demos.convergence.demo_3methods_comparison_no_bc_2d() for the example demo.

Parameters

opts (OptionsGeneral) –

Returns

err_abs_jm_vs_jv : absolute error Jacobi matrix form vs Jacobi vector form
err_rel_jm_vs_jv : relative error Jacobi matrix form vs Jacobi vector form in percent
err_abs_jv_vs_pk : absolute error Jacobi vector form vs Poisson kernel
err_rel_jv_vs_pk : relative error Jacobi vector form vs Poisson kernel in percent
err_abs_jm_vs_pk : absolute error Jacobi matrix form vs Poisson kernel
err_rel_jm_vs_pk : relative error Jacobi matrix form vs Poisson kernel in percent
jacobi_solution_matrix_form : Jacobi solution in matrix form
jacobi_matrix_form_residuals : Jacobi solution in matrix form residuals
jacobi_solution_vector_form : Jacobi solution in vector form
jacobi_vector_form_residuals : Jacobi solution in vector form residuals
poisson_solution : Poisson kernel solution
poisson_residual : Poisson kernel solution residuals
data_padded : generate sample data

src.functions.analytics.compare_jacobi_poisson_neumann_edge_mirrored_correction_2d(opts)

Compare ground truth Jacobi solution to Poisson filters solution (only for inverse) using tiled mirror method to treat the Neumann boundary. This is the core principle of the mirror marching algorithm proposed in the paper to do the boundary treatment.

The mirrored data tiles are automatically computed based on the target Jacobi iteration to make sure there is enough padding for the base data matrix when doing the convolution. See demos.boundary.demo_wall_neumann_tiled_mirrored_2d() for an example demo.

Parameters

opts (OptionsGeneral) –

Returns

jacobi_solution : ground truth solution
poisson_solution : Poisson filter solution
data_domain : created data matrix
data_mirrored : expanded data matrix used in the Poisson filters method
safe_rank : computed safe rank after kernel reduction
err_abs : absolute error
err_rel : relative error

src.functions.analytics.compute_adaptive_truncation_factors(opts, filters_iterations=None)

Given a list of target iterations, collect info about the effective maximum number of rank and filter size required for each iteration, for either 2D or 3D case.

In case of filters_iterations=None force generate the filters, for a range of target iterations [1, max_itr], where max_itr is the same as target iteration set in the options.

See demos.convergence.demo_adaptive_truncation_info_2d() for the example demo.

Warning

Filters are assumed to be already truncated upon load or generation.

Parameters

opts (OptionsGeneral) –
filters_iterations (List[int]) – list of target iterations. If filters are already available use this parameter, else (None) force compute them.

Returns

output will be a (max_itr, 3) shape, with format (max_effective_rank, max_filter_size, actual_filter_size). If in dynamic truncation mode return the analysis, else a zero vector with the size of max_itr

src.functions.analytics.generate_collision_mask(M, path, load_shape_boundary)

Create a 2D collision masks with the same size as the input data matrix. Walls are marked as both contours and solids. Do cell marking for solid/contour/data for complex collider shapes, if loaded.

It is optional to load a complex shape, in which case we do proper sampling to best match the resolutions of the loaded image and the domain matrix.

The mask marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

M (ndarray) – 2D data domain
path (str) – path to load the image of the complex object
load_shape_boundary (bool) – load a complex object shape from file

Returns

a collision mask with +1 flags for solid, and a contour mask for the contours around the solid object with +2 flags.

src.functions.analytics.make_special_data_central_square(size, radius, value)

Make a 2D data matrix with a non-zero box at the center and zero everywhere else.

Parameters

size (tuple) – base size of the square matrix
radius (float) – not a circle radius, but the span length of the non-zeros in the central box. Take it as base square length.
value (float) – non-zero value to be assigned to the pattern

Returns

2d data matrix with a non-zero box pattern at the center

src.functions.analytics.prepare_padded_input_matrix(opts)

Cook up a data matrix (2D or 3D), either with randomized or constant values, with proper dynamic wall padding.

We expand the matrix to account for convolution shrinkage. This helps with comparing the Jacobi and the Poisson kernel solutions as the two might end up having different sizes.

Note

Check out helper.commons.DataMatrixMode() for choices of making the data matrix.

Parameters

opts (OptionsGeneral) –

Returns

data_padded: padded (expanded) data matrix
data_domain: original un-padded data matrix
padding_size: computed padding size due to convolution shrinkage

src.functions.analytics.solve_poisson_2d(M, poisson_kernel, opts, skip_margin=0)

General solver for 2D Poisson equation, for either inverse or forward Poisson.

The function supports solving the Poisson equation with the full kernel (input), reduced kernel (computed here), or separable Poisson filters (computed here), depending on how parameters are set in the input options. Truncation parameters will be retrieved from options as well.

Note

To set rank, order, and other relevant parameters you need to pack OptionsKernel and OptionsReduction in OptionsGeneral and send it to this function. These dataclasses are in helper.commons.py. To see how to pack options, look at main demos, or see helper.commons.generic_options().

Parameters

M (ndarray) – input 2D domain matrix
poisson_kernel (ndarray) – precomputed 2D Poisson kernel (inverse or forward)
opts (OptionsGeneral) –
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Default= 0)

Returns

solution to the 2D Poisson equation, computed safe rank

src.functions.analytics.subtract_contour_from_solid(solid, contour)

Given a 2D solid object mask and its corresponding contour, extract the interior part of the solid maks by subtracting the two masks.

Parameters

solid (ndarray) – 2d input matrix
contour (ndarray) – 2d input matrix

Returns

new 2d matrix with marked interior cells by 1, else 0

src.functions.decompositions module

author: Shahin (Amir Hossein) Rabbani
contact: shahin.rab@gmail.com
copyright: See License

Technicals

While in 2D the eigenvalues are sorted, it is not guaranteed to be true in 3D. This is due to the lack of consensus on rank definition in tensor decomposition, and the fact that current symmetric CP algorithm used by our method, is an iterative method that does its best to deal with a NP-hard problem.

This difference between 2D and 3D eigenvalues properties has an important implication in practice:

In 2D it is safe to generate, say, rank-8 filters and only use, for example, the first 3 to capture a significant portion of the variance. The contribution of the 4th rank is expected to be less than the first 3 ranks.

In 3D we might get different filters if we first get, say, rank-8 filters and choose the first 3, compared to when we directly compute rank-3 filters and use all of them. In both cases we are interested in the first 3 ranks, but in the latter case results are expected to be numerically more stable. There is no guarantee to get a decreasing variance contribution with higher ranks in 3D, just like what one would expect in 2D.

Poisson filter computation from decomposition components

Check out get_filters_from_svd_components_2d() to see how we do filter computation. The same method applied to 3D using Sym-CP and with an additional set of vectors for the 3rd dimension.

Absorbing eigenvalues in 2D and 3D

2D : Because the Poisson kernel is always square with non-negative eigenvalues, the singular values obtained from SVD coincide with the eigenvalues. This is not true in general, and is only true in our case given the aforementioned properties of the kernel. To absorb the eigenvalues in filters, we scale each of the two ranked filters by \(\sqrt{\sigma_{r}}\), the eigenvalue corresponding to rank \(r\).
3D : Symmetric-CP decomposition of 3D Poisson kernel gives a set of factors and cores, where cores can be taken as the best approximation of true tensor eigenvalues (this should always be taken with a grain of salt for tensor eigen decomposition). To absorb the eigenvalues in filters, we scale each of the three ranked filters by \(\sqrt[3]{\sigma_{r}}\), the eigenvalue corresponding to rank \(r\).

Note

Safe Rank

The Poisson kernel is always half-rank, meaning its rank is equal to half of its size along any of its dimensions.
We extensively use a safe rank to protect the user from asking for invalid desired ranks in 2D. A safe rank is the minimum of the Poisson kernel actual maximum possible rank and the desired rank input by the user. The main function to do this is rank_safety_clamp_2d() in 2D. There is no need for a rank clamp in 3D because of how CP decomposition works.
In our functions ‘safe’ simply means the number of ranked filters we would like to include in our convolutions, and use it to help with setting up loops etc.
While in 2D safe means a rank that does not exceed the actual rank of the kernel, in 3D it is different. Due to the lack of a clear definition of rank in tensor decomposition, we can pretty much ask for any rank in the CP-vew eigen decomposition and always get something that “works”. This vagueness of a proper rank definition in 3D is fundamental, which partially contributes to the fact that tensor decomposition is NP-hard.

src.functions.decompositions.absorb_cores_3d(cores, factors)

Absorb the core weights (eigenvalues) into filters (eigenvectors). Scale each filter by the cubic root of the corresponding core weight to get the Poisson filters.

Warning

The absorption of eigenvalues only works when doing the tensor decomposition using Symmetric CP.
Without this step we do not really have Poisson filters, only eigenvectors.
This is based on (filter, rank) order, which is the opposite of 2D filter order.

Note

Symmetric-CP decomposition of 3D Poisson kernel gives a set of factors and cores, where cores can be taken as the best approximation of true tensor eigenvalues (this should always be taken with a grain of salt for tensor eigen decomposition as there is no consensus on rank definition in 3D and the tensor eigen decomposition is NP-hard). To absorb the eigenvalues in filters, we scale each of the three ranked filters by \(\sqrt[3]{\sigma_{r}}\), the eigenvalue corresponding to rank \(r\).

Parameters

cores (ndarray) – eigenvalues
factors (ndarray) – eigenv ectors

Returns

3D Poisson filters

src.functions.decompositions.compute_low_rank_kernel_from_filters_2d(hor, ver, safe_rank)

Compute a low rank kernel from the given Poisson filters. Filters can be in the original form or can be truncated.

Parameters

hor (ndarray) – 2d array: horizontal filters as 2-tuple (rank, filter)
ver (ndarray) – 2d array: vertical filters as 2-tuple (rank, filter)
safe_rank (int) – cumulative rank, i.e. maximum rank to be included. must be safe.

Returns

square low rank kernel

src.functions.decompositions.compute_nth_kernel_mode_2d(hor, ver, rank)

Reconstruct the rank-1 matrix corresponding to the \(r\) -th rank (mode) from horizontal and vertical Poisson filters for a desired rank.

Warning

Singular values must be already absorbed into the filters.

Parameters

hor (ndarray) – 2d array: horizontal filters as 2-tuple (rank, filter)
ver (ndarray) – 2d array: vertical filters as 2-tuple (rank, filter)
rank (int) – desired safe rank. rank needs to be equal or less than the original decomposed matrix.

Returns

\(r\) -th mode square matrix constructed from the \(r\) -th pair of filters

src.functions.decompositions.compute_nth_kernel_mode_3d(ranked_filter)

Reconstruct the rank-1 tensor corresponding to the \(r\) -th rank (mode) from the Poisson filter of a certain rank.

Note

All 3 dimensions use the same filter.

Warning

Core values (eigenvalues) must be already absorbed into the filters.

Parameters: ranked_filter (ndarray) – single Poisson filter, 1d array
Returns: \(r\) -th mode tensor constructed from the \(r\) -th ranked filter

src.functions.decompositions.compute_separable_filters_trunc_adaptive_2d(U, S, VT, rank, trunc_threshold, preserve_shape=True)

Compute the Poisson filters from SVD components and truncate them based on a fixed cut-off threshold.

Note

It does not guarantee the returned filters have the same size due to adaptive truncation unless preserve_shape=True (Default).

Parameters

U (ndarray) – \(\mathbf{U}\) in \(\mathbf{USV^T}\)
S (ndarray) – \(\mathbf{S}\) in \(\mathbf{USV^T}\)
VT (ndarray) – \(\mathbf{V^T}\) in \(\mathbf{USV^T}\)
rank (int) – desired rank. It will be safely clamp if larger than the input matrix actual rank.
trunc_threshold – truncation threshold (absolute value)
preserve_shape – if True keep the original shape and fill them with zeros (Default= True)

Returns

stacked array of truncated horizontal and vertical separable filters for \(n\) ranks, in the form of a 2-tuple (rank, 1d filter). The clamped safe rank is also returned.

src.functions.decompositions.compute_separable_filters_trunc_percent_2d(U, S, VT, rank, trunc_factor)

Compute the Poisson filters from SVD components and truncate a certain percentage of them.

Since filters are symmetrical, the truncation is applied to both sides of the filters.

Parameters

U (ndarray) – \(\mathbf{U}\) in \(\mathbf{USV^T}\)
S (ndarray) – \(\mathbf{S}\) in \(\mathbf{USV^T}\)
VT (ndarray) – \(\mathbf{V^T}\) in \(\mathbf{USV^T}\)
rank (int) – desired rank. It will be safely clamp if larger than the input matrix actual rank.
trunc_factor (float) – truncation percentage in [0, 1]

Returns

stacked array of truncated horizontal and vertical separable filters for \(n\) ranks, in the form of a 2-tuple (rank, 1d filter). The clamped safe rank is also returned.

src.functions.decompositions.compute_separable_filters_truncated_2d(U, S, VT, rank, trunc_method, trunc_factor, preserve_shape)

Compute and truncate 2D Poisson filters.

Truncation is either based on cutting off a certain percentage of the filter, or using a fixed cut-off threshold (adaptive truncation).

Parameters

U (ndarray) – \(\mathbf{U}\) in \(\mathbf{USV^T}\)
S (ndarray) – \(\mathbf{S}\) in \(\mathbf{USV^T}\)
VT (ndarray) – \(\mathbf{V^T}\) in \(\mathbf{USV^T}\)
rank (int) – desired rank. Does not have to be a safe rank
trunc_method (TruncationMode) – either PERCENTAGE or FIXED_THRESHOLD
trunc_factor (float) – percentage in [0, 1] if PERCENTAGE, fixed cut-off threshold if FIXED_THRESHOLD
preserve_shape (bool) – if True keep the original shape and fill them with zeros, else return the shrunk filter

Returns

stacked array of truncated horizontal and vertical separable filters for \(n\) ranks, in the form of a 2-tuple (rank, 1d filter). The clamped safe rank is also returned.

src.functions.decompositions.demo_test_svd_absorb()

Sanity check: Testing the eigenvalue absorption into the separable filters in 2D.

We expect the low-rank matrix obtained from the reduced space exactly matching the reconstruction from the outer product of horizontal and vertical filters. Here by “matching exactly” we mean machine precision.

Because the Poisson kernel is always square with non-negative eigenvalues, the singular values obtained from SVD coincide with the eigenvalues. This is not true in general, and is only true in our case given the special properties of the Poisson kernel, like symmetry.

To absorb the eigenvalues in filters, we scale each of the two ranked filters by \(\sqrt{\sigma_{r}}\), the eigenvalue corresponding to rank \(r\).

See get_filters_from_svd_components_2d() for more details.

src.functions.decompositions.get_filters_from_svd_components_2d(rank, U, S, VT)

Compute separable Poisson filters using the Canonical Polyadic (CP) view of the kernel matrix.

We do not need to explicitly perform Canonical Polyadic Decomposition (CPD) in 2D. Instead, we can decompose the matrix into a summation of rank-1 matrices (as CPD does) by manipulating the components obtained from Singular Value Decomposition (SVD) of the same matrix. This involves reconstructing each rank-1 matrix by multiplying columns of \(U\) and rows of \(V^T\) scaled by the squared root of the corresponding singular value.

Given the SVD of a square kernel \(\mathbf{F}\)

\[\mathbf{F} \approx \mathbf{USV^T}\]

where:

The columns of \(\mathbf{U}\) are the left-singular vectors of \(\mathbf{F}\)
The columns of \(\mathbf{V}\) are the right-singular vectors of \(\mathbf{F}\)
The values along the diagonal of \(\mathbf{S}\) are the singular values of \(\mathbf{F}\)

The CP view is obtained by

\[\mathbf{F} \approx \displaystyle\sum_{j=1}^n s_j \mathbf{u}_j \otimes \mathbf{v}_j\]

where \(\mathbf{u}_j\) and \(\mathbf{v}_j\) are the \(j\) -th columns of the left- and right-singular vectors, \(s_j\) is the \(j\) -th singular value, and \(\otimes\) is the outer product. The summation gives \(n\) rank-1 matrices, each of which meets the rank-1 separability condition to compute our filters, where \(n\) can be the maximum rank of \(\mathbf{F}\) or less.

With the above formulation the \(j\) -th component corresponds to the \(r\) -th rank in the CP view. This means the \(\mathbf{u}_j\) and \(\mathbf{v}_j\) can be taken as the separable filters needed to reconstruct the \(r\) -th mode \(\mathbf{F}_r\), scaled by its singular value \(s_r\) (which is the same as its eigenvalue, as explained in the next section)

\[\begin{split}\mathbf{F_r} = s_r (\mathbf{u_r} \otimes \mathbf{v_r})=s_r \left[\begin{array}{c} u_{r_1} \\ u_{r_2} \\ \vdots \\ u_{r_n} \end{array}\right] \otimes \left[\begin{array}{c} v_{r_1} \\ v_{r_2} \\ \vdots \\ v_{r_n} \end{array}\right] = s_r \left[\begin{array}{cccc} u_{r_1} v_{r_1} & u_{r_1} v_{r_2} & \cdots & u_{r_1} v_{r_n} \\ u_{r_2} v_{r_1} & u_{r_2} v_{r_2} & \cdots & u_{r_2} v_{r_n} \\ \vdots & \vdots & \ddots & \vdots \\ u_{r_n} v_{r_1} & u_{r_n} v_{r_2} & \cdots & u_{r_n} v_{r_n} \end{array}\right]\end{split}\]

where \(\mathbf{F_r}\) is the \(r\) -th rank-1 square matrix component of \(\mathbf{F}\) in the CP view. As the final step to get the Poisson filters we need to absorb the singular values into the filters.

Absorbing singular values into separable filters

Because the Poisson kernel is always square with non-negative eigenvalues, the singular values obtained from SVD coincide with the eigenvalues. This is not true in general, and is only true in our case given the aforementioned properties of the kernel. To absorb the eigenvalues in filters, we scale each of the two ranked filters by \(\sqrt{\sigma_{r}}\), where \(\sigma_r = s_j = s_r\), the \(r\) -th singular value.

The \(r\) -th vertical \(f_v\) and horizontal \(f_h\) Poisson filters reconstructing the \(r\) -th mode are

\[\begin{split}\mathbf{f_{v_r}}=\sqrt{\sigma_{r}} \left[\begin{array}{c} u_{r_1} \\ u_{r_2} \\ \vdots \\ u_{r_n} \end{array}\right], \; \; \; \mathbf{f_{h_r}}=\sqrt{\sigma_{r}} \left[\begin{array}{c} v_{r_1} \\ v_{r_2} \\ \vdots \\ v_{r_n} \end{array}\right]\end{split}\]

which are then used to compute \(\mathbf{F_r} = \mathbf{f_{v_r}} \otimes \mathbf{f_{h_r}}\)

Note

In theory, the horizontal and vertical Poisson filters are the same due to the symmetry in the full Poisson kernel. However, and in practice, SVD might give singular vectors with opposing signs. To avoid any problem, we save and apply each of the filters separately.

This is not a problem in 3D because of the way Sym-CP does the decomposition.

Parameters

rank (int) – desired rank. Does not have to be a safe rank
U (ndarray) – \(\mathbf{U}\) in \(\mathbf{USV^T}\)
S (ndarray) – \(\mathbf{S}\) in \(\mathbf{USV^T}\)
VT (ndarray) – \(\mathbf{V^T}\) in \(\mathbf{USV^T}\)

Returns

horizontal and vertical Poisson filters

src.functions.decompositions.get_max_rank_2d(itr, zero_init)

Get the maximum possible rank of a Poisson kernel with the given target Jacobi iteration.

Parameters

itr (int) – target Jacobi iteration
zero_init (bool) – if we are zero-initializing the Jacobi solution, in which case the corresponding Poisson kernel will have a smaller size. Typical value is True because we usually start from \(x=0\), i.e. no warm starting.

Returns

maximum possible rank

src.functions.decompositions.poisson_components_2d(opts)

Generate 2d Poisson kernel and decompose it using SVD to get \(\mathbf{USV^T}\)

Note

To set rank, order, and other relevant parameters you need to pack OptionsKernel and OptionsReduction in OptionsGeneral and send it to this function. These dataclasses are in helper.commons.py. To see how to pack options, look at main demos, or see helper.commons.generic_options().

Parameters: opts (OptionsGeneral) – parameters bundle
Returns: \(\mathbf{U}\), \(\mathbf{S}\), and \(\mathbf{V^T}\) components, as well as the low rank kernel with the same size as the input with the desired rank

src.functions.decompositions.poisson_compute_modes_trim_filters_3d(opts, rank_filter_reorder, filter_trim_zeros=False, preserve_shape=True)

Compute 3D ranked modes.

Parameters

opts (OptionsGeneral) – parameters bundle
rank_filter_reorder (bool) – if True change 2-tuple order from default (filters, ranks) to (ranks, filters), which is more consistent with 2D filters. Typical value is True.
filter_trim_zeros (bool) – if True do not preserve shape, cut the zeros on each side of the filter (Default= False)
preserve_shape (bool) – if True keep the original shape and fill them with zeros, else return the shrunk filter (Default= True)

Returns

3D ranked modes, Poisson filters, a low rank reconstructed kernel along with the full kernel

src.functions.decompositions.poisson_decomposition_components_3d(opts)

Compute 3D cores, factors, low rank kernel and full kernel.

Note

We use Symmetric CP as the main decomposition method. Tucker decomposition is also available for experimentation.

Parameters: opts (OptionsGeneral) – parameters bundle
Returns: reduced low rank kernel, cores and factors after decomposition, as well as the full (i.e. unreduced) kernel

src.functions.decompositions.poisson_filters_2d(opts)

Generate 2D Poisson filters.

First generate the kernel then compute the filters using singular value decomposition (SVD). Check out get_filters_from_svd_components_2d() to see how we do filter computation.

Note

To set rank, order, and other relevant parameters you need to pack OptionsKernel and OptionsReduction in OptionsGeneral and send it to this function. These dataclasses are in helper.commons.py. To see how to pack options, look at main demos, or see helper.commons.generic_options().

Parameters: opts (OptionsGeneral) – parameters bundle
Returns: horizontal and vertical filters

src.functions.decompositions.poisson_filters_3d(opts, rank_filter_reorder, preserve_shape=True)

Generate 3D Poisson filters (with truncation if desired).

First generate the kernel then compute the filters using Eigenvalue Decomposition of a 3D kernel tensor. See poisson_decomposition_components_3d() for decomposition methods.

Note

Check out get_filters_from_svd_components_2d() to learn about the basics of Poisson filter computation in 2D. For 3D filters apply the same principles except for an extra filter in 3D dimension. The outer product \(\otimes\) in 2D becomes tensor product in 3D.

Note

To set rank, order, and other relevant parameters you need to pack OptionsKernel and OptionsReduction in OptionsGeneral and send it to this function. These dataclasses are in helper.commons.py. To see how to pack options, look at main demos, or see helper.commons.generic_options().

Parameters

opts (OptionsGeneral) – parameters bundle
rank_filter_reorder (bool) – if True change 2-tuple order from default (filters, ranks) to (ranks, filters), which is more consistent with 2D filters. Typical value is True.
preserve_shape – if True keep the original shape and fill them with zeros, else return the shrunk filter (Default= True)

Returns

3 Poisson filters for horizontal, vertical and depth (fiber) passes. The filters can be already truncated if truncation parameters are set in opts
generated Poisson kernel tensor
reduced (low rank) reconstructed Poisson kernel corresponding to the desired rank specified in opts
safe rank

src.functions.decompositions.poisson_svd_2d(P, rank=1)

Singular value decomposition (SVD) of the 2D Poisson kernel.

Decompose, reduce and reconstruct a low rank kernel by truncating the singular values based on the given desired rank.

Note

Useful links

Parameters

P (ndarray) – input Poisson kernel matrix
rank (int) – desired rank. Must be a safe rank, i.e. the minimum of the Poisson kernel actual maximum possible rank and the desired rank input by the user. See notes and rank_safety_clamp_2d()

Returns

\(\mathbf{U}\), \(\mathbf{S}\), and \(\mathbf{V^T}\) components in svd(P) = \(\mathbf{USV^T}\), as well as the low rank kernel with the same size as the input and reduced to have maximum a rank as the desired input rank.

src.functions.decompositions.poisson_symcp_3d(kernel, rank)

Decompose a full 3D Poisson kernel using symmetric CP (Canonical Polyadic) decomposition.

Warning

For a kernel of shape (1, 1, 1) corresponding to a target iteration of 1, Sym-CP does not work. To avoid bugs or crashes, we return a set of zero \(1 \times 1\) vectors (as many as target rank), with the first factor being the only non-zero vector and with the value equal to the only scalar value we have from the (1, 1, 1) kernel. Cores will be all ones. The reduced kernel will be the same as the original full kernel (there is no reduction).

Parameters

rank (int) – target decomposition rank
kernel (ndarray) – full 3D Poisson kernel

Returns

reduced (lowe rank) kernel and decomposition cores and factors

src.functions.decompositions.poisson_tucker_3d(kernel, rank)

Decompose a full 3D Poisson functions using tensor Tucker decomposition.

Parameters

rank (int) – target decomposition rank
kernel (ndarray) – full 3D Poisson kernel

Returns

reduced (lowe rank) kernel and decomposition cores and factors

src.functions.decompositions.rank_safety_clamp_2d(itr, rank, zero_init)

Clamp the given desired rank based on the maximum possible rank of a Poisson kernel with the given target Jacobi iteration.

The Poisson kernel is always half-rank, meaning its rank is equal to half of its size along any of its dimensions. A safe rank is the minimum of the Poisson kernel actual maximum possible rank and the desired rank.

Parameters

itr (int) – target Jacobi iteration
rank (int) – desired rank
zero_init (bool) – if we are zero-initializing the Jacobi solution, in which case the corresponding Poisson kernel will have a smaller size. Typical value is True because we usually start from \(x=0\), i.e. no warm starting.

Returns

safe rank

src.functions.generator module

author: Shahin (Amir Hossein) Rabbani
contact: shahin.rab@gmail.com
copyright: See License

Technicals

We use a recursive function to generate the Poisson kernel. Due to similarities between sub-problems, we enhance the recursion performance by memoization, i.e. saving the solution and reusing it multiple times. See poisson_kernel_2d() and poisson_kernel_3d().

Notes

Choosing the kernel type UNIFIED vs STANDARD:

UNIFIED
The UNIFIED type produces the exact same parametric kernel for both INVERSE and FORWARD Poisson equations. While for the FORWARD version the kernel already supports warm starting, for the INVERSE version it only works with a zero initial guess. Also, we will need to negate the input to make it work for INVERSE applications.

Upside: Same kernel generator for both INVERSE and FORWARD. No pre-assumption about the type of the Poisson equation we are solving.

Downside: for INVERSE Poisson (e.g. poisson pressure) you would need to negate the input divergence before feeding it to the solver.

STANDARD (Safer)
The STANDARD type has an alternating sign for \(\alpha\) in the parametric kernel, so the two kernels for FORWARD and INVERSE are not exactly the same, but it supports warm starting for both FORWARD and particularly INVERSE. An example of warm starting is when if you want to initialize the first guess of the pressure in Poisson-pressure solve from RHS (divergence).

Upside: Supports warm starting (but within the constraints of Poisson filters applications - see paper). Less confusion in applying filters and convergence analysis.

Downside: Needs a separate call to each FORWARD and INVERSE Poisson kernel generators, which is not a major inconvenience.

STANDARD is the recommended default type.

Application Notes:

In addition to the target Jacobi iteration, \(\alpha\) and \(\beta\) are the two major players when generating the Poisson kernel. Things like the unit cell size \(dx\), time step \(dt\) and diffusivity factor \(\kappa\) (where \(dt\) and \(\kappa\) are only used by the FORWARD Poisson kernel), are captured by \(\alpha\) and \(\beta\). Check out compute_alpha() and compute_beta() to see how they are set. In practice, you would only need to provide \(dx\), \(dt\) and \(\kappa\), and the generator automatically computes \(\alpha\) and \(\beta\) based on the problem dimension and the Poisson equation type.
\(\alpha\) has a negative sign in the original Jacobi equation for pressure. To make the kernel generation consistent for both FORWARD and INVERSE we treat \(\alpha\) as positive for both cases. So for example to solve Poisson-pressure (inverse Poisson) you should multiply the rhs divergence by \(-1\) to account for the sign change.
Warm starting: often found as zero_init parameter in some functions, Poisson filters have limited support for warm starting the linear solve (as discussed in the paper). See i_want_to_warm_start() and is_ok_to_zero_init().

Tip

If you decide to change/improve things:

You can generate Poisson kernels with \(\alpha=1\) then just scale \(B\) in \(L*X = B\) (the input matrix) with the actual \(\alpha\) you are interested in. This means if doing only INVERSE you can factorize \(\alpha\) from the kernel generation and change it in real time in your application (multiplying filter convolution results by \(\alpha\)). However, in the FORWARD case, and due to the presence of \(\alpha\) in the way \(\beta\) is computed, we can’t do this factorization.
In the current implementation, we decided to explicitly include \(\alpha\) in the kernel to have a uniform formulation that works for both FORWARD and INVERSE.

src.functions.generator.compute_alpha(dx, dt, kappa, solver_type)

Compute \(\alpha\) that is used in kernel generation, for both 2D and 3D (see poisson_kernel_2d() and poisson_kernel_3d() for application).

In forward and inverse Poisson equations we have different equations:

Forward: \(\alpha = +{(dx)}^2/{(\kappa . dt)}\)
Inverse: \(\alpha = -{(dx)}^2\)

where:

\(dx:\) cell size
\(dt:\) time step (if forward Poisson is intended)
\(\kappa\): diffusivity (if forward Poisson is intended)

Parameters

dx (float) – cell size
dt (float) – time step
kappa (float) – diffusivity
solver_type (PoissonSolverType) – inverse or forward

Returns

\(\alpha\)

src.functions.generator.compute_beta(alpha, solver_type, dim)

Compute \(\beta\) that is used in kernel generation, for both 2D and 3D (see poisson_kernel_2d() and poisson_kernel_3d() for application).

\(\beta\) values based on dimension and Poisson solver type:

2D Inverse: \(\beta = 4\)
2D Forward: \(\beta = 4 + \alpha\)
3D Inverse: \(\beta = 6\)
3D Forward: \(\beta = 6 + \alpha\)

Parameters

alpha (float) – see compute_alpha()
solver_type (PoissonSolverType) – inverse or forward
dim (SolverDimension) – solver dimension, 2D or 3D

Returns

\(\beta\)

src.functions.generator.get_alpha_sign(solver_type, kernel_type)

Compute \(\alpha\) sign that is used in kernel generation, for both 2D and 3D (see poisson_kernel_2d() and poisson_kernel_3d() for application).

In forward and inverse Poisson equations we have different equations, and hence different \(\alpha\) signs:

Forward: \(\alpha = +{(dx)}^2/{(\kappa . dt)}\)
Inverse: \(\alpha = -{(dx)}^2\)

where:

\(dx:\) cell size
\(dt:\) time step (if forward Poisson is intended)
\(\kappa\): diffusivity (if forward Poisson is intended)

Parameters

solver_type (PoissonSolverType) – inverse or forward
kernel_type (PoissonKernelType) – unified or standard

Returns

sign of alpha

src.functions.generator.get_kernel_size(itr, zero_init)

Get the expected Poisson kernel size along each dimension for a given target Jacobi iteration (order) and based on whether we want to warm start or not.

Note

The kernel is always either square (2D matrix) or cube (3D tensor). The kernel actual size for a given target Jacobi iteration \(i\) is:

zero_init=True : \({(2*i-1)}^d\)
zero_init=False : \({(2*i+1)}^d\)

for dimension \(d=2,3\). This function only returns the base size and not the kernel actual size.

Parameters

itr (int) – target Jacobi iteration
zero_init (bool) – see is_ok_to_zero_init()

Returns

kernel size along each dimension

src.functions.generator.get_kernel_size_from_half_size(half_size)

Given a half Poisson filter size compute the actual kernel size.

Parameters: half_size (int) – half filter size excluding the element in the middle.
Returns: kernel size

src.functions.generator.i_want_to_warm_start(warm_start, solver_type)

Automatically generate parameters for the kernel and solver if we want to warm start the solution to \(Ax=b\).

Note

If you are doing inverse Poisson and

want to avoid negating your input matrix \(B\) (in the matrix-form \(LX=B\)), or
want to warm start the solution with the input matrix

then it is better to use warm start.

src.functions.mathlib module

author: Shahin (Amir Hossein) Rabbani
contact: shahin.rab@gmail.com
copyright: See License

Technicals

The Poisson solver has a linear formulation like in the typical \(Ax=b\) but with a difference of replacing vectors \(x\) and \(b\) with matrices, and replacing matrix-vector multiplication with a convolution operator. We use a proper notation with a difference naming convention:

\(Ax = b \leftrightarrow L*X = B\)

where:

\(L\): Laplacian operator
\(*\): convolution operator
\(X\): given or unknown matrix, depending on the Poisson equation we are solving for.
\(B\): given or unknown rhs matrix, depending on the Poisson equation we are solving for.

LHS and RHS are both in matrix-forms, i.e. all matrices/tensors have the same dimension.

Note

We call \(Ax=b\) the vector-form and \(L*X=B\) the matrix-form.

Note

The Poisson kernel (and subsequently its separable filters) are computed based on the matrix-form, so for an infinite domain without any boundary treatment the solution using the Poisson kernel and the one from the matrix-form must match to the machine precision.

The linear solver works with two versions of the Poisson equation, namely inverse and forward :

Inverse : An example is Poisson-pressure where \(B\) is the input divergence and \(X\) is the output pressure.

Forward : An example is diffusion where \(X\) is the input density and the output \(B\) is the diffused quantity.

Depending on what we are solving for, we have different setups for input/output:

Inverse Poisson equation : Given \(M\) as input in \(L*X = M\), obtain solution \(X\) by implicitly approximating \(L^{-1}\) in \(X = L^{-1}*M\).

Forward Poisson equation : Given \(M\) as input, obtain the solution \(B\) to the diffusion equation, i.e. perform \(L*M = B\) , where \(B\) is the output.

The solution is computed using an ‘implicit’ finite difference method, Just like Jacobi.

Note

Note how the input matrix \(M\) changes role in \(L*X=B\) based on the type of the Poisson equation setup.

We use multi-rank Poisson filters for a given rank. Given the Poisson equation in the matrix-form \(L*X=B\), the Poisson kernel \(L\) (in forward setup) and its inverse \(L^{-1}\) (in inverse setup) are already baked into the Poisson filters. Just provide the input data matrix and the corresponding filters matching the formulation setup you are interested in.

In general, we can replace \(L\) and \(L^{-1}\) with a unified kernel \(F\) that operates on a data domain matrix (or tensor), and perform Eigen decomposition on \(F\) to get the Poisson filters.

Convolution Order of Poisson Filters

Poisson filters are used in cascaded convolutions to get the solution.

In 3D, we have

\(F * M \approx \displaystyle\sum_{r=1}^n f_{v_r} * (f_{h_r} * (f_{d_r} * M))\)

where

\(F\) - Full Poisson kernel (either \(L\) or \(L^{-1}\))
\(M\) - Input data field
\(f_v\) - Vertical filter
\(f_h\) - Horizontal filter
\(f_d\) - Depth (fiber) filter
double subscript \(_r\) means the filter corresponding the current rank
\(\displaystyle\sum_{r=1}^n\) is multi-rank summation (i.e. modal solutions)

The convolution order goes from the inner brackets to outer brackets, meaning first we need to convolve \(M\) with the fiber filter, then convolve the results with the horizontal and vertical filters.

For multi-rank convolution we have separate and independent convolutions passes on \(M\), then sum up the results. The summation comes from the Canonical Polyadic Decomposition (CPD) view in our matrix/tensor decomposition setup (different in 2D and 3D), which makes it possible to have rank-1 kernel convolutions to get modal solutions taking care of different frequencies in the data domain.

Warning

DO NOT feed input data matrix \(M\) in outer bracket convolution.

ALWAYS use the results of the previous convolution pass to do the next one.

Also see solve_poisson_separable_filters_2d() and solve_poisson_separable_filters_wall_aware_3d() for the convolution order in 2D and 3D.

Note

The ground truth Jacobi solver we use in our code uses the same matrix-form setup as the Poisson Filters/Kernel setup. It is a straightforward practice to establish the connection between Jacobi matrix-form and the more widely known and used version, Jacobi vector-form .

src.functions.mathlib.apply_adaptive_truncation_1d(array_1d, safe_rank, cut_off, preserve_shape)

Adaptive truncate all ranked filters using a fixed threshold as the cut-off value, assuming symmetrical filters and values are sorted sideways from largest (center) to smallest (tales).

Parameters

array_1d (ndarray) – input filter
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel
cut_off (float) – truncation threshold (absolute value)
preserve_shape (bool) – if True keep the original shape and fill them with zeros

Returns

truncated (smaller) array. If preserving the shape, keep the shape and insert zeros in the truncated parts.

src.functions.mathlib.apply_laplacian_2d(X)

Convolve a 2d domain with a 2d Laplacian operator in the matrix from. This is based on a 3x3 base Laplacian kernel when finite differencing (the same as the Laplacian base kernel in Jacobi).

Note

matrix from of \(Ax=b\) is \(L*X=B\)

where

\(L\): Laplacian operator
\(*\): convolution operator
\(X\): given or unknown matrix, depending on the Poisson equation we are solving for.
\(B\): given or unknown rhs matrix, depending on the Poisson equation we are solving for.

All matrices have the same dimension in the matrix form.

Warning

Each dimension must be at least 3 to allow for fetching the marginal wall cells.

Parameters: X (ndarray) – input 2d matrix
Returns: \(L*X\)

src.functions.mathlib.apply_laplacian_3d(X)

Convolve a 3d domain with a 3d Laplacian operator in the matrix from. This is based on a 3x3x3 base Laplacian kernel when finite differencing (the same as the Laplacian base kernel in Jacobi).

Note

matrix from of \(Ax=b\) is \(L*X=B\)

where

\(L\): Laplacian operator
\(*\): convolution operator
\(X\): given or unknown tensor, depending on the Poisson equation we are solving for.
\(B\): given or unknown rhs tensor, depending on the Poisson equation we are solving for.

All tensors have the same dimension in the matrix form.

Warning

Each dimension must be at least 3 to allow for fetching the marginal wall cells.

Parameters: X (ndarray) – input 3d tensor
Returns: \(L*X\)

src.functions.mathlib.bivariate_gaussian_normalized(x, y, x0, y0, ef_rad)

Generate a normalized 2d gaussian.

Parameters

x (ndarray) – x values
y (ndarray) – y values
x0 (float) – mean x
y0 (float) – mean y
ef_rad (float) – effective radius, full-width-half-maximum

Returns

normalized 2d Gaussian

src.functions.mathlib.compute_abs_rel_error(M1, M2)

Compute absolute and relative errors, with M1 being the reference matrix

Parameters

M1 (ndarray) – reference matrix
M2 (ndarray) – test matrix

Returns

absolute error, relative error in %

src.functions.mathlib.compute_conv_padding(convolution_kernel_size)

Compute the padding required to compensate for shrinkage due to convolution.

Parameters: convolution_kernel_size (int) – size of a square convolution kernel
Returns: required padding size
Return type: int

src.functions.mathlib.compute_gradient_2d(M, grad_scale=2.0, half_dx=0.5)

Compute the gradient of a scalar field. The output would be a 2d vector field.

Warning

This computes the gradient for the central cell. If, for example, we are interested in a MACGrid setup we need to be careful how to interpret this.

Note

The default values are recommended: grad_scale = 2.0, half_dx = 0.5

Parameters

M (ndarray) – input scalar field
grad_scale (float) – gradient scale (Default= 2.0)
half_dx (float) – half the cell size (Default= 0.5)

Returns

2d vector field of gradients

src.functions.mathlib.compute_l1_norm_error(M1, M2)

Compute L1 norm of the difference between two matrices.

Parameters

M1 (ndarray) – input matrix
M2 (ndarray) – input matrix

Returns

error scalar

src.functions.mathlib.compute_norm(M, method)

Different norms of the matrix.

Parameters

M (ndarray) – input matrix
method (NormOptions) – frobenius, mse, infinite

Returns

error scalar

src.functions.mathlib.compute_residual_poisson_operator(X, B, solver_dimension)

Compute L2-norm residual \(r= L*X - B\) using 2d or 3d Laplacian \(L\) in the matrix form (equivalent to \(r=Ax-b\) in the vector form).

Warning

This is based on the matrix form setup, excluding wall boundaries by one cell, and is only valid for inverse Poisson equation. The forward Poisson version needs to be implemented.

Warning

This function is not safe if not excluding wall boundaries. Use this function in subdomains to properly exclude the wall boundaries.

Parameters

X (ndarray) – input in \(L*X=B\)
B (ndarray) – input in \(L*X=B\)
solver_dimension (SolverDimension) – 2d or 3d

Returns

residual scalar

src.functions.mathlib.compute_residual_poisson_tensor_operator(X, B, solver_dimension)

Compute residual matrix/tensor \(r= L*X - B\) using 2d or 3d Laplacian \(L\) in the matrix form (equivalent to \(r=Ax-b\) in the vector form).

Warning

This is based on the matrix form setup, excluding wall boundaries by one cell, and is only valid for inverse Poisson equation. The forward Poisson version needs to be implemented.

Warning

This function is not safe if not excluding wall boundaries. Use this function in subdomains to properly exclude the wall boundaries.

Parameters

X (ndarray) – input in \(L*X=B\)
B (ndarray) – input in \(L*X=B\)
solver_dimension (SolverDimension) – 2d or 3d

Returns

residual matrix/tensor

src.functions.mathlib.compute_residual_subdomain(X, B, sub_shape, solver_dimension)

Compute subdomain L2-norm residual \(r= L*X - B\) using 2d or 3d Laplacian \(L\) in the matrix form (equivalent to \(r=Ax-b\) in the vector form).

Warning

This is based on the matrix form setup, excluding wall boundaries by one cell, and is only valid for inverse Poisson equation. The forward Poisson version needs to be implemented.

Parameters

X (ndarray) – input in \(L*X=B\)
B (ndarray) – input in \(L*X=B\)
sub_shape (2-tuple or 3-tuple) – subdomain shape, 2d or 3d
solver_dimension (SolverDimension) – 2d or 3d

Returns

residual (2norm scalar)

src.functions.mathlib.compute_residual_subdomain_2d(X, B, sub_shape)

Compute 2d subdomain L2-norm residual \(r= L*X - B\) using 2d Laplacian \(L\) in the matrix form (equivalent to \(r=Ax-b\) in the vector form).

Warning

This is based on the matrix form setup, excluding wall boundaries by one cell, and is only valid for inverse Poisson equation. The forward Poisson version needs to be implemented.

Parameters

X (ndarray) – input in \(L*X=B\)
B (ndarray) – input in \(L*X=B\)
sub_shape (2-tuple) – subdomain shape

Returns

residual (2norm scalar)

src.functions.mathlib.compute_residual_subdomain_3d(X, B, sub_shape)

Compute 3d subdomain L2-norm residual \(r= L*X - B\) using 2d Laplacian \(L\) in the matrix form (equivalent to \(r=Ax-b\) in the vector form).

Warning

This is based on the matrix form setup, excluding wall boundaries by one cell, and is only valid for inverse Poisson equation. The forward Poisson version needs to be implemented.

Parameters

X (ndarray) – input in \(L*X=B\)
B (ndarray) – input in \(L*X=B\)
sub_shape (3-tuple) – subdomain shape

Returns

residual (2norm scalar)

src.functions.mathlib.compute_ssim(test_image, ref_image)

Compute SSIM - Structural Similarity Index.

Measure the perceptual similarity and difference between two images, as well as the gradient and the mean SSIM. Check out the formal definition and these examples.

How to compute the SSIM difference:

\(S\) is the local similarity index. It is usually expected to be in [0, 1] but it can also be [-1, 1] where negative values mean the same structure but with inverted values (it is due to a cross product). Since we do not really care about inversion we take the absolute value of \(S\) to make it [0, 1]. The computed difference, working with everything normalized, gives 0 when the two matrices are exactly the same (because \(S\) will be 1), and 1 when they are completely different (\(S\) will be 0).

Suggested vragne for plotting: [-1, 1].

Parameters

test_image (ndarray) –
ref_image (ndarray) –

Returns

ssim_image - structural similarity indices. The full SSIM image. This is only returned if full is set to True
ssim_diff - structural difference
ssim_grad - structural similarity gradient
ssim_mean - mean ssim scalar value

src.functions.mathlib.construct_laplacian_nd_vector_friendly(dim, base_size, positive_definite=False, neumann=False, singularity_value=1)

Construct the Laplacian matrix for the vector from of \(Ax=b\) where \(A\) is the Laplacian.

Note

The system size in \(Ax=b\) is the size of the vectors \(x\) and \(b\) in the vector form, and can be achieved from the size of the Laplacian matrix:

2D : \(\sqrt{dim(A)}\)
3D : \(\sqrt[3]{dim(A)}\)

Parameters

dim (SolverDimension) – D2 or D3
base_size (int) – for a square matrix, number of rows or columns (or depth in 3D). This will be used to compute the Laplacian size. e.g. base_size = 5 : 2D Laplacian \(5^2\), 3D Laplacian \(5^3\)
positive_definite (bool) – just a flag to flip the sign of the matrix to have positive values on the diagonals. Default is negative diagonals (False), to make it consistent with the Poisson kernel generator. In our case where we are just interested in convergence properties this sign flip does not change anything
neumann (bool) – with Neumann boundary condition (Default= False)
singularity_value (float) – regularization factor for treating the ill-conditioned matrix (Default= 1)

Returns

square/cubic Laplacian matrix/tensor in 2D/3D

src.functions.mathlib.convergence_rate(defects)

Compute the convergence for a multi-grid setup. Inspired by this work.

Parameters

defects –

Returns

average convergence rate
convergence rates

src.functions.mathlib.convolve_filter_1d(M, filter_1d, orientation, skip_margin)

Convolve a 2d matrix with a 1d filter.

Parameters

M (ndarray) – 2d input matrix. It does not have to be square
filter_1d (ndarray) – convolutional filter
orientation (ConvolutionOrientation) – horizontal or vertical filter
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Typical value= 0)

Returns

convolved matrix with shrunk size in the given orientation.

src.functions.mathlib.convolve_filter_obj_aware_2d(M, filter_1d, orientation, collision_mask)

Convolve a 2d data domain with a 1d filter of given orientation. Neumann boundary treatment around in-domain solid object.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

Parameters

M (ndarray) – input 2d matrix to be convolved on
filter_1d (ndarray) – convolutional filter
orientation (ConvolutionOrientation) – filter orientation, horizontal or vertical
collision_mask (ndarray) – 2d object solid mask as in-domain collider

Returns

convolved matrix with shrunk size in the given orientation

src.functions.mathlib.convolve_filter_wall_aware_2d(M, filter_1d, orientation)

Convolve a 2d data domain with a 1d filter of given orientation. Neumann boundary treatment around domain walls.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

Parameters

M (ndarray) – input 2d matrix to be convolved on
filter_1d (ndarray) – convolutional filter
orientation (ConvolutionOrientation) – filter orientation, horizontal or vertical

Returns

convolved matrix with shrunk size in the given orientation

src.functions.mathlib.convolve_filter_wall_aware_3d(M, filter_1d, orientation)

Convolve a 3d data domain with a 1d filter of given orientation. Neumann boundary treatment around domain walls.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

Parameters

M (ndarray) – input 3d matrix (tensor) to be convolved on
filter_1d (ndarray) – convolutional filter
orientation (ConvolutionOrientation) – filter orientation, horizontal, vertical, or fiber (depth)

Returns

convolved tensor with shrunk size in the given orientation

src.functions.mathlib.convolve_kernel_2d(M, kernel, skip_margin=0)

Convolve the input 2d matrix with a kernel.

Parameters

M (ndarray) – input matrix. It does not have to be square.
kernel (ndarray) – convolutional kernel
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Default= 0)

Returns

convolved matrix with shrunk size.

src.functions.mathlib.convolve_kernel_single_pass_2d(M, kernel)

Convolve the kernel with a matrix of equal size to give one scalar output. No sliding windows is used.

Parameters

M (ndarray) – input matrix
kernel (ndarray) – convolutional kernel matrix

Returns

scalar output

src.functions.mathlib.expand_with_padding(M, pad_size, pad_value, dim, opts_boundary=None)

Add equal padding to each side of the matrix with a fixed value (wall padding).

Parameters

M (ndarray) – 2d or 3d matrix
pad_size (int) –
pad_value (int) –
dim (SolverDimension) – 2D or 3D
opts_boundary (OptionsBoundary) – Optional. If available, enforce padding on specific walls only.

Returns

new padded matrix if pad_size != 0, else return the original matrix

src.functions.mathlib.expand_with_padding_2d(M, pad_size, pad_value, opts_boundary_detailed_wall=None)

Add equal padding to each side of the matrix with a fixed value (wall padding).

Parameters

M (ndarray) – 2d matrix
pad_size (int) –
pad_value (int) –
opts_boundary_detailed_wall (OptionsBoundary) – Optional. If available, enforce padding on specific walls only.

Returns

new padded matrix if pad_size != 0, else return the original matrix

src.functions.mathlib.expand_with_padding_3d(M, pad_size, pad_value, opts_boundary_detailed_wall=None)

Add equal padding to each side of the matrix with a fixed value (wall padding).

Parameters

M (ndarray) – 3d matrix
pad_size (int) –
pad_value (int) –
opts_boundary_detailed_wall (OptionsBoundary) – Optional. If available, enforce padding on specific walls only.

Returns

new padded matrix if pad_size != 0, else return the original matrix

src.functions.mathlib.extract_frame_values(M, skip_margin, thickness)

Extract the elements of the outer frame of the 2d matrix.

The rest of the matrix will be zero. The frame starts from the skip_margin.

Parameters

M (ndarray) – 2d input matrix to extract from
skip_margin (int) – offset from sides
thickness (int) – thickness of the frame in terms of number of elements

Returns

new 2d matrix with the extracted frame values.

src.functions.mathlib.extract_from_center_2d(M, sub_shape)

Extract a sub-matrix from the center.

Warning

Both the input matrix and the extraction shape must indicate a symmetrical shape with odd number for rows and columns.

Parameters

M (ndarray) – input 2d matrix
sub_shape (tuple) – tuple (rows, cols) of the sub matrix

Returns

sub-matrix

src.functions.mathlib.extract_from_center_3d(M, sub_shape)

Extract a sub-matrix from the center.

Warning

Both the input matrix and the extraction shape must indicate a symmetrical shape with odd number for rows and columns.

Parameters

M (ndarray) – input 3d matrix
sub_shape (tuple) – tuple (rows, cols, depth) of the sub matrix

Returns

sub-matrix

src.functions.mathlib.find_max_rank_and_filter_size(ranked_filters, safe_rank)

Return maximum rank and filter size required based on excluding the zero-out elements of their values set during adaptive truncation.

Note

Assuming filters are already adaptively truncated, meaning their small values are set to zero given a fixed truncation threshold. All filters have the same size, preserving their original size.

Parameters

ranked_filters (ndarray) – 2d array of the shape (ranks, filters)
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel.

Returns

max_effective_rank : maximum rank required
max_filter_size: maximum filter size required (max of all ranked filters)
actual_filter_size: filter size if there was no truncation

src.functions.mathlib.frobenius_norm(M)

Frobenius norm of the matrix.

Parameters: M (ndarray) – input matrix
Returns: error scalar

src.functions.mathlib.generate_bivariate_gaussian(size, ef_rad, center=None)

Generate a normalized 2d gaussian with maximum 1 and minimum 0.

The effective radius (ef_rad) will be automatically computed to get \(\mu_x\) and \(\mu_y\).

Parameters

size (int) – size of the 2d matrix
ef_rad (float) – effective radius, full-width-half-maximum
center (ndarray) – if None, the center of the matrix, else movable center (Default= None)

Returns

2d Gaussian

src.functions.mathlib.get_1d_array_half_size(length, trunc_factor_percent, enforce_min_size=True)

Given 1d array length and a truncation factor (%), compute the new half array size. Half size excludes the element in the middle. Minimum filter/kernel size should be 3, a half size of 0 means a convolution of a 1x1 element.

Parameters

length (int) – size of the 1d array
trunc_factor_percent (float) – truncation factor in [0, 1]
enforce_min_size (bool) – this is specially used for convolutional stuff. If True, minimum filter/kernel size should be 3, which is achieved by avoiding a half size of 0. This is necessary to remain consistent with the Jacobi 3x3 base kernel.

Returns

half_size - new half size
middle_index - the index to the middle element

src.functions.mathlib.get_effective_nonzero_1darray(arr)

Find where non-zero elements are.

Parameters: arr – input array
Returns: array with only nonzero elements

src.functions.mathlib.get_kernel_effective_size(opts)

Get the actual convolution kernel size, whether it is full or reduced.

Parameters: opts (OptionsGeneral) –
Returns: effective kernel size

src.functions.mathlib.get_kernel_trimmed_size(kernel_size, trunc_factor_percent)

Compute a new kernel size based on a truncation factor % and the original kernel size (assume a square kernel).

Parameters

kernel_size (int) – original kernel size
trunc_factor_percent (float) – truncation factor in [0, 1]

Returns

new kernel size

src.functions.mathlib.get_sub_domain_shape(data_original, num_exclude_cell, dim)

Compute the size of the subdomain data given the numbers of cells to exclude from each side.

Parameters

data_original (ndarray) – original data matrix, either 2d or 3d
num_exclude_cell (int) – number of cells to exclude from each side
dim (SolverDimension) – 2D or 3D

Returns

2d or 3d shape

Return type

tuple

src.functions.mathlib.get_truncate_indices(arr, cut_off)

Find indices of array values larger than or equal to the cut_off value.

We get the indices that should be kept. For a symmetrical array with larger values around the center, this means an index list of the sub-array spanning a range around the center.

Parameters

arr (ndarray) – input 1darray
cut_off (float) – truncation threshold.

Returns

cut - the exact truncation index, cut_indices: all indices that should be kept
cut_indices - the index list of all cut values

Warning

Returns None for both if the whole array is subject to truncation (nothing will be left).

src.functions.mathlib.inf_norm(M)

Infinite norm of the matrix.

Parameters: M (ndarray) – input matrix
Returns: error scalar

src.functions.mathlib.interp_1darray(arr, resolution, symmetric_x_axis=True, kind='cubic')

Interpolate a 1d array.

Parameters

arr (1darray) – input
resolution (int) – 1 no interpolation, 1> increased resolution
kind (str) – interpolation function (Default=cubic); see scipy.interpolate for options.
symmetric_x_axis (bool) – if True add points equally to each side of the array center

Returns

up-sampled indices and array

src.functions.mathlib.is_active_domain(mask, i, j)

Check if a certain position of the mask is active (data) domain.

The marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

mask (ndarray) – mask matrix
i (int) – row index
j (int) – column index

Returns

True if the position is part of the active (data) domain

Return type

bool

src.functions.mathlib.is_contour_2d(mask, i, j)

Check if a certain position of the mask is contour.

The marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

mask (ndarray) – mask matrix
i (int) – row index
j (int) – column index

Returns

True if the position is part of the contour

Return type

bool

src.functions.mathlib.is_solid_2d(mask, i, j)

Check if a certain position of a mask is solid.

The marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

mask (ndarray) – mask matrix
i (int) – row index
j (int) – column index

Returns

True if it is solid or contour

Return type

bool

src.functions.mathlib.is_wall_2d(M, i, j)

Check if a certain position of a 2d domain matrix falls in walls.

Walls are 1-cell wide paddings around the domain.

Parameters

M (ndarray) – domain matrix
i (int) – row index
j (int) – column index

Returns

True if it is inside the wall

Return type

bool

src.functions.mathlib.is_wall_3d(M, i, j, k)

Check if a certain position of a 3d domain matrix falls in walls.

Walls are 1-cell wide paddings around the domain.

Parameters

M (3darray) – domain matrix
i (int) – row index
j (int) – column index
k (int) – depth index

Returns

True if it is inside the wall

Return type

bool

src.functions.mathlib.make_4_tiles_mirrored_2d(M, is_odd)

Given the input matrix make even or odd tile block arrangement where the top left block is the

See make_4_tiles_mirrored_odd_2d() and make_4_tiles_mirrored_even_2d() for more details.

Parameters

M (ndarray) – input 2d matrix
is_odd (bool) – even or odd arrangement

Returns

even or odd mirrored matrix

src.functions.mathlib.make_4_tiles_mirrored_even_2d(M)

Given the input matrix make 4 tiles block even arrangement where the top left block is the input, the top right block is the mirror of the first block and the two bottom blocks are the mirrors of the top blocks.

For example, given

\[\begin{split}\begin{bmatrix} 1 & 2 \\ a & b \end{bmatrix}\end{split}\]

its 4-tiles even structure looks like this:

\[\begin{split}\begin{bmatrix} \mathbf{1} & \mathbf{2} & 2 & 1 \\ \mathbf{a} & \mathbf{b} & b & a \\ a & b & b & a \\ 1 & 2 & 2 & 1 \end{bmatrix}\end{split}\]

Parameters: M (ndarray) – input 2d matrix
Returns: even mirrored matrix

src.functions.mathlib.make_4_tiles_mirrored_odd_2d(M)

Given the input matrix make 4 tiles block odd arrangement where the top left block is the input, the top right block is the mirror of the first block and the two bottom blocks are the mirrors of the top blocks.

For example, given

\[\begin{split}\begin{bmatrix} 1 & 2 \\ a & b \end{bmatrix}\end{split}\]

its 4-tiles even structure looks like this:

\[\begin{split}\begin{bmatrix} \mathbf{1} & \mathbf{2} & 1 \\ \mathbf{a} & \mathbf{b} & a \\ 1 & 2 & 1 \end{bmatrix}\end{split}\]

Parameters: M (ndarray) – input 2d matrix
Returns: odd mirrored matrix

src.functions.mathlib.make_9_tiles_mirrored_2d(M)

Given the input matrix make 9 tiles arrangement, where the central block is the input, and each of the remaining 8 tiles are the mirrored of their neighbours. Each edge separating two tiles acts as the axis the mirrored image is flipped. The output therefore has 3x the size of the original matrix.

For example for a 2x2 matrix

\[\begin{split}\begin{bmatrix} 1 & 2 \\ a & b \end{bmatrix}\end{split}\]

its 9-tiles structure looks like this, with the matrix being the central block:

\[\begin{split}\begin{bmatrix} b & a & a & b & b & a \\ 2 & 1 & 1 & 2 & 2 & 1 \\ 2 & 1 & \mathbf{1} & \mathbf{2} & 2 & 1 \\ b & a & \mathbf{a} & \mathbf{b} & b & a \\ b & a & a & b & b & a \\ 2 & 1 & 1 & 2 & 2 & 1 \end{bmatrix}\end{split}\]

Parameters: M (ndarray) – input 2d matrix
Returns: 9-tiles matrix, with 3 times the size of the input matrix

src.functions.mathlib.mark_contour_2d(mask, i, j)

Mark a certain position of the mask to be contour.

The marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

mask (ndarray) – mask matrix
i (int) – row index
j (int) – column index

src.functions.mathlib.mark_solid_2d(mask, i, j)

Mark a certain position of the mask to be solid.

The marking convention:

data domain = 0
solid = 1
contour = 2 (still part of solid, but with a special flag)

Parameters

mask (ndarray) – mask matrix
i (int) – row index
j (int) – column index

src.functions.mathlib.mre(x, y=None, abs_err=False, axis=None)

Mean relative error of two matrices. If only one matrix is given, returns the average of L1 norm of x.

Parameters

x (ndarray) – input array/matrix
y (ndarray) – input array/matrix, optional (Default= None)
abs_err (bool) – if True use the absolute values of the inputs to compute the error (Default= False)
axis (int, 2-tuple of ints) – specifies the axis along which to compute the vector norms (Default= None)

Returns

error scalar

src.functions.mathlib.ms_norm(M)

Mean squared norm of the matrix. Element-wise.

Parameters: M (ndarray) – input matrix
Returns: error scalar

src.functions.mathlib.mse(x, y=None, axis=None)

Mean squared error of two matrices. If only one is given, returns the L2-norm of x.

Parameters

x (ndarray) – input array/matrix
y (ndarray) – input array/matrix, optional (Default= None)
axis (int, 2-tuple of ints) – specifies the axis along which to compute the vector norms (Default= None)

Returns

error scalar

src.functions.mathlib.normalize_range(M, symmetric=False)

Normalize the matrix values.

Parameters

M (ndarray) – input matrix
symmetric (bool) – [-1, +1] if True, else [0, +1], optional (Default= True)

Returns

normalized matrix

src.functions.mathlib.outer_product_2d(a, b)

Outer product of two vectors.

Parameters

a – first vertical vector
b – second horizontal vector

Returns

\(a \times b\)

src.functions.mathlib.post_process_obj_boundary_enforcement_2d(M, collision_mask, opts)

Boundary enforcement to make sure domain values give the right edge gradient (Neumann). Usually used as an optional last step clean up after the Poisson equation solve.

Note

Complex Object Boundary

This is not perfect. For corner cells on the object there is a race condition as what value we should use to copy from. This is because there are more than one cell in the active domain to copy from. A remedy would be to take average. Nonetheless, because we don’t really compute the gradient, and we are just interested in the active domain values, we skip improving this step.

This is for the collocated grid.

For MACGrid we do not really have to worry about the interior cells inside the object when solving, for instance, for the pressure. The pressure gradient on the edge cells is already ensured to be zero during the Jacobi iterations, so we do not have to explicitly set the object cell values.

Parameters

M (ndarray) – input 2d matrix
collision_mask (ndarray) – 2d object solid mask as in-domain collider
opts (OptionsGeneral) – general options (contains boundary enforcement type, but only supports Neumann for now)

Returns

M treated domain

src.functions.mathlib.post_process_wall_boundary_enforcement_2d(M, opts)

Boundary enforcement to make sure domain values give the right edge gradient (Neumann). Usually used as an optional last step clean up after the Poisson equation solve. Walls only.

Parameters

M (ndarray) – input 2d matrix
opts (OptionsGeneral) – general options (contains boundary enforcement type, but only supports Neumann for now)

Returns

M treated domain

src.functions.mathlib.post_process_wall_boundary_enforcement_3d(M, opts)

Boundary enforcement to make sure domain values give the right edge gradient (Neumann). Usually used as an optional last step clean up after the Poisson equation solve. Walls only.

Parameters

M (ndarray) – input 3d matrix
opts (OptionsGeneral) – general options (contains boundary enforcement type, but only supports Neumann for now)

Returns

M treated domain

src.functions.mathlib.set_frame_boundary_from_matrix(M, thickness, source)

Set the boundary values of the given matrix from another matrix.

The extracted boundary elements can have larger than one-cell thickness.

Parameters

M (ndarray) – 2d target matrix to edit
thickness (int) – thickness of the boundary
source (ndarray) – 2d source matrix to extract the boundary values from

Returns

modified matrix

src.functions.mathlib.set_frame_fixed_value(M, skip_margin, thickness, value, left=True, right=True, top=True, bottom=True)

Set the elements of the outer frame of the 2d matrix to the given value.

The frame starts from the skip_margin. The values in the skip_margin and inside the matrix (excluding frame) will remain intact.

Parameters

M (ndarray) – input 2d matrix
skip_margin (int) – offset from sides
thickness (int) – thickness of the frame in terms of number of elements
value (float) – the values inside the frame (fixed for all)
left (bool) – apply the left side of the frame
right (bool) – apply the right side of the frame
top (bool) – apply the top side of the frame
bottom (bool) – apply the bottom side of the frame

Returns

overwritten input matrix with the frame values

src.functions.mathlib.set_obj_bound_2d(M, here_index, collision_mask, bound_type)

Boundary condition enforcement on a specific cell of the solid object inside the domain.

Parameters

M (ndarray) – input 2d matrix
here_index (Vector2DInt) – cell index
collision_mask (ndarray) – 2d object solid mask as in-domain collider
bound_type (BoundaryType) – boundary enforcement type. Currently only supports Neumann on cell edge

Returns

M treated domain

src.functions.mathlib.set_obj_bound_neumann_edge_2d(M, collision_mask, here_index)

Enforce Pure Neumann boundary condition on a specific cell of the solid object inside the domain.

Set the wall values from the adjacent interior domain cell. Zero gradient on the edge, i.e. the boundary is defined on the edge separating the interior and the exterior cells.

If current cell is obstacle, do nothing. If neighbour cell is wall, consider that the lateral solid value has the same value as ‘here’ to enforce Neumann boundary (\(dp/dn=0\)).

Note

Complex Object Boundary

This is not perfect. For corner cells on the object there is a race condition as what value we should use to copy from. This is because there are more than one cell in the active domain to copy from. A remedy would be to take average. Nonetheless, because we don’t really compute the gradient, and we are just interested in the active domain values, we skip improving this step.

This is for the collocated grid.

For MACGrid we do not really have to worry about the interior cells inside the object when solving, for instance, for the pressure. The pressure gradient on the edge cells is already ensured to be zero during the Jacobi iterations, so we do not have to explicitly set the object cell values.

Parameters

M (ndarray) – input 2d matrix
collision_mask (ndarray) – 2d object solid mask as in-domain collider
here_index (Vector2DInt) – cell index

Returns

M treated domain

src.functions.mathlib.set_padding_2d(M, padding_size, padding_value)

Set the padding region of the input 2d matrix with the padding value.

Parameters

M (ndarray) – input 2d matrix
padding_size (int) – thickness of the padding region
padding_value (int) – const value to be put in the padding region

Returns

padded matrix

Return type

ndarray

src.functions.mathlib.set_padding_3d(M, padding_size, padding_value)

Set the padding region of the input 3d matrix with the padding value.

Parameters

M (ndarray) – input 3d matrix
padding_size (int) – thickness of the padding region
padding_value (int) – const value to be put in the padding region

Returns

padded matrix

Return type

ndarray

src.functions.mathlib.set_single_cell_wall_bound_2d(M, here_index, bound_type)

Boundary condition enforcement on a specific cell of the domain walls.

Parameters

M (ndarray) – input 2d matrix
here_index (Vector2DInt) – cell index
bound_type (BoundaryType) – boundary enforcement type. Currently only supports Neumann on cell edge

Returns

M treated domain

src.functions.mathlib.set_single_cell_wall_bound_3d(M, here_index, bound_type)

Boundary condition enforcement on a specific cell of the domain walls.

Parameters

M (ndarray) – input 3d matrix
here_index (Vector3DInt) – cell index
bound_type (BoundaryType) – boundary enforcement type. Currently only supports Neumann on cell edge

Returns

M treated domain

src.functions.mathlib.set_single_cell_wall_neumann_edge_2d(M, here_index)

Enforce Pure Neumann boundary condition on a specific cell of domain walls.

Set the wall values from the adjacent interior domain cell. Zero gradient on the edge, i.e. the boundary is defined on the edge separating the interior and the wall cells.

If current cell is wall, do nothing. If neighbour cell is obstacle, consider that the lateral solid value has the same value as ‘here’ to enforce Neumann boundary (\(dp/dn=0\)).

e.g. when left wall is between \(x_0\) and \(x_1\), then we impose \(x_0 = x_1\), where \(x_0\) is inside the wall and \(x_1\) is inside the domain.

Note

We do not care much about what values corner cells are getting because in the fluids context only the gradient between the wall and the fluid (the cell edges) matters. The gradient between solid wall cells should not affect the projection step as it does not make sense to correct the fluid velocity inside the wall.

Parameters

M (ndarray) – input 2d matrix
here_index (Vector2DInt) – cell index

Returns

M treated domain walls

src.functions.mathlib.set_single_cell_wall_neumann_edge_3d(M, here_index)

Enforce Pure Neumann boundary condition on a specific cell of domain walls.

Set the wall values from the adjacent interior domain cell. Zero gradient on the edge, i.e. the boundary is defined on the edge separating the interior and the wall cells.

If current cell is wall, do nothing. If neighbour cell is obstacle, consider that the lateral solid value has the same value as ‘here’ to enforce Neumann boundary (\(dp/dn=0\)).

e.g. when left wall is between \(x_0\) and \(x_1\), then we impose \(x_0 = x_1\), where \(x_0\) is inside the wall and \(x_1\) is inside the domain.

Note

We do not care much about what values corner cells are getting because in the fluids context only the gradient between the wall and the fluid (the cell edges) matters. The gradient between solid wall cells should not affect the projection step as it does not make sense to correct the fluid velocity inside the wall.

Parameters

M (ndarray) – input 3d matrix
here_index (Vector3DInt) – cell index

Returns

M treated domain walls

src.functions.mathlib.set_submatrix_zero_2d(M, skip_margin)

Set the sub-matrix to zero.

The values inside the skip margin will remain intact. If skip_margin is zero, this means no skip margin, which means set everything to zero.

Parameters

M (ndarray) – 2d input matrix
skip_margin (int) – how many elements from each side should keep their values

Returns

modified matrix

src.functions.mathlib.set_submatrix_zero_3d(M, skip_margin)

Set the sub-matrix to zero.

The values inside the skip margin will remain intact. If skip_margin is zero, this means no skip margin, which means set everything to zero.

Parameters

M (ndarray) – 3d input matrix
skip_margin (int) – how many elements from each side should keep their values

Returns

modified matrix

src.functions.mathlib.set_wall_bound_2d(M, bound_type)

Boundary condition enforcement on the domain walls.

Parameters

M (ndarray) – input 2d matrix
bound_type (BoundaryType) – boundary enforcement type. Currently only supports Neumann.

Returns

M treated domain walls

src.functions.mathlib.set_wall_bound_3d(M, bound_type)

Boundary condition enforcement on the domain walls.

Parameters

M (ndarray) – input 3d matrix
bound_type (BoundaryType) – boundary enforcement type. Currently only supports Neumann.

Returns

M treated domain walls

src.functions.mathlib.set_wall_bound_neumann_center_2d(M)

Enforce Pure Neumann boundary condition on the cell centers.

Set the wall values from the adjacent interior domain cell. Zero gradient in the center of the interior cell. To achieve zero gradient in the center of first interior cell we need to equate the wall cell and the second interior cell.

e.g. when left wall is \(x_1\), then we impose \(x_0 = x_2\), where \(x_1\) is the actual boundary and \(x_0\) is inside the wall.

Note

We do not care much about what values corner cells are getting because in the fluids context only the gradient between the wall and the fluid (the cell edges) matters. The gradient between solid wall cells should not affect the projection step as it does not make sense to correct the fluid velocity inside the wall.

Parameters: M (ndarray) – input 2d matrix
Returns: M treated domain walls

src.functions.mathlib.set_wall_bound_neumann_edge_2d(M)

Enforce Pure Neumann boundary condition on the wall cell edges.

Set the wall values from the adjacent interior domain cell. Zero gradient on the edge, i.e. the boundary is defined on the edge separating the interior and the wall cells.

e.g. when left wall is between \(x_0\) and \(x_1\), then we impose \(x_0 = x_1\), where \(x_0\) is inside the wall and \(x_1\) is inside the domain.

Note

We do not care much about what values corner cells are getting because in the fluids context only the gradient between the wall and the fluid (the cell edges) matters. The gradient between solid wall cells should not affect the projection step as it does not make sense to correct the fluid velocity inside the wall.

Parameters: M (ndarray) – input 2d matrix
Returns: M treated domain walls

src.functions.mathlib.set_wall_bound_neumann_edge_3d(M)

Enforce Pure Neumann boundary condition on the wall cell edges.

Set the wall values from the adjacent interior domain cell. Zero gradient on the edge, i.e. the boundary is defined on the edge separating the interior and the wall cells.

e.g. when left wall is between \(x_0\) and \(x_1\), then we impose \(x_0 = x_1\), where \(x_0\) is inside the wall and \(x_1\) is inside the domain.

Note

We do not care much about what values corner cells are getting because in the fluids context only the gradient between the wall and the fluid (the cell edges) matters. The gradient between solid wall cells should not affect the projection step as it does not make sense to correct the fluid velocity inside the wall.

Parameters: M (ndarray) – input 3d matrix
Returns: M treated domain walls

src.functions.mathlib.solve_jacobi_matrix_form_no_boundary_2d(M, opts, do_residuals=False, is_subdomain_residual=True, subdomain_shape=None)

Solve the Poisson equation with Jacobi in the matrix form for forward and inverse Poisson equations in 2D with no boundary treatment.

This version of Jacobi solver perfectly matches the results of Poisson filters when boundary treatment is ignored, in other words, for infinite domains.

Note

We use a general Jacobi setup with flexible \(\alpha\) and \(\beta\) instead of fixed values. This allows to adjust the weights based on the type of the Poisson equation.

Warning

The computed residual is in the matrix form and is only valid for inverse Poisson setup.

Note

We need one extra cell trimming from sides because of using pure Poisson filters in comparison. We could have easily added the incomplete Laplacian kernels for the corners and the cells next to the wall, but this will make it inconsistent with the Jacobi function used to generate the Poisson kernels.

Parameters

M (ndarray) – input 2d matrix; if inverse setup, \(M=B\) in \(L*X=B\), and if forward setup, \(M=X\) in in \(L*X=B\), with \(B\) being the unknown
opts (OptionsGeneral) – general options
do_residuals (bool) – collect residuals at each iteration (Default= False)
is_subdomain_residual (bool) – compute the residual only for a subdomain (Default= True)
subdomain_shape (tuple) – subdomain size to be used in computing the residuals (Default= None)

Returns

Out - solution
residuals - residual per iteration

src.functions.mathlib.solve_jacobi_single_padding_obj_collision_2d(M, opts, collision_mask=None)

Solve the Poisson equation with Jacobi in the matrix form for forward and inverse Poisson equations in 2D with wall and in-domain boundary treatment of solid objects with complex shapes.

Note

We use a general Jacobi setup with flexible \(\alpha\) and \(\beta\) instead of fixed values. This allows to adjust the weights based on the type of the Poisson equation.

Note

This function only supports Neumann boundary treatment on the cell edges. The walls are single cell padding.

Parameters

M (ndarray) – input 2d matrix; if inverse setup, \(M=B\) in \(L*X=B\), and if forward setup, \(M=X\) in in \(L*X=B\), with \(B\) being the unknown
opts (OptionsGeneral) – general options (contains number of iterations along with many other variables)
collision_mask – 2d object solid mask as in-domain collider

Returns

Out - solution

src.functions.mathlib.solve_jacobi_single_padding_only_wall_2d(M, sub_shape_residual, opts)

Solve the Poisson equation with Jacobi in the matrix form for forward and inverse Poisson equations in 2D with wall Neumann boundary treatment.

Note

We use a general Jacobi setup with flexible \(\alpha\) and \(\beta\) instead of fixed values. This allows to adjust the weights based on the type of the Poisson equation.

Note

This function only supports Neumann boundary treatment on the cell edges. The walls are single cell padding.

Parameters

M (ndarray) – input 2d matrix; if inverse setup, \(M=B\) in \(L*X=B\), and if forward setup, \(M=X\) in in \(L*X=B\), with \(B\) being the unknown
opts (OptionsGeneral) – general options (contains number of iterations along with many other variables)
sub_shape_residual (2-tuple) – subdomain shape used in computing the residual

Returns

Out - solution
residuals - residual per iteration

src.functions.mathlib.solve_jacobi_single_padding_only_wall_3d(M, sub_shape_residual, opts)

Solve the Poisson equation with Jacobi in the matrix form for forward and inverse Poisson equations in 3D with wall Neumann boundary treatment. Walls only.

Note

We use a general Jacobi setup with flexible \(\alpha\) and \(\beta\) instead of fixed values. This allows to adjust the weights based on the type of the Poisson equation.

Note

This function only supports Neumann boundary treatment on the cell edges. The walls are single cell padding.

Parameters

M (ndarray) – input 3d tensor; if inverse setup, \(M=B\) in \(L*X=B\), and if forward setup, \(M=X\) in in \(L*X=B\), with \(B\) being the unknown
opts (OptionsGeneral) – general options (contains number of iterations along with many other variables)
sub_shape_residual (3-tuple) – subdomain shape used in computing the residual

Returns

Out - solution
residuals - residual per iteration

src.functions.mathlib.solve_jacobi_vector_form(A, b, max_iterations, is_subdomain_residual, sub_domain_shape, dim, warm_start=False)

Solving \(Ax=b\) in vector form.

Works for both 2D and 3D as long as proper \(A\) and \(b\) is fed in the vector form.

Parameters

A (ndarray) – n x n (flat Laplacian)
b (ndarray) – n x 1 (rhs)
max_iterations (int) – stop when exceeding max iteration
is_subdomain_residual (bool) – if computing the residual only for a subdomain.
sub_domain_shape (tuple) – the shape of the subdomain if doing subdomain residual. Must be odd size for proper matrix extraction.
dim (SolverDimension) – 2D or 3D, dimension of the problem
warm_start – this is loosely based on just copying \(b\), currently often makes the convergence worse (Default= False)

Returns

x - solution
residuals - residual per iteration

src.functions.mathlib.solve_poisson_full_kernel_2d(M, kernel, skip_margin=0)

Solve the Poisson equation using 2d full Poisson kernel.

Make a convolution pass excluding the boundary, reading it but not updating it. The Poisson kernel is always square shape.

Parameters

M (ndarray) – input scalar field (matrix)
kernel (ndarray) – Poisson square convolutional kernel
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Default= 0)

Returns

solution as scalar field (matrix)

src.functions.mathlib.solve_poisson_separable_filters_2d(M, filter_hor, filter_ver, safe_rank, skip_margin=0)

Poisson’s equation solve in the reduced space using convolution of separable filters with no boundary treatment.

Warning

No boundary condition treatment.

See solve_poisson_separable_filters_obj_aware_2d() and solve_poisson_separable_filters_wall_aware_2d() for the version with boundary treatment.

This uses multi-rank Poisson filters for a given rank. Given the Poisson equation in the matrix form \(L*X=B\), the Poisson kernel \(L\) (in forward setup) and its inverse \(L^{-1}\) (in inverse setup) are already baked into the Poisson filters. Just provide the input data matrix and the corresponding filters matching the formulation setup you are interested in.

Note

Order of convolution:

The convolution order using separable filters in 3D is \(F * M \approx \displaystyle\sum_{r=1}^n f_{v_r} * (f_{h_r} * M)\)

where

\(F\) - Full Poisson kernel (either \(L\) or \(L^{-1}\))
\(M\) - Input data field
\(f_h\) - Horizontal filter extracted from \(V^T\) in \(SVD(F) = U.S.V^T\), with \(S\) values absorbed
\(f_v\) - Vertical filter extracted from \(U\) in \(SVD(F) = U.S.V^T\), with \(S\) values absorbed
double subscript \(_r\) means the filter corresponding the current rank
\(\displaystyle\sum_{r=1}^n\) is multi-rank summation (i.e. modal solutions)

Filters are obtained from Eigen decomposition of \(F\) using Singular Value Decomposition (SVD) in the Canonical Polyadic Decomposition (CPD) view. The convolution order goes from the inner bracket to outer bracket, meaning first we need to convolve \(M\) with the horizontal filter, then convolve the results with the vertical filter.

For multi-rank convolution we have separate and independent convolutions passes on \(M\), then sum up the results. The summation comes from the CPD view in tensor decomposition, which makes it possible to have rank-1 kernel convolutions to get modal solutions taking care of different frequencies in the data domain.

Warning

DO NOT feed input data matrix \(M\) in outer bracket convolution.

ALWAYS use the results of the previous convolution pass to do the next one.

Parameters

M (ndarray) – input 2d matrix to be convolved on
filter_hor (ndarray) – horizontal Poisson filters
filter_ver (ndarray) – vertical Poisson filters
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel. It is the minimum of Poisson kernel \(L\) rows and columns. You can also use src.functions.decompositions.rank_safety_clamp_2d() to set this
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Default= 0)

Returns

the solution to Poisson’s equation

src.functions.mathlib.solve_poisson_separable_filters_from_components_2d(M, U, S, VT, rank, trunc_method, trunc_factor, skip_margin=0)

2D Poisson solve in a reduced space from SVD components. Filters are computed on the fly and truncated before application.

Note

see solve_poisson_separable_filters_2d() for details on filter extraction and application.

Parameters

M (ndarray) – input 2d matrix to be convolved on
U (ndarray) – \(U\) in \(USV^T\)
S (ndarray) – \(S\) in \(USV^T\)
VT (ndarray) – \(V^T\) in \(USV^T\)
rank (int) – desired rank. It will be safely clamped if larger than the input matrix actual rank.
trunc_factor (float) – if the truncation method is PERCENTAGE then a value in [0, 1], else fixed floating point cut off threshold (FIXED_THRESHOLD)
trunc_method (TruncationMode) – PERCENTAGE or FIXED_THRESHOLD (adaptive truncation)
skip_margin (int) – number of lateral elements to skip in the convolution. This helps with saving computation time when having redundant padding (Default= 0)

Returns

the solution to Poisson’s equation

src.functions.mathlib.solve_poisson_separable_filters_obj_aware_2d(M, filter_hor, filter_ver, safe_rank, collision_mask, opts, individual_rank=None)

Poisson’s equation solve in the reduced space using convolution of separable filters with Neumann boundary treatment around in-domain solid object.

Note

See solve_poisson_separable_filters_2d() for explanation of doing Poisson filter convolution.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

Parameters

M (ndarray) – input 2d matrix to be convolved on
filter_hor (ndarray) – horizontal Poisson filters
filter_ver (ndarray) – vertical Poisson filters
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel. It is the minimum of Poisson kernel \(L\) rows and columns. You can also use src.functions.decompositions.rank_safety_clamp_2d() to set this
collision_mask (ndarray) – 2d object solid mask as in-domain collider
opts (OptionsGeneral) – general options
individual_rank (int) – if not None, only return the result of convolving with this single rank

Returns

the solution to Poisson’s equation

src.functions.mathlib.solve_poisson_separable_filters_wall_aware_2d(M, filter_hor, filter_ver, safe_rank, opts)

Poisson’s equation solve in the reduced space using convolution of separable filters with Neumann boundary treatment around the domain walls (no in-domain solid object treatment).

Note

See solve_poisson_separable_filters_2d() for explanation of doing Poisson filter convolution.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

Parameters

M (ndarray) – input 2d matrix to be convolved on
filter_hor (ndarray) – horizontal Poisson filters
filter_ver (ndarray) – vertical Poisson filters
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel. It is the minimum of Poisson kernel \(L\) rows and columns. You can also use src.functions.decompositions.rank_safety_clamp_2d() to set this
opts (OptionsGeneral) – general options

Returns

the solution to Poisson’s equation

src.functions.mathlib.solve_poisson_separable_filters_wall_aware_3d(M, filters_1d, safe_rank, opts)

Poisson’s equation solve in the reduced space using convolution of separable filters with Neumann boundary treatment around the domain walls (no in-domain solid object treatment) - 3D.

We use Mirror Marching algorithm to enforce Neumann boundary condition. See paper.

We use multi-rank Poisson filters for a given rank. Given the Poisson equation in the matrix form \(L*X=B\), the Poisson kernel \(L\) (in forward setup) and its inverse \(L^{-1}\) (in inverse setup) are already baked into the Poisson filters. Just provide the input data matrix and the corresponding filters matching the formulation setup you are interested in.

Note

Order of convolution:

The convolution order using separable filters in 3D is \(F * M \approx \displaystyle\sum_{r=1}^n f_{v_r} * (f_{h_r} * (f_{d_r} * M))\)

where

\(F\) - Full Poisson kernel (either \(L\) or \(L^{-1}\))
\(M\) - Input data field (tensor)
\(f_v\) - Vertical filter
\(f_h\) - Horizontal filter
\(f_d\) - Depth (fiber) filter
double subscript \(_r\) means the filter corresponding the current rank
\(\displaystyle\sum_{r=1}^n\) is multi-rank summation (i.e. modal solutions)

Filters are obtained from Tensor Eigen Decomposition of \(F\) using Symmetric-CP (Symmetric Canonical Polyadic Decomposition). The convolution order goes from the inner brackets to outer brackets, meaning first we need to convolve \(M\) with the fiber filter, then convolve the results with the horizontal and vertical filters.

For multi-rank convolution we have separate and independent convolutions passes on \(M\), then sum up the results. The summation comes from the Canonical Polyadic Decomposition (CPD) view in tensor decomposition, which makes it possible to have rank-1 kernel convolutions to get modal solutions taking care of different frequencies in the data domain.

Warning

DO NOT feed input data matrix \(M\) in outer bracket convolution.

ALWAYS use the results of the previous convolution pass to do the next one.

Warning

Safe Rank

While in 2D ‘safe’ means a rank that does not exceed the actual rank of the kernel, in 3D it is different. Due to the lack of a clear definition of rank in tensor decomposition, we can pretty much ask for any rank in the CP-vew eigen decomposition and always get something that “works”. This vagueness of a proper rank definition in 3D is partially linked to the fact that tensor decomposition is NP-hard.

In this function ‘safe’ simply means the number of ranked filters we would like to include in our convolutions, and use it to help with setting up loops etc.

Parameters

M (ndarray) – input 3d matrix (tensor) to be convolved on
filters_1d (ndarray) – ranked (but not sorted) filters. Same filter is used for all 3 convolutional orientations: horizontal, vertical, and fiber (depth)
safe_rank (int) – the number of ranked filters to be used in convolution
opts (OptionsGeneral) – general options

Returns

the solution to Poisson’s equation

src.functions.mathlib.trim(M, size)

Trim the input 2d matrix, which will shrink equally on each side and dimension.

Parameters

M (ndarray) – input 2d matrix
size (int) – trimming size

Returns

trimmed matrix

src.functions.mathlib.truncate_filters(truncation_method, truncation_value, safe_rank, filters_1d, preserve_shape=True)

Truncate filters either by PERCENTAGE or FIXED_THRESHOLD (adaptive truncation).

Poisson filters are symmetrical, hence truncation automatically implies symmetrically removing values from sides.

Parameters

truncation_method (TruncationMode) – PERCENTAGE or FIXED_THRESHOLD (adaptive truncation)
truncation_value (float) – if the truncation method is PERCENTAGE then a value in [0, 1], else fixed floating point cut off threshold (FIXED_THRESHOLD)
safe_rank (int) – desired input rank. ‘safe’ means a rank that does not exceed the actual rank of the kernel
filters_1d (ndarray) –
preserve_shape (bool) – if True keep the original shape and fill them with zeros (Default= True)

Returns

a new smaller array with truncated elements

src.functions.mathlib.truncate_fixed_threshold_1darray(arr, cut_off, preserve_shape=True)

Adaptive truncation. Cut any elements smaller than a fixed absolute value, assuming a symmetrical filter and values are sorted sideways from largest (center) to smallest (tales).

Only works with symmetrical filters, and only cuts the outer parts of the filters.

[….outer part…. cut || …….inner part……. || cut ….outer part….]

Cut elements < cut_off

Keep elements >= cut_off

Checking the absolute values against the absolute threshold value.

Parameters

arr (ndarray) – input 1d array
cut_off (float) – truncation threshold (absolute value)
preserve_shape (bool) – if True keep the original shape and fill them with zeros (Default= True)

Returns

truncated (smaller) array. If preserving the shape, keep the shape and insert zeros in the truncated parts.

src.functions.mathlib.truncate_percent_filter_1d(arr, trunc_percent)

Truncate a percentage of a 1d array symmetrically.

It finds the middle, then throws out the truncation % from both ends. The results will be exactly symmetrical for odd sized arrays. For arrays with an even size the middle will be the ceiling (the larger index of the pair in the middle). For a 100% truncation the output is the middle element.

If truncation results in fractional elimination of an element, we assume ceiling and still eliminate that element.

Parameters

arr (ndarray) – input 1d array
trunc_percent (float) – truncation factor in [0, 1]

Returns

a new smaller array with truncated elements

src.functions.mathlib.zero_out_small_values(M, threshold=1e-09)

Zero out all the elements whose absolute values are below the threshold

Parameters

M (ndarray) – input array
threshold (float) – non-negative

Returns

same as input with zeros in place of small values