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.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 the origin, 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 origin happens to be. 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.
- 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.
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.