skshape.numerics package

This module includes numerical integration procedures to approximate the integral of a given function on the given mesh. The functions can be integrated using fixed Gaussian quadrature on the elements of the mesh, or using adaptive integration either by refining mesh elements or increasing order of quadrature based on error criteria.

skshape.numerics.integration.L2_error(f1, f2, mesh, quadrature_degree=None)[source]

L2 error between f1, f2 on mesh using fixed Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the L2 error between f1 and f2 on mesh numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the maximum of the degree of the functions f1 and f2 is used if their degrees are available, otherwise quadrature degree is 2.

Parameters

f1 (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
f2 (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

L2_error_est – The estimated L2 error between f1 and f2 on given mesh.

Return type

float

skshape.numerics.integration.L2_norm(f, mesh, quadrature_degree=None)[source]

L2 norm of f on mesh using fixed Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the L2 norm of f on mesh numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the degree of the function f is used if available, otherwise quadrature degree is 2.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

L2_norm_est – The estimated L2 norm of f on given mesh.

Return type

float

class skshape.numerics.integration.Quadrature[source]

Bases: object

degree()[source]

max_degree()[source]

max_order()[source]

order()[source]

points()[source]

weights()[source]

class skshape.numerics.integration.Quadrature1d(degree=None, order=None)[source]

Bases: skshape.numerics.integration.Quadrature

integrate(fct)[source]

iterpoints()[source]

max_degree()[source]

max_order()[source]

class skshape.numerics.integration.Quadrature2d(degree=None, order=None)[source]

Bases: skshape.numerics.integration.Quadrature

integrate(fct)[source]

iterpoints()[source]

max_degree()[source]

max_order()[source]

skshape.numerics.integration.adapt_integrate(f, mesh, tol=None)[source]

Adaptively integrate a function on a given mesh by refinement.

This procedure takes a function and a mesh, and estimates the integral of f on mesh by numerical quadrature. This is done adaptively. Quadrature error in each element is estimated by comparing the values of low-order and high-order quadratures. If the total error is higher than tol, the elements with highest error are refined. The refinements continue until total error is below the prescribed tol.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration. Some elements of the mesh might be refined into smaller elements.
tol (float) – A tolerance value for the total error in integration. The default value for tol is None. If the user does not specify a tolerance. This function tries to estimate a tol based on minimum scale computed from f.diameter() and f.resolution(). If this fails, tol is set to 100*eps, where eps is the machine epsilon for floating point numbers.

Returns

integral – The computed integral of f on mesh.

Return type

float

skshape.numerics.integration.adapt_integrate_tolerance(f, mesh, quad0=None, quad1=None)[source]

Estimates a minimum tolerance for adaptive integration for (f, mesh).

This procedure takes a function and a mesh, and estimates a minimum tolerance for adaptive integration. The function is used to specify a minimum length scale, which we do not want to overresolve in adaptive integration. For example, if the function is based on an ImageFunction, the ImageFunction provides its resolution and diameter, from which we can compute the pixel size. The variations in function value at a scale of the pixel size determines the minimum tolerance; we should not use many elements with adaptive integration to resolve variations at this scale. The function is sampled with small random elements and the average error is scaled to the mesh size to obtain the tolerance.

Parameters

f (MeshFunction or function_like) – A function of X (a dxn array of n coordinate points in d-space) or an instance of a MeshFunction. It should have member functions: f.resolution(), f.diameter(), both of which are used to figure out the smallest scale.
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), whose ref_element is used to define the random elements.
quad0 (Quadrature) – The less accurate quadrature. Default value is None, in which case a quadrature of order=2 is used.
quad1 (Quadrature) – The more accurate quadrature. Default value is None, in which case a quadrature of order=quad0.order()+1 is used.

Returns

tol – The estimated minimum tolerance for the given functions and mesh.

Return type

float

skshape.numerics.integration.adapt_order_integrate(f, mesh, tol=None)[source]

Adaptively integrate a function by increasing quadrature order.

This procedure takes a function and a mesh, and estimates the integral of f on mesh using numerical quadrature. This is done adaptively. Quadrature error in each element is estimated by comparing the values of low-order and high-order quadratures. If the total error is higher than tol, a higher order quadrature is applied to the elements with highest error in order to integrate those elements with more accuracy. This continues until total error is below the prescribed tol.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration. Some elements of the mesh might be refined into smaller elements.
tol (float) – A tolerance value for the total error in integration.

Returns

integral – The computed integral of f on mesh.

Return type

float

skshape.numerics.integration.double_integrate(f, mesh1, mesh2, quadrature_degree=None)[source]

Double integral of function on two meshes with Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the double integral of f on mesh1 and mesh2 numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the degree of the function f is used if available, otherwise quadrature degree is 2.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh1 (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the first domain of integration.
mesh2 (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the second domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

integral – The computed integral of f on given mesh.

Return type

float

skshape.numerics.integration.integrate(f, mesh, quadrature_degree=None)[source]

Integrate a function on a mesh using fixed Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the integral of f on mesh numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the degree of the function f is used if available, otherwise quadrature degree is 2.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

integral – The computed integral of f on given mesh.

Return type

float

skshape.numerics.integration.max_error(f1, f2, mesh, quadrature_degree=None)[source]

Max error between f1, f2 on mesh using fixed Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the maximum error between f1 and f2 on mesh numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the maximum of the degree of the functions f1 and f2 is used if their degrees are available, otherwise quadrature degree is 2.

Parameters

f1 (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
f2 (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

max_error_est – The estimated L2 error between f1 and f2 on given mesh.

Return type

float

skshape.numerics.integration.max_norm(f, mesh, quadrature_degree=None)[source]

Max norm of f on mesh using fixed Gaussian quadrature.

This procedure takes a function and a mesh, and estimates the maximum norm of f on mesh numerically using Gaussian quadrature. The degree of the quadrature can be specified. If it is not given, then the degree of the function f is used if available, otherwise quadrature degree is 2.

Parameters

f (MeshFunction or function_like) – An instance of MeshFunction or a function of X where X is a dxn array of n coordinate points in d-space).
mesh (Mesh) – An instance of a mesh (possibly Curve, CurveHierarchy, Domain2d), that gives the domain of integration.
quadrature_degree (int, optional) – The degree of Gaussian quadrature to be used for integration.

Returns

max_norm_est – The estimated maximum norm of f on given mesh.

Return type

float

skshape.numerics.marking module

Marking strategies for adaptation of elements based on errors.

This module includes function implementing different marking strategies to decide which elements of a mesh to refine or to coarsen based on the element errors computed.

Summary of marking strategies:

	refinement criterion	coarsening criterion
maximum strategy	error > gamma_r * tol	error < gamma_c * tol
fixed ratio	given ratio of sorted err	no coarsening
fixed el tol	error > tol	error < gamma_c * tol
equidistribution (n = len(error))	error > gamma_r * tol/n	error < gamma_c * tol/n

skshape.numerics.marking.equidistribution_marking(error, parameters, mask=None)[source]

Marking errors above gamma_r*tol/n for refinement (some for coarsening).

Marks positions with error above gamma_r*tol/n for refinement and those with error below gamma_c*tol/n for coarsening, where n is len(error). If gamma_c is not defined, then the array of coarsening indices is empty.

Parameters

error (NumPy array) – A float array of error values.
parameters (dict) – A parameter dictionary of the form ‘’{ ‘tolerance’:tol, ‘gamma_refine’:gamma_r, ‘gamma_coarsen’:gamma_c }’’
mask (NumPy array) – A boolean array of the same size as error or an array of integer indices indicating which errors should be checked. Its default value is None, in which case all the error values are checked.

Returns

refine_indices (Numpy array) – An array of integers denoting the error positions marked for refinement.
coarsen_indices (Numpy array) – An array of integers denoting the error positions marked for coarsening.

skshape.numerics.marking.fixed_el_tol_marking(error, parameters, mask=None)[source]

Marking errors above a fixed tol for marking and some for coarsening.

Marks positions with error above a fixed tol for refinement and those with error below tol*gamma_c for coarsening. If gamma_c is not defined, then the array of coarsening indices is empty.

Parameters

error (NumPy array) – A float array of error values.
parameters (dict) – A parameter dictionary of the form ‘’{ ‘tolerance’:tol, ‘gamma_coarsen’:gamma_c }’’
mask (NumPy array) – A boolean array of the same size as error or an array of integer indices indicating which errors should be checked. Its default value is None, in which case all the error values are checked.

Returns

refine_indices (NumPy array) – An array of integers denoting the error positions marked for refinement.
coarsen_indices (Numpy array) – An array of integers denoting the error positions marked for coarsening.

skshape.numerics.marking.fixed_ratio_marking(error, parameters, mask=None)[source]

Marking a fixed ratio of the errors in descending order.

Marks a fixed ratio of error positions when the errors are sorted in decreasing order, for example, the highest 0.2 (20%) of all errors.

Parameters

error (NumPy array) – A float array of error values.
parameters (dict) – A dictionary of the form ‘’{ ‘tolerance’:tol, ‘ratio’:ratio }’’. Errors < ratio*tol are marked for refinement.
mask (NumPy array) – A boolean array of the same size as error or an array of integer indices indicating which errors should be checked. Its default value is None, in which case all the error values are checked.

Returns

refine_indices (NumPy array) – An array of integers denoting the error positions marked for refinement.
coarsen_indices (NumPy array) – An empty array of integers (none marked for coarsening).

skshape.numerics.marking.maximum_strategy_marking(error, parameters, mask=None)[source]

Marking the indices of the error vector using the maximum strategy.

Marks the indices of the error vector using the maximum strategy:

the indices of errors larger than tol*gamma_r are returned for refinement.
the indices of errors smaller than tol*gamma_c are returned for coarsening.

If gamma_c is not defined, then the array of coarsening indices is empty.

Parameters

error (NumPy array) – A float array of error values.
parameters (dict) – A parameter dictionary of the form ‘’{ ‘tolerance’:tol, ‘gamma_refine’:gamma_r, ‘gamma_coarsen’:gamma_c }’’ tol is the error tolerance, gamma_r is a ratio for refinement, gamma_c is a ratio for coarsening.
mask (NumPy array) – A boolean Numpy array of the same size as error or an array of integer indices indicating which errors should be checked. Its default value is None, in which case all the error values are checked.

Returns

refine_indices (NumPy array) – An array of integers denoting the error positions marked for refinement.
coarsen_indices (NumPy array) – An array of integers denoting the error positions marked for coarsening.

skshape.numerics.matrix module

class skshape.numerics.matrix.BlockDiagMatrix(*args, **kwargs)[source]

Bases: skshape.numerics.matrix.Matrix

copy()[source]

factorize()[source]

rmatvec(v)[source]

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters: x ({matrix, ndarray}) – An array with shape (M,) or (M,1).
Returns: y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.
Return type: {matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

size()[source]

solve(v, factorization=None)[source]

submatrices()[source]

to_sparse(format='lil')[source]

class skshape.numerics.matrix.Matrix(*args, **kwargs)[source]

Bases: scipy.sparse.linalg.interface.LinearOperator

copy()[source]

matvec(v)[source]

Matrix-vector multiplication.

Performs the operation y=A*x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters: x ({matrix, ndarray}) – An array with shape (N,) or (N,1).
Returns: y – A matrix or ndarray with shape (M,) or (M,1) depending on the type and shape of the x argument.
Return type: {matrix, ndarray}

Notes

This matvec wraps the user-specified matvec routine or overridden _matvec method to ensure that y has the correct shape and type.

norm()[source]

rmatvec(v)[source]

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters: x ({matrix, ndarray}) – An array with shape (M,) or (M,1).
Returns: y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.
Return type: {matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

size()[source]

solve(d)[source]

class skshape.numerics.matrix.MatrixWithRankOneUpdates(*args, **kwargs)[source]

Bases: skshape.numerics.matrix.Matrix

copy()[source]

core_matrix()[source]

rmatvec(x)[source]

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters: x ({matrix, ndarray}) – An array with shape (M,) or (M,1).
Returns: y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.
Return type: {matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

size()[source]

solve(b)[source]

update_vectors(copy=True)[source]

class skshape.numerics.matrix.TridiagMatrix(*args, **kwargs)[source]

Bases: skshape.numerics.matrix.Matrix

The TridiagMatrix class to store tridiagonal matrices of the form

( b(0) c(0)             a(0) )
( a(1) b(1) c(1)             )
(      a(2)  .  .            )
(            .  .  .         )
(               .  .  c(n-1) )
( c(n)            a(n) b(n)  )

It is used for implementing efficient linear algebra for tridiagonal matrices.

The following methods are available:

set_diag(d, band_no)[source]

get_diag(d, band_no=None)[source]

rmatvec(v)[source]

factorize()[source]

solve(d, L=None, M=None)[source]

to_sparse(format='lil')[source]

norm(norm_type='Frobenius')[source]

copy()[source]

size()[source]

copy()[source]

factorize()[source]

get_diag(band_no=None)[source]

Returns the selected diagonal of the tridiag matrix

Tridiag.getDiag( d, band_no ) returns the values in the band specified with band_no, which should be one of -1,0,1. The default value of band_no is None. In that case, all three diagonals are returned.

norm(norm_type='Frobenius')[source]

Calculates the norm of the matrix.

The norm_type should be one 0, 2, inf or ‘Frobenius’ (the default value).

rmatvec(v)[source]

Adjoint matrix-vector multiplication.

Performs the operation y = A^H * x where A is an MxN linear operator and x is a column vector or 1-d array.

Parameters: x ({matrix, ndarray}) – An array with shape (M,) or (M,1).
Returns: y – A matrix or ndarray with shape (N,) or (N,1) depending on the type and shape of the x argument.
Return type: {matrix, ndarray}

Notes

This rmatvec wraps the user-specified rmatvec routine or overridden _rmatvec method to ensure that y has the correct shape and type.

set_diag(d, band_no)[source]

Set the diagonal band to given d array.

Tridiag.setDiag( d, band_no ) sets the band specified with band_no to the values given in d. The argument band_no should be one of -1,0,1.

size()[source]: Tridiag.size() returns the size of the tridiagonal matrix

solve(d, L=None, M=None)[source]

to_sparse(format='lil')[source]

skshape.numerics.matrix.fast_tridiag_solve(A, d)[source]

skshape.numerics.matrix.tridiag_LU(a, b, c)[source]

Computes the LU decomposition of the triadiagonal matrix.

tridiag_LU( a, b, c ) computes the LU decomposition of the tridiagonal matrix

    ( b(0) c(0)             a(0) )
    ( a(1) b(1) c(1)             )
A = (      a(2)  .  .            )
    (            .  .  .         )
    (               .  .  c(n-1) )
    ( c(n)            a(n) b(n)  )

from the given three bands a,b,c and returns two vectors l,m that give the LU decomposition

A = L U

where

L = Id + diag(-1,l), U = diag(0,m) + diag(1,r)

skshape.numerics.matrix.tridiag_solve(A, d)[source]

Efficient solve of a tridiagonal linear system.

tridiag_solve( A, d ) takes a tridiagonal matrix A and right hand side vector d and returns the solution x of the linear system A x = d.

skshape.numerics.matrix.tridiag_solve_from_LU(l, m, c, d)[source]

Computes the solution of the linear system from LU decomposition.

tridiag_solve_from_LU( l, m, c, d ) takes the LU decomposition l,m (which should be calculated by tridiagLU), the third band c (band_no = 1), and right hand side d of the tridiagonal linear system and returns its solution. See documentation of tridiagLU for more informations on the arguments.

skshape.numerics package

Subpackages

Submodules

skshape.numerics.function module

skshape.numerics.integration module

skshape.numerics.marking module

skshape.numerics.matrix module

Module contents

Table of Contents