geoml.math
Arrays in, arrays out: the geometry the containers and meshes are built on, and the TensorFlow helpers the models use. Nothing here holds a container.
Geometry
Rotations and angles, the triangulated-surface predicates behind
Surface3D and the assignments, the sub-block lattice
arithmetic, and cell declustering.
- geoml.math.geometry.fan_triangulation(faces)[source]
Splits faces, given as vertex indices, into triangles.
A DXF 3DFACE has four corners and a MESH face may have more, neither of which a Surface3D has room for. A fan from the first corner is the standard split, and is exact for the convex faces a triangulated surface is made of.
- Parameters:
faces (sequence) – One sequence of vertex indices per face, of any length.
- Returns:
triangles (array) – An (n, 3) array of vertex indices.
- geoml.math.geometry.vertex_normals(points, triangles)[source]
The unit normal at each vertex, from the triangles meeting there.
A DXF file carries no normals and Surface3D keeps one per vertex, as marching_cubes hands them over. The cross product of a triangle’s edges has twice the triangle’s area for its length, so summing the face normals before normalizing weights each one by its area — which keeps a large face from being outvoted by the slivers around it.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
- Returns:
normals (array) – An (n, 3) array of unit vectors, one per vertex.
- geoml.math.geometry.weld(points, triangles, precision=6)[source]
Merges vertices sitting at the same place, remapping the triangles.
Whether a surface is closed is a question about its edges, and an edge is only shared if the triangles meeting along it say so with the same two indices. Plenty of meshes are closed in space while indexing every triangle’s corners separately — pyvista.Cylinder is one — and welding is what lets the seams be seen for what they are.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are matched to.
- Returns:
points (array) – The distinct vertices.
triangles (array) – The triangles, indexing into them.
- geoml.math.geometry.drop_degenerate_faces(points, triangles, precision=6)[source]
Removes the faces of a triangulation that bound nothing.
Two kinds go: a triangle whose corners are not three distinct places, which has no area, and a face carrying a twin, which encloses no volume with it. A contour of an unstructured grid emits both wherever the surface passes exactly through a cell corner — the level set pinches to a point there, and the marching cubes case that covers it writes the slivers out anyway.
They matter because they are read as a winding failure. Both make an edge appear twice the same way round, so reversed_edges counts them and mesh3d returns a plain Mesh3D for a surface that is closed and consistent everywhere it has area.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are welded to, as in weld.
- Returns:
points (array) – The vertices, exactly as they came.
triangles (array) – The faces that survived, in the order they came, still indexing the original vertices.
Notes
Coincidence is judged on welded indices and reported on the caller’s, so a vertex list is never reordered or renumbered. That matters: weld sorts, and from_dxf promises a round trip that returns a file’s own vertices in its own order. Vertices left unused by a dropped face stay where they are, costing a row and confusing nothing — every measure that cares welds for itself.
How many copies of a twinned face to drop is worked out on the surface, not predicted. Dropping every copy is right for a zero-thickness flap, whose edges are either its own alone or already carried by the surface it lies against; it is wrong for a wall that happens to be recorded twice, where it would leave each edge bounding one face — an open edge where there is no hole. And the two cannot be told apart one group at a time, because the groups couple: in a doubled patch — a membrane several faces wide, which is what decimation makes of a collapsed thin feature — the middle face sees every neighbour as a twin, and any rule that scores its edges in isolation drops it while its neighbours stay, tearing a hole down the middle of the patch.
So duplicated faces start out dropped, and copies are put back one at a time wherever the surface as it stands is left with an edge on exactly one face — a boundary the mesh did not have. Each pass resurrects at least one group or stops, so the loop is bounded by the number of groups, and each resurrection can heal the edges of the next: the rim of a doubled patch comes back first, and the middle on the pass after, once the rim’s return has left its edges half-supported.
- geoml.math.geometry.split_touching_edges(points, triangles, precision=6)[source]
Splits every edge where two pieces of a surface touch.
A surface touching itself along a line – a level set meeting its own closing cap edge-on, two bodies meeting along an edge – comes out welded as an edge four triangles share. It is neither closed nor consistently wound in any reading, and no repair of the winding settles it. Around such an edge the four triangles alternate in the direction they walk it and bound two wedges of inside between two of outside; each triangle is paired with its neighbour across an inside wedge, and each vertex of the edge gets a copy for every piece that meets there, so the pieces touch rather than share.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices, consistently wound where the pieces are.
precision (int) – Decimal places the coordinates are welded to, as in weld.
- Returns:
points, triangles (arrays) – The very arrays given where no edge is shared by more than two triangles. Otherwise the welded vertices, a copy appended for each extra piece meeting at a vertex, and the triangles indexing them. A copy sits exactly on its original; moving the pieces apart is the caller’s business.
Notes
An edge shared by some other number of triangles, or by four that do not alternate, is left as it is.
- geoml.math.geometry.edge_defects(points, triangles, precision=6)[source]
Both edge counts a mesh is judged on, from one welding.
open_edges and reversed_edges ask two questions of the same welded triangulation, and every mesh is built asking both. Welding is the expensive half — a rounding and a unique over the vertices — so doing it once for the pair is worth the one extra function: measured on 27 656 triangles, the two separately cost 9 ms and 7 ms of which 4 ms was each one’s welding, and together they cost 12 ms.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are welded to.
- Returns:
open_count (int) – Edges belonging to a single triangle: none on a closed mesh.
reversed_count (int) – Edges walked the same way round by both triangles sharing them: none where the winding is consistent.
See also
open_edges,reversed_edges,one
- geoml.math.geometry.open_edges(points, triangles, precision=6)[source]
How many of a surface’s edges belong to a single triangle.
None on a closed body, where every edge is shared by two faces; at least the outline on a sheet. It is what tells the two apart, and so which questions a surface can answer — a body has no elevation above a location, and a sheet has no inside. Vertices are welded first, so a mesh that is closed in space counts as closed however its corners are indexed.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are welded to.
- Returns:
count (int) – The number of edges belonging to one triangle only.
- geoml.math.geometry.reversed_edges(points, triangles, precision=6)[source]
How many edges the triangles sharing them walk the same way round.
None where the winding is consistent: two triangles meeting along an edge traverse it in opposite directions, which is what makes “outward” mean one thing over a whole closed surface. Any at all and some triangle faces the wrong way, which an inside/outside test reads as a hole in the body — quietly, and only in the region the offending faces bound.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are welded to.
- Returns:
count (int) – The number of edges traversed more than once in the same direction.
- geoml.math.geometry.area(points, triangles)[source]
The surface area of a triangulation.
Meaningful whether or not the surface closes, unlike its volume.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
- Returns:
area (float)
- geoml.math.geometry.components(points, triangles, precision=6)[source]
Labels the triangles by the connected piece of surface they belong to.
A boolean operation readily answers with a surface in several pieces — an ore body cut in two, a shell around a cavity — and each piece is a body in its own right. Vertices are welded first, since pieces that touch only through unwelded corners are one piece in space.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
precision (int) – Decimal places the coordinates are welded to.
- Returns:
count (int) – How many pieces there are.
labels (array) – One piece number per triangle.
- geoml.math.geometry.single_valued(points, triangles, tolerance=1e-09)[source]
Whether a surface stands at one height over each (x, y).
True where every triangle projects onto the ground the same way round. A fold or an overhang turns some of them over, and a closed body turns its whole underside over, so both are caught. Triangles standing vertically project to nothing at all and are allowed: a cliff is single valued everywhere except along the line of its face.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
tolerance (float) – Projected areas this much smaller than the largest count as nothing.
- Returns:
single_valued (bool)
- geoml.math.geometry.signed_volume(points, triangles)[source]
The volume a closed surface encloses, negative if it is wound inwards.
Each triangle forms a tetrahedron with a fixed point, whose signed volume is a sixth of the determinant of its corners; over a closed surface those add up to what it encloses, wherever the point happens to be. The point is the vertices’ own centre rather than the origin, whose tetrahedra at mine-grid coordinates are vast and cancel one another down to the answer within rounding. The sign is the useful part: it says which way the triangles face taken together, which is what an inside/outside test must know and cannot learn from any one of them.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices, of a closed surface.
- Returns:
volume (float) – Positive where the triangles face outwards, negative where they face in. Meaningless for a surface that is not closed.
- geoml.math.geometry.sheet_interpolator(points, triangles)[source]
Prepares a sheet to be asked its elevation.
The sheet must be single valued — checking that is the caller’s business, and open_edges is what tells a sheet from a body. matplotlib takes a folded triangulation without complaint, answering with whichever of its sheets it happens to find.
- Parameters:
points (array) – An (n, 3) array of vertex coordinates.
triangles (array) – An (m, 3) array of vertex indices.
- Returns:
interpolator (matplotlib.tri.LinearTriInterpolator) – To be handed to sheet_elevation.
- geoml.math.geometry.sheet_elevation(interpolator, coordinates)[source]
The sheet’s height over each location, NaN past its edge.
- Parameters:
interpolator (matplotlib.tri.LinearTriInterpolator) – From sheet_interpolator.
coordinates (array) – An (n, 2) or (n, 3) array; only the first two columns are read.
- Returns:
elevation (array) – One height per location, NaN where the sheet does not reach.
- geoml.math.geometry.inside_solid(mesh, coordinates)[source]
Whether each location falls within a closed body, asking VTK.
- Parameters:
mesh (pyvista.PolyData) – The body, which must be watertight — checking that is the caller’s business, and open_edges is what tells it.
coordinates (array) – An (n, 3) array of locations.
- Returns:
inside (array) – One boolean per location.
- class geoml.math.geometry.ConcaveHull(points, length)[source]
Bases:
objectThe alpha shape of a point set: a concave hull at a chosen length.
The Delaunay triangulation’s simplices are kept where their circumradius is below length, so the hull follows the data at that scale – it fills the interior between drill fences closer than length apart and leaves out a notch wider than that, where a convex hull would bridge the notch and a ball around each sample would leave the interior out. Built by concave_hull.
- points
The (n_points, n_dim) input.
- Type:
array
- length
The circumradius the simplices were kept under.
- Type:
float
- kept
One boolean per Delaunay simplex.
- Type:
array
- simplices
The kept simplices, as vertex indices into points.
- Type:
array
- boundary
The facets of kept simplices that face an unkept one or nothing – (n_facets, n_dim) vertex indices: segments in 2D, triangles in 3D.
- Type:
array
- geoml.math.geometry.concave_hull(points, length)[source]
The concave hull of a point set at a length scale – an alpha shape.
- Parameters:
points – (n_points, n_dim) sample locations, 2D or 3D.
length – The largest circumradius a Delaunay simplex may have and still be part of the hull: a length in the coordinates’ units, so a hull “at 100 m” spans gaps in the data narrower than that and stops at wider ones.
- Returns:
ConcaveHull – With contains(coordinates) for the inside test and boundary for the facets.
- Raises:
ValueError – If the points do not span their space (fewer than n_dim + 1 of them, or all on a line or plane), which leaves nothing to triangulate.
- geoml.math.geometry.bounding_box(points)[source]
Computes a point set’s bounding box and its diagonal.
- Parameters:
points (array) – A set of coordinates.
- Returns:
bbox (array-like) – Array with the box’s minimum and maximum values in each direction.
d (float) – The box’s diagonal length.
- geoml.math.geometry.declustering_weights(coordinates, values=None, cell=None, origins=4, n_sizes=24)[source]
Cell-declustering weights, one per location.
Samples are rarely laid down evenly. Drilling follows the ore, so the interesting ground is crowded and the rest is sparse, and every statistic that treats the samples as equal votes then describes the sampling rather than the field. Cell declustering is the classical repair: lay a lattice over the data, and split one vote among the samples sharing a cell, so a crowded cell speaks once rather than twenty times.
- Parameters:
coordinates – (n_data, n_dim) sample locations.
values – Sample values, needed only to choose cell automatically.
cell – Cell side. Chosen automatically when absent, which needs values.
origins – How many shifted lattices to average the weights over. Where the lattice starts is arbitrary and on a small sample it moves the answer.
n_sizes – Cell sides tried when choosing one.
- Returns:
weights (array) – (n_data,), summing to n_data, so that a set of evenly spread samples comes back at one apiece.
cell (float) – The size used, whether given or chosen.
Notes
The automatic choice follows the usual practice (Deutsch & Journel’s declus): sweep the cell size and keep the one whose declustered mean departs furthest from the naive one. Both extremes of the sweep return the naive mean – a cell below the sample spacing gives every point its own vote, and one larger than the domain puts them all in one cell – so the departure has an interior maximum, and taking it by absolute value handles clustering in high and in low values alike without being told which happened.
References
Deutsch, C. V., & Journel, A. G. (1998). GSLIB: Geostatistical Software Library and User’s Guide (2nd ed.). Oxford University Press.
- geoml.math.geometry.sub_block_index(discretization)[source]
Which sub-block sits where, as integer counts along each axis.
Axis 0 varies fastest, the order _blockdata has always used and the one the likelihood’s noise is indexed by. Both the sub-block offsets and, in a BlockSet3D, the children of a split are built from this, so sub-block j of a block and child j of that same block are the same corner of it.
- geoml.math.geometry.unit_sub_grid(discretization)[source]
Sub-block offsets from a block’s centre, as fractions of its size.
The same layout _blockdata builds, but divided through by the block so that one array serves every size. Scaling it per block is the whole of what a variable-size block model has to do differently when it fans out.
- geoml.math.geometry.trilinear_weights(discretization)[source]
What each of a block’s eight corners is worth at each sub-block centre.
A corner carries what the blocks meeting there say, so reading the corners at the sub-blocks is how a child learns the shape running across its parent. The layout is symmetric about the centre, so the weights average to an eighth apiece and a correction built from them cancels over the children – which is what keeps a block’s own estimate the mean of the children standing in for it.
- geoml.math.geometry.dyadic_overlaps(origin, size)[source]
The cells of a dyadic tiling that sit inside a larger one.
Cells aligned to their own power-of-two size — origin a multiple of it, per axis — are either disjoint or nested, never partially overlapping, so overlap detection collapses to an ancestor lookup: a cell overlaps iff the cube one of the present sizes above it is itself present. Built for reading foreign octrees, whose files are not ours to trust.
- Parameters:
origin (array) – (n, 3) integer cell origins, each a multiple of its cell’s size.
size (array) – (n,) integer cell sizes, powers of two.
- Returns:
offending (array) – Indices of the cells contained in some larger present cell; empty when the tiling is sound.
- geoml.math.geometry.dyadic_complement(origin, size, shape, coarsest)[source]
The gaps of a dyadic tiling, as the largest aligned cells that fit.
What makes a partial octree full again: a foreign file usually carries only the cells inside a domain of interest, and the always-full design wants the rest present and marked rather than absent. Walks the box top down — a candidate cube that is a present cell, or sits inside one, is covered; one holding a present descendant splits into its eight children; one holding nothing is a gap, emitted at that size.
- Parameters:
origin (arrays) – The present cells, as in dyadic_overlaps — already checked not to overlap, since a nested pair here would be read as coverage.
size (arrays) – The present cells, as in dyadic_overlaps — already checked not to overlap, since a nested pair here would be read as coverage.
shape (array) – (3,) the box to fill, in cells; each axis a multiple of coarsest.
coarsest (int) – The largest cell size to emit, a power of two; the recursion starts on the grid of cubes this size.
- Returns:
gap_origin, gap_size (arrays) – The cells that complete the tiling; empty when it already is one.
- geoml.math.geometry.grow(corners, marked, rings)[source]
Add rings of neighbouring blocks, through the corners blocks share.
Sparse on purpose: dilating a mask over the base lattice would cost a cell for every one the model exists to avoid carrying.
- geoml.math.geometry.point_normals(points, k=12, orient='concave')[source]
Unit normals of a scattered point set lying on a surface.
Each normal is the direction of least spread among a point’s nearest neighbours (Hoppe et al. 1992), which fixes it up to sign. The sign is made consistent over the whole set by propagating it along a minimum spanning tree of the neighbour graph, from one root per connected piece, so that the field can serve as a gradient constraint.
- Parameters:
points – (n, 2) or (n, 3).
k – Neighbours per point.
orient – “concave” (default): the root’s normal points toward the side the surface curves to, read from where its neighbours’ centroid falls relative to the tangent plane; the root is the most curved point of its piece. A surface that reads as flat has no such side, which is warned about, and the sign is then whatever the root drew. A vector instead orients every normal to have a positive component along it.
- Returns:
normals (array) – (n, d) unit normals.
- Raises:
ValueError – In three dimensions, when the points read as a single line, which cannot fix a normal: give the normals, or points off the line.
References
Hoppe, H., DeRose, T., Duchamp, T., McDonald, J. and Stuetzle, W. (1992). Surface reconstruction from unorganized points. SIGGRAPH 1992, 71-78.
Implicit surfaces
A parameter-free radial basis function interpolant through points on a surface and their normals; the surface is its zero level set.
Parameter-free radial basis function interpolants for implicit surfaces.
An implicit surface is the zero level set of a scalar field fitted to points on it, with value zero, and to normals, as its gradient. The field is a polyharmonic radial basis function expansion with a linear drift: no range, no trainable parameter, one linear solve. The Hermite form takes the normals as gradient constraints (Macêdo et al. 2011; Hillier et al. 2014); the off-surface form of Carr et al. (2001) turns each normal into two displaced points instead, which is what lets the thin-plate spline and the linear basis, not twice differentiable at the origin, be used.
In this package’s terms the Hermite system is the posterior mean of a Gaussian process with noise-free value and gradient observations, which is the potential-field method of Lajaunie et al. (1997) written with a conditionally positive definite basis instead of a covariance. None of it is original here; the references are on HermiteRBF.
- geoml.math.rbf.solve_hermite(centres, values, gradient_centres=None, gradients=None, basis='cubic')[source]
The interpolant’s weights, as tensors, for centres given as tensors.
The same symmetric system HermiteRBF solves, assembled and solved in TensorFlow so that it can sit inside a traced graph – what a fault network needs, whose older surfaces are refitted on coordinates the younger faults’ trainable slips restore. Dense, so for centres in the hundreds; the standalone class thins a large set first.
- Parameters:
centres – [n, d] float64 tensor, in the frame the weights will be used in.
values – [n].
gradient_centres – [m, d] locations and the gradient constraints at them, or None for values alone.
gradients – [m, d] locations and the gradient constraints at them, or None for values alone.
basis – As in HermiteRBF; gradient constraints need “cubic”.
- Returns:
alpha, beta, drift – [n], [m, d] (None without gradients) and [d + 1].
- geoml.math.rbf.field(u, centres, alpha, gradient_centres, beta, drift, basis='cubic')[source]
Value and gradient of an interpolant at u, given its weights.
Both [n, d] float64 tensors in the frame of centres; the functional core HermiteRBF evaluates through, exposed for a fault network.
- class geoml.math.rbf.HermiteRBF(points, values=None, normals=None, basis='cubic', transform=None, max_error=None, offset=None, k=12)[source]
Bases:
objectA scalar field through scattered points, with their normals as its gradient: the implicit surface is its zero level set.
Polyharmonic radial basis function interpolant with a linear drift and no free parameter. The cubic basis takes the normals as gradient constraints (the Hermite form). The thin-plate and linear bases are not twice differentiable at the origin, so with them each normal becomes two points displaced along it, at offset and -offset, carrying those values (the off-surface form). Both give a field that is negative behind the normals and positive ahead of them, close to a signed distance near the surface.
- Parameters:
points – The observations, (n, d), on the surface unless values says otherwise.
values – One value per point; zero by default, which is a point on the surface.
normals – The gradient constraints. One of: an (n, d) array aligned with the points, rows of NaN where a point carries none; a pair (locations, vectors) of (m, d) arrays, for normals measured away from the points; None, to derive one for every point by geometry.point_normals, oriented toward the surface’s concavity; or False, to fit values alone. A single normal is enough to pin the field – nothing is required per point.
basis – “cubic” (r^3, the default, Hermite), “thin_plate” (r^2 log r) or “linear” (r), the last two through off-surface points.
transform – A fixed geoml.transform object the coordinates go through before distances are measured, for anisotropy. Given normals are mapped through its Jacobian; derived ones are found in the transformed space directly.
max_error – A geometric error budget, in the units of the field. When given, the centres are chosen greedily, starting from a subset and adding the worst-fitting points until every point is within the budget, so a large redundant set is fitted with a fraction of its points.
offset – The displacement of the off-surface points, for the thin-plate and linear bases; half the median spacing between neighbouring points by default.
k – Neighbours used to derive normals when none are given.
- n_centres
How many points carry the field.
- max_residual
The largest value residual over all points after the fit.
References
Carr, J. C., Beatson, R. K., Cherrie, J. B., Mitchell, T. J., Fright, W. R., McCallum, B. C. and Evans, T. R. (2001). Reconstruction and representation of 3D objects with radial basis functions. SIGGRAPH 2001, 67-76.
Macêdo, I., Gois, J. P. and Velho, L. (2011). Hermite radial basis functions implicits. Computer Graphics Forum 30(1), 27-42.
Hillier, M. J., Schetselaar, E. M., de Kemp, E. A. and Perron, G. (2014). Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions. Mathematical Geosciences 46, 931-953.
Lajaunie, C., Courrioux, G. and Manuel, L. (1997). Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation. Mathematical Geology 29, 571-584.
Duchon, J. (1977). Splines minimizing rotation-invariant semi-norms in Sobolev spaces. In: Constructive Theory of Functions of Several Variables, Springer, 85-100.
Wendland, H. (2005). Scattered Data Approximation. Cambridge University Press.
- property n_centres
- property kept
Indices of the points that carry the field, into the points given.
All of them without max_error; the greedy choice with it, which is what to hand on when a thinned set must be recorded exactly, as a fault transform records its observations.
- to_working(x)[source]
Coordinates of the caller’s frame into the fitted frame, as a tensor: what solve_hermite and field take.
- property gradient_points
Where the gradient constraints sit, (m, d) in the caller’s frame; empty without any.
- property working_gradients
The gradient constraints in the fitted frame, as a (locations, vectors) pair of (m, d) arrays, or None.
- property scale
The factor a gradient in the fitted frame is divided by to be one in the caller’s frame.
TensorFlow helpers
TensorFlow helpers in everyday use. The larger numerical machinery (solvers, Lanczos, Kronecker products) is geoml.math.linalg.
- geoml.math.tf.silence_retracing_notices(silence=True)[source]
Whether to drop TensorFlow’s retracing notice for geoML’s own graphs.
On by default, and installed when geoml is imported. Call it with False to hear them, which is worth doing if a prediction seems to be spending its time compiling rather than computing.
- Parameters:
silence (bool) – True to drop the notices, False to let them through.
See also
geoml.models.GPOptionsjit_predict, the other tracing-related knob.
- geoml.math.tf.pairwise_dist(mat_a, mat_b)[source]
Computes pairwise distances between each elements of matrix and each elements of mat_b.
Args: mat_a, [m,d] matrix mat_b, [n,d] matrix
Returns: dist, [m,n] matrix of pairwise distances
code from https://gist.github.com/mbsariyildiz/34cdc26afb630e8cae079048eef91865
Interpolation
- class geoml.math.interpolate.CubicSpline[source]
Bases:
_CubicSpline- interpolate(x, y, xnew, grad=False)[source]
Optimized cubic spline interpolation.
- Args:
x: [N, B] tensor of knot x-coordinates. Must be sorted. y: [N, B] tensor of knot y-coordinates. xnew: [M, B] tensor of new x-coordinates to interpolate at. grad: bool, if True, returns the derivative dy/dx at xnew.
- invert(x, y, ynew, steps=12)[source]
Solves interpolate(x, y, t) == ynew for t.
Interpolating the knots the other way round – interpolate(y, x, ynew) – is the obvious inverse and is not one: the inverse of a cubic is not a cubic, so the two curves meet at the knots and part between them. This solves the actual polynomial instead, by Newton from that same swapped-spline estimate, which converges quadratically because the estimate is already close.
- Parameters:
x – Knot coordinates, [N, B], with x sorted and y monotone in it.
y – Knot coordinates, [N, B], with x sorted and y monotone in it.
ynew – Values to invert, [M, B].
steps – Newton iterations. Each is a Horner pass over values already gathered – no search, no gather – so the default is generous on purpose. Measured over six samples of lognormal data at knot counts from 11 to 161, the worst error falls as 1e-3 at three steps, 1e-7 at eight and 5e-11 at twelve, and twelve is where it stops improving. Fewer than that is a false economy; more buys nothing.
- Returns:
t, shaped like ynew.
Notes
Accuracy is limited by the transform, not by the iteration. A normal-score fit spends half its intervals nearly flat – 40% of them have a span below 1e-4 even at the default five knots per arm – and where the forward map compresses by a factor s, a residual at machine precision comes back magnified by 1/s. The 5e-11 floor above is exactly that: 2.2e-16 / 1e-5. It is the information the forward map discarded, and no solver recovers it. What this replaced left 1e-1.
The interval is located once, in y: a monotone map puts t in the interval whose index ynew occupies among the y-knots, so the bracket cannot move and no step needs to search again. Each iterate is clipped back into that bracket, which is what keeps the method from wandering where the spline is nearly flat and the Newton step is consequently enormous.
The last correction is taken from a detached iterate, so the derivative reported for this operation is the implicit one, dt/dynew = 1 / f’(t), exactly – rather than whatever differentiating the unrolled iteration would produce.
- class geoml.math.interpolate.MonotonicCubicSpline[source]
Bases:
CubicSplineImplementation of the monotonic spline algorithm by Steffen (1990).
- class geoml.math.interpolate.CubicConv2DSeparable(grid)[source]
Bases:
_Interpolator- make_interpolation_matrix(coordinates, derivative=-1)[source]
Generates a sparse matrix for interpolating from a regular grid to a new set of positions in one dimension.
- Parameters:
coordinates (array-like) – Coordinates to interpolate on.
derivative (int) – Direction to derivate on (-1 for no derivative).
- Returns:
interp (InterpolationMatrix) – The interpolator object.