From d335e218fafd933cf5d9ddc952d876f847a6a05e Mon Sep 17 00:00:00 2001 From: Pietro Fanti Date: Thu, 18 Jun 2026 11:19:39 +0200 Subject: [PATCH 1/3] docs(torch): add proper docs for the torch module and make mentions in the README more concise --- README.md | 52 ++++------------------------- docs/api/python.rst | 25 ++++++++++++++ docs/quickstart/examples_python.rst | 37 ++++++++++++++++++++ 3 files changed, 68 insertions(+), 46 deletions(-) diff --git a/README.md b/README.md index 06c5478..ca24646 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,6 @@ - [Input & Output (C++ and Python)](#input--output-c-and-python) - [Minimal Python Example](#minimal-python-example) - [Minimal C++ Example](#minimal-c-example) - - [Minimal PyTorch Example (Differentiable)](#minimal-pytorch-example-differentiable) - [Installation](#installation) - [With conda](#with-conda) - [With pip](#with-pip) @@ -233,51 +232,12 @@ Similarly to Python, the C++ implementation also provides mesh checking capabili > [!TIP] > For reference, have a look at [the main method](./src/main.cpp) of the C++ executable. -### Minimal PyTorch Example (Differentiable) - -Besides the C++/pybind11 interface above, `polyhedral_gravity.torch` provides a pure -PyTorch re-implementation of the same Tsoulis line-integral formula. It is fully -differentiable (autograd) with respect to vertex positions and density, and runs on -both CPU and CUDA. On CPU it is slower than the compiled C++ implementation, but with -a CUDA GPU it can be faster when evaluating many points at once, depending on hardware -and problem size — see the [benchmark notebook](notebooks/polyhedral-gravity-torch-benchmark.ipynb). - -It requires PyTorch, which is **not** a dependency of the base package and must be -installed separately: - -```bash -pip install torch -``` - -```python -import torch -from polyhedral_gravity.torch import evaluate - -cube_vertices = torch.tensor( - [[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], - [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]], - dtype=torch.float64, requires_grad=True, -) -cube_faces = torch.tensor( - [[1, 3, 2], [0, 3, 1], [0, 1, 5], [0, 5, 4], [0, 7, 3], [0, 4, 7], - [1, 2, 6], [1, 6, 5], [2, 3, 6], [3, 7, 6], [4, 5, 6], [4, 6, 7]], -) -cube_density = torch.tensor(1.0, requires_grad=True) -computation_points = torch.tensor([[0.0, 0.0, 2.0]]) - -potential, acceleration, tensor = evaluate( - cube_vertices, cube_faces, cube_density, computation_points, -) - -# Gradients flow back through the analytic formula, e.g. for shape/density optimization -potential.sum().backward() -print(cube_vertices.grad, cube_density.grad) -``` - -> [!TIP] -> See the [PyTorch interface notebook](notebooks/polyhedral-gravity-torch.ipynb) for a -> worked example, and the [benchmark notebook](notebooks/polyhedral-gravity-torch-benchmark.ipynb) -> for a CPU/GPU performance comparison against the C++ implementation. +> [!NOTE] +> An optional, differentiable `polyhedral_gravity.torch` submodule is also available, +> useful if you need autograd gradients (e.g. for shape/density optimization) or want +> to run on a CUDA GPU — see the +> [PyTorch interface documentation](https://esa.github.io/polyhedral-gravity-model/api/python.html#pytorch-interface-differentiable) +> for details and examples. ## Installation diff --git a/docs/api/python.rst b/docs/api/python.rst index b86ec20..1e58c6e 100644 --- a/docs/api/python.rst +++ b/docs/api/python.rst @@ -72,3 +72,28 @@ Below is the list of the available attributes: Specifies the logging level fixed at compile time. This corresponds to the value defined by the ``POLYHEDRAL_GRAVITY_LOGGING_LEVEL`` C++ variable. + + +PyTorch Interface (Differentiable) +----------------------------------- + +The optional :code:`polyhedral_gravity.torch` submodule provides a pure-PyTorch, +autograd-differentiable re-implementation of :code:`evaluate(..)`. It is not +imported by default and requires PyTorch to be installed separately +(:code:`pip install torch`). See :ref:`examples-python` for a usage example, +and the `PyTorch interface notebook `__ +and `benchmark notebook `__ +for further details. + +.. py:function:: polyhedral_gravity.torch.evaluate(vertices, faces, density, computation_points, gravitational_constant=6.67430e-11) + + Gravitational potential, acceleration, and gradient tensor for a polyhedron, + differentiable w.r.t. vertex positions and density. + + :param vertices: ``(N, 3)`` vertex positions [m] (``torch.Tensor``). + :param faces: ``(F, 3)`` triangle vertex indices, integer dtype (``torch.Tensor``). + :param density: Constant density [kg/m^3]. + :param computation_points: ``(Q, 3)`` evaluation positions [m] (``torch.Tensor``). + :param gravitational_constant: Gravitational constant, defaults to ``6.67430e-11``. + :returns: Tuple of ``potential (Q,)``, ``acceleration (Q, 3)``, and gradient ``tensor (Q, 6)`` + (ordered as ``[Vxx, Vyy, Vzz, Vxy, Vxz, Vyz]``). diff --git a/docs/quickstart/examples_python.rst b/docs/quickstart/examples_python.rst index c463f65..f714a74 100644 --- a/docs/quickstart/examples_python.rst +++ b/docs/quickstart/examples_python.rst @@ -227,3 +227,40 @@ Have a look at the example below to see how to use the :code:`GravityEvaluable` # as fast as the following (find the runtime details below), which returns # a list of triplets comprising potential, acceleration, tensor results = evaluable(computation_points, parallel=True) + + +PyTorch Interface (Differentiable) +----------------------------------- + +In addition to the C++-backed :code:`evaluate(..)`, the optional +:code:`polyhedral_gravity.torch` submodule provides a pure-PyTorch +re-implementation of the same Tsoulis (2012) line-integral formula. It is +differentiable with respect to both the polyhedron's vertex positions and its +density, and runs on CPU or CUDA. This is useful for gradient-based +optimization, e.g. recovering a shape or density distribution from gravity +measurements. It requires PyTorch to be installed separately +(:code:`pip install torch`); it is not a hard dependency of +:code:`polyhedral-gravity-model`. + +.. code-block:: python + + import torch + from polyhedral_gravity.torch import evaluate as torch_evaluate + + vertices = torch.tensor(..., dtype=torch.float64, requires_grad=True) # (N, 3) + faces = torch.tensor(...) # (F, 3) int + density = torch.tensor(1.0, requires_grad=True) + computation_points = torch.tensor(...) # (Q, 3) + + potential, acceleration, tensor = torch_evaluate( + vertices, faces, density, computation_points, + ) + + # Gradients flow back through the analytic formula, e.g. for shape/density optimization + potential.sum().backward() + print(vertices.grad, density.grad) + +.. tip:: + See the `PyTorch interface notebook `__ + for a worked example, and the `benchmark notebook `__ + for a CPU/GPU performance comparison against the C++ implementation. From 4ecb54bec3030f89dbe1b253bb63a0726c8ddd8b Mon Sep 17 00:00:00 2001 From: Pietro Fanti Date: Thu, 18 Jun 2026 11:20:38 +0200 Subject: [PATCH 2/3] fix(torch): rename variables to respect PEP8 conventions --- python/polyhedral_gravity/torch.py | 216 ++++++++++++++---------- test/python/test_gravity_model_torch.py | 94 ++++++----- 2 files changed, 177 insertions(+), 133 deletions(-) diff --git a/python/polyhedral_gravity/torch.py b/python/polyhedral_gravity/torch.py index b9f3839..e1cc7ae 100644 --- a/python/polyhedral_gravity/torch.py +++ b/python/polyhedral_gravity/torch.py @@ -11,8 +11,7 @@ from torch import Tensor except ImportError as _err: raise ImportError( - "polyhedral_gravity.torch requires PyTorch. " - "Install it with: pip install torch" + "polyhedral_gravity.torch requires PyTorch. Install it with: pip install torch" ) from _err import math @@ -25,18 +24,18 @@ def evaluate( faces: Tensor, density: float, computation_points: Tensor, - G: float = G_SI, + gravitational_constant: float = G_SI, ) -> tuple[Tensor, Tensor, Tensor]: """ Gravitational potential, acceleration, and gradient tensor for a polyhedron, differentiable w.r.t. vertex positions and density. Args: - vertices: (N, 3) vertex positions [m]. - faces: (F, 3) triangle vertex indices, integer dtype. - density: Constant density [kg/m^3]. - computation_points: (Q, 3) evaluation positions [m]. - G: Gravitational constant, defaults to 6.67430e-11. + vertices: (N, 3) vertex positions [m]. + faces: (F, 3) triangle vertex indices, integer dtype. + density: Constant density [kg/m^3]. + computation_points: (Q, 3) evaluation positions [m]. + gravitational_constant: Gravitational constant, defaults to 6.67430e-11. Returns: potential: (Q,) gravitational potential [m^2/s^2]. @@ -45,46 +44,50 @@ def evaluate( components ordered as [Vxx, Vyy, Vzz, Vxy, Vxz, Vyz]. """ - EPS = 1e-30 + eps = 1e-30 # face vertex coords: (F, 3, 3) - v = vertices[faces] # [f, vertex_idx, xyz] + face_vertices = vertices[faces] # [f, vertex_idx, xyz] # shift so each query point is at the origin: (Q, F, 3, 3) - v_p = v[None] - computation_points[:, None, None, :] + face_vertices_rel = face_vertices[None] - computation_points[:, None, None, :] - v0 = v_p[:, :, 0] # (Q, F, 3) - v1 = v_p[:, :, 1] - v2 = v_p[:, :, 2] + v0 = face_vertices_rel[:, :, 0] # (Q, F, 3) + v1 = face_vertices_rel[:, :, 1] + v2 = face_vertices_rel[:, :, 2] # ------------------------------------------------------------------ # - # Edge vectors G_pq and outward unit normal N_p + # Edge vectors G_pq and outward unit normal (face_normal) # ------------------------------------------------------------------ # - G0 = v1 - v0 # (Q, F, 3) - G1 = v2 - v1 - G2 = v0 - v2 + edge0 = v1 - v0 # (Q, F, 3) + edge1 = v2 - v1 + edge2 = v0 - v2 - G_segs = torch.stack([G0, G1, G2], dim=2) # (Q, F, 3, 3) [q,f,seg,xyz] - V_segs = torch.stack([v0, v1, v2], dim=2) # (Q, F, 3, 3) segment start vertices + edges = torch.stack([edge0, edge1, edge2], dim=2) # (Q, F, 3, 3) [q,f,seg,xyz] + segment_start_points = torch.stack( + [v0, v1, v2], dim=2 + ) # (Q, F, 3, 3) segment start vertices - cross_01 = torch.linalg.cross(G0, G1) # (Q, F, 3) - N_p = cross_01 / cross_01.norm(dim=-1, keepdim=True).clamp(min=EPS) + cross_01 = torch.linalg.cross(edge0, edge1) # (Q, F, 3) + face_normal = cross_01 / cross_01.norm(dim=-1, keepdim=True).clamp(min=eps) # ------------------------------------------------------------------ # - # Segment unit normals n_pq = (G_pq x N_p) / |...| + # Segment unit normals n_pq = (G_pq x face_normal) / |...| # ------------------------------------------------------------------ # - N_p_e = N_p[:, :, None].expand_as(G_segs) # (Q, F, 3, 3) - n_pq = torch.linalg.cross(G_segs, N_p_e) - n_pq = n_pq / n_pq.norm(dim=-1, keepdim=True).clamp(min=EPS) + face_normal_expanded = face_normal[:, :, None].expand_as(edges) # (Q, F, 3, 3) + n_pq = torch.linalg.cross(edges, face_normal_expanded) + n_pq = n_pq / n_pq.norm(dim=-1, keepdim=True).clamp(min=eps) # ------------------------------------------------------------------ # # Signed plane distance, plane normal orientation, projection P' # ------------------------------------------------------------------ # - hp_signed = (N_p * v0).sum(dim=-1) # (Q, F) = N_p . v0 + hp_signed = (face_normal * v0).sum(dim=-1) # (Q, F) = N_p . v0 h_p = hp_signed.abs() - sigma_p = hp_signed.sign() # sigma_p + sigma_p = hp_signed.sign() # sigma_p - Pp = hp_signed[..., None] * N_p # (Q, F, 3) P' = hp_signed * N_p + plane_projection = ( + hp_signed[..., None] * face_normal + ) # (Q, F, 3) P' = hp_signed * face_normal # ------------------------------------------------------------------ # # Segment normal orientations sigma_pq and edge distances h_pq @@ -95,116 +98,155 @@ def evaluate( # gradient graph; the equivalent norm-based formula (||P''-P'||) shares # the same value but its graph misses that contribution. # ------------------------------------------------------------------ # - Pp_e = Pp[:, :, None].expand_as(V_segs) # (Q, F, 3, 3) - h_pq_raw = -(n_pq * (Pp_e - V_segs)).sum(dim=-1) # (Q, F, 3) + plane_projection_expanded = plane_projection[:, :, None].expand_as( + segment_start_points + ) # (Q, F, 3, 3) + h_pq_raw = -(n_pq * (plane_projection_expanded - segment_start_points)).sum( + dim=-1 + ) # (Q, F, 3) # Threshold relative to edge length: floating-point noise in # normalised-vector arithmetic can produce |h_pq_raw| / edge_scale ~ 1e-16 # for algebraically-zero distances. Zeroing these before .sign() prevents # misclassifying the singularity case (all_inside vs on_edge/on_vertex), - # which would corrupt singA by 2-4x. - G_norm = G_segs.norm(dim=-1) # (Q, F, 3) - edge_scale = G_norm.amax(dim=-1, keepdim=True) # (Q, F, 1) + # which would corrupt singularity_scalar by 2-4x. + edge_length = edges.norm(dim=-1) # (Q, F, 3) + edge_scale = edge_length.amax(dim=-1, keepdim=True) # (Q, F, 1) h_pq_signed = torch.where( h_pq_raw.abs() < 1e-10 * edge_scale, torch.zeros_like(h_pq_raw), h_pq_raw, ) - sigma_pq = h_pq_signed.sign() # (Q, F, 3) - h_pq = h_pq_signed.abs() # (Q, F, 3) + sigma_pq = h_pq_signed.sign() # (Q, F, 3) + h_pq = h_pq_signed.abs() # (Q, F, 3) # ------------------------------------------------------------------ # # Projection parameter t along each edge (needed for s1, s2) # ------------------------------------------------------------------ # - G_sq = (G_segs * G_segs).sum(dim=-1).clamp(min=EPS) # (Q, F, 3) - t_pq = ((Pp_e - V_segs) * G_segs).sum(dim=-1) / G_sq # (Q, F, 3) + edge_length_sq = (edges * edges).sum(dim=-1).clamp(min=eps) # (Q, F, 3) + t_pq = ((plane_projection_expanded - segment_start_points) * edges).sum( + dim=-1 + ) / edge_length_sq # (Q, F, 3) # ------------------------------------------------------------------ # # Distances from P (origin) to segment endpoints # ------------------------------------------------------------------ # - V_segs_end = torch.roll(V_segs, -1, dims=2) # v_end = v_{(q+1)%3} - l1 = V_segs.norm(dim=-1) # (Q, F, 3) |v_start| - l2 = V_segs_end.norm(dim=-1) # (Q, F, 3) |v_end| + segment_end_points = torch.roll( + segment_start_points, -1, dims=2 + ) # v_end = v_{(q+1)%3} + l1 = segment_start_points.norm(dim=-1) # (Q, F, 3) |v_start| + l2 = segment_end_points.norm(dim=-1) # (Q, F, 3) |v_end| # signed 1-D distances along the edge: s1 = -t*|G|, s2 = (1-t)*|G| - s1 = -t_pq * G_norm - s2 = (1.0 - t_pq) * G_norm + s1 = -t_pq * edge_length + s2 = (1.0 - t_pq) * edge_length # ------------------------------------------------------------------ # # Transcendental expressions LN_pq and AN_pq (Tsoulis Eqs 14-15) # ------------------------------------------------------------------ # - ln_num = (s2 + l2).clamp(min=EPS) - ln_den = (s1 + l1).clamp(min=EPS) - LN = torch.log(ln_num / ln_den) + ln_num = (s2 + l2).clamp(min=eps) + ln_den = (s1 + l1).clamp(min=eps) + log_term = torch.log(ln_num / ln_den) - h_p_e = h_p[:, :, None].expand_as(h_pq) # (Q, F, 3) + h_p_e = h_p[:, :, None].expand_as(h_pq) # (Q, F, 3) degenerate = (h_p_e < 1e-15) | (h_pq < 1e-15) - AN = torch.atan(h_p_e * s2 / (h_pq * l2).clamp(min=EPS)) \ - - torch.atan(h_p_e * s1 / (h_pq * l1).clamp(min=EPS)) - AN = torch.where(degenerate, torch.zeros_like(AN), AN) + arctan_term = torch.atan(h_p_e * s2 / (h_pq * l2).clamp(min=eps)) - torch.atan( + h_p_e * s1 / (h_pq * l1).clamp(min=eps) + ) + arctan_term = torch.where(degenerate, torch.zeros_like(arctan_term), arctan_term) # ------------------------------------------------------------------ # # Per-face sums S1, S2, singularity correction (Tsoulis Eqs 11-12) # ------------------------------------------------------------------ # - sum1 = (h_pq_signed * LN).sum(dim=-1) # (Q, F) - sum2 = (sigma_pq * AN).sum(dim=-1) # (Q, F) - sum1_tensor = (n_pq * LN[..., None]).sum(dim=-2) # (Q, F, 3) for gradient tensor + sum1 = (h_pq_signed * log_term).sum(dim=-1) # (Q, F) + sum2 = (sigma_pq * arctan_term).sum(dim=-1) # (Q, F) + sum1_tensor = (n_pq * log_term[..., None]).sum( + dim=-2 + ) # (Q, F, 3) for gradient tensor - # singA correction - three cases detected from the sigma_pq pattern: + # Singularity correction - three cases detected from the sigma_pq pattern: # all +1 -> P' strictly inside face -> -2*pi*h_p # two +1, one 0 -> P' on a face edge -> -pi*h_p # one +1, two 0 -> P' at a face vertex -> -theta_v*h_p - n_pos = (sigma_pq > 0).sum(dim=-1) # (Q, F) - n_zero = (sigma_pq == 0).sum(dim=-1) # (Q, F) + n_pos = (sigma_pq > 0).sum(dim=-1) # (Q, F) + n_zero = (sigma_pq == 0).sum(dim=-1) # (Q, F) all_inside = n_pos == 3 - on_edge = (n_pos == 2) & (n_zero == 1) - on_vertex = (n_pos == 1) & (n_zero == 2) - - G_hat = G_segs / G_norm[..., None].clamp(min=EPS) # (Q, F, 3, 3) - cos_t0 = -(G_hat[..., 0, :] * G_hat[..., 2, :]).sum(dim=-1) # (Q, F) - cos_t1 = -(G_hat[..., 1, :] * G_hat[..., 0, :]).sum(dim=-1) - cos_t2 = -(G_hat[..., 2, :] * G_hat[..., 1, :]).sum(dim=-1) + on_edge = (n_pos == 2) & (n_zero == 1) + on_vertex = (n_pos == 1) & (n_zero == 2) + + edge_unit_vector = edges / edge_length[..., None].clamp(min=eps) # (Q, F, 3, 3) + cos_t0 = -(edge_unit_vector[..., 0, :] * edge_unit_vector[..., 2, :]).sum( + dim=-1 + ) # (Q, F) + cos_t1 = -(edge_unit_vector[..., 1, :] * edge_unit_vector[..., 0, :]).sum(dim=-1) + cos_t2 = -(edge_unit_vector[..., 2, :] * edge_unit_vector[..., 1, :]).sum(dim=-1) theta_0 = torch.acos(cos_t0.clamp(-1.0, 1.0)) theta_1 = torch.acos(cos_t1.clamp(-1.0, 1.0)) theta_2 = torch.acos(cos_t2.clamp(-1.0, 1.0)) - theta_v = (theta_2 * (sigma_pq[..., 0] > 0).to(h_p.dtype) + - theta_0 * (sigma_pq[..., 1] > 0).to(h_p.dtype) + - theta_1 * (sigma_pq[..., 2] > 0).to(h_p.dtype)) # (Q, F) - - sing_A = torch.where(all_inside, -2.0 * math.pi * h_p, - torch.where(on_edge, -math.pi * h_p, - torch.where(on_vertex, -theta_v * h_p, - torch.zeros_like(h_p)))) + theta_v = ( + theta_2 * (sigma_pq[..., 0] > 0).to(h_p.dtype) + + theta_0 * (sigma_pq[..., 1] > 0).to(h_p.dtype) + + theta_1 * (sigma_pq[..., 2] > 0).to(h_p.dtype) + ) # (Q, F) + + singularity_scalar = torch.where( + all_inside, + -2.0 * math.pi * h_p, + torch.where( + on_edge, + -math.pi * h_p, + torch.where(on_vertex, -theta_v * h_p, torch.zeros_like(h_p)), + ), + ) - sing_B = torch.where(all_inside[..., None], -2.0 * math.pi * sigma_p[..., None] * N_p, - torch.where(on_edge[..., None], -math.pi * sigma_p[..., None] * N_p, - torch.where(on_vertex[..., None], -theta_v[..., None] * sigma_p[..., None] * N_p, - torch.zeros_like(N_p)))) # (Q, F, 3) + singularity_vector = torch.where( + all_inside[..., None], + -2.0 * math.pi * sigma_p[..., None] * face_normal, + torch.where( + on_edge[..., None], + -math.pi * sigma_p[..., None] * face_normal, + torch.where( + on_vertex[..., None], + -theta_v[..., None] * sigma_p[..., None] * face_normal, + torch.zeros_like(face_normal), + ), + ), + ) # (Q, F, 3) - A_p = sum1 + h_p * sum2 + sing_A # (Q, F) per-face scalar sum + face_sum = sum1 + h_p * sum2 + singularity_scalar # (Q, F) per-face scalar sum # ------------------------------------------------------------------ # # Sum over faces -> potential and acceleration # ------------------------------------------------------------------ # - V_contributions = sigma_p * h_p * A_p # (Q, F) - g_contributions = N_p * A_p[..., None] # (Q, F, 3) + potential_terms = sigma_p * h_p * face_sum # (Q, F) + g_contributions = face_normal * face_sum[..., None] # (Q, F, 3) - potential = (G * density / 2.0) * V_contributions.sum(dim=-1) - acceleration = -G * density * g_contributions.sum(dim=-2) + potential = (gravitational_constant * density / 2.0) * potential_terms.sum(dim=-1) + acceleration = -gravitational_constant * density * g_contributions.sum(dim=-2) # ------------------------------------------------------------------ # # Gradient tensor (Tsoulis Eq. 13) # ------------------------------------------------------------------ # - subSum = sum1_tensor + N_p * (sigma_p * sum2)[..., None] + sing_B # (Q, F, 3) + tensor_subterm = ( + sum1_tensor + face_normal * (sigma_p * sum2)[..., None] + singularity_vector + ) # (Q, F, 3) - # diagonal components: Vxx, Vyy, Vzz -> N_p[i] * subSum[i] - diag = N_p * subSum # (Q, F, 3) + # diagonal components: Vxx, Vyy, Vzz -> face_normal[i] * tensor_subterm[i] + diag = face_normal * tensor_subterm # (Q, F, 3) # off-diagonal: Vxy=N[0]*sub[1], Vxz=N[0]*sub[2], Vyz=N[1]*sub[2] - off_N = torch.stack([N_p[..., 0], N_p[..., 0], N_p[..., 1]], dim=-1) - off_sub = torch.stack([subSum[..., 1], subSum[..., 2], subSum[..., 2]], dim=-1) - offdiag = off_N * off_sub # (Q, F, 3) + offdiag_normal = torch.stack( + [face_normal[..., 0], face_normal[..., 0], face_normal[..., 1]], dim=-1 + ) + offdiag_subterm = torch.stack( + [tensor_subterm[..., 1], tensor_subterm[..., 2], tensor_subterm[..., 2]], dim=-1 + ) + offdiag = offdiag_normal * offdiag_subterm # (Q, F, 3) - tensor = G * density * torch.cat([diag, offdiag], dim=-1).sum(dim=-2) # (Q, 6) + tensor = ( + gravitational_constant + * density + * torch.cat([diag, offdiag], dim=-1).sum(dim=-2) + ) # (Q, 6) return potential, acceleration, tensor diff --git a/test/python/test_gravity_model_torch.py b/test/python/test_gravity_model_torch.py index e5c4617..d3981b1 100644 --- a/test/python/test_gravity_model_torch.py +++ b/test/python/test_gravity_model_torch.py @@ -110,13 +110,13 @@ def _rel_err_acc(mine, ref): return np.linalg.norm(mine - ref, axis=-1) / (np.linalg.norm(ref, axis=-1) + 1e-30) -def _esa_batch(poly, pts): - Vs, gs = [], [] +def _cpp_batch(poly, pts): + potentials, gs = [], [] for p in pts: - V, g, _ = evaluate(poly, p) - Vs.append(V) + potential, g, _ = evaluate(poly, p) + potentials.append(potential) gs.append(g) - return np.array(Vs), np.array(gs) + return np.array(potentials), np.array(gs) # --------------------------------------------------------------------------- @@ -124,7 +124,7 @@ def _esa_batch(poly, pts): # --------------------------------------------------------------------------- def _make_body(b): - esa_poly = Polyhedron( + cpp_poly = Polyhedron( (b["vertices"], b["faces"]), density=1.0, normal_orientation=NormalOrientation.OUTWARDS, @@ -133,7 +133,7 @@ def _make_body(b): **b, "verts_t": torch.tensor(b["vertices"], dtype=torch.float64), "faces_t": torch.tensor(b["faces"], dtype=torch.long), - "esa_poly": esa_poly, + "cpp_poly": cpp_poly, } @@ -162,18 +162,18 @@ class TestAccuracy: def test_exterior_points(self, body): pts = body["exterior_points"] - V, g, _ = torch_evaluate( + potential, g, _ = torch_evaluate( body["verts_t"], body["faces_t"], 1.0, torch.tensor(pts, dtype=torch.float64), ) - V_np = V.detach().numpy() + potential_np = potential.detach().numpy() g_np = g.detach().numpy() for i, p in enumerate(pts): - esa_V, esa_g, _ = evaluate(body["esa_poly"], p) - assert _rel_err_pot(V_np[i], esa_V) < REL_TOL, \ + cpp_potential, cpp_g, _ = evaluate(body["cpp_poly"], p) + assert _rel_err_pot(potential_np[i], cpp_potential) < REL_TOL, \ f"body={body['id']} point={p}: potential rel error too large" - assert _rel_err_acc(g_np[i], np.array(esa_g)) < REL_TOL, \ + assert _rel_err_acc(g_np[i], np.array(cpp_g)) < REL_TOL, \ f"body={body['id']} point={p}: acceleration rel error too large" @@ -182,25 +182,25 @@ class TestSingularityCases: def test_edge_singularity(self, body_with_singularities): body = body_with_singularities - V, g, _ = torch_evaluate( + potential, g, _ = torch_evaluate( body["verts_t"], body["faces_t"], 1.0, torch.tensor(body["edge_singular_points"], dtype=torch.float64), ) for i, p in enumerate(body["edge_singular_points"]): - esa_V, esa_g, _ = evaluate(body["esa_poly"], p) - assert _rel_err_pot(V.detach().numpy()[i], esa_V) < REL_TOL - assert _rel_err_acc(g.detach().numpy()[i], np.array(esa_g)) < REL_TOL + cpp_potential, cpp_g, _ = evaluate(body["cpp_poly"], p) + assert _rel_err_pot(potential.detach().numpy()[i], cpp_potential) < REL_TOL + assert _rel_err_acc(g.detach().numpy()[i], np.array(cpp_g)) < REL_TOL def test_vertex_singularity(self, body_with_singularities): body = body_with_singularities - V, g, _ = torch_evaluate( + potential, g, _ = torch_evaluate( body["verts_t"], body["faces_t"], 1.0, torch.tensor(body["vertex_singular_points"], dtype=torch.float64), ) for i, p in enumerate(body["vertex_singular_points"]): - esa_V, esa_g, _ = evaluate(body["esa_poly"], p) - assert _rel_err_pot(V.detach().numpy()[i], esa_V) < REL_TOL - assert _rel_err_acc(g.detach().numpy()[i], np.array(esa_g)) < REL_TOL + cpp_potential, cpp_g, _ = evaluate(body["cpp_poly"], p) + assert _rel_err_pot(potential.detach().numpy()[i], cpp_potential) < REL_TOL + assert _rel_err_acc(g.detach().numpy()[i], np.array(cpp_g)) < REL_TOL class TestBatched: @@ -218,12 +218,12 @@ def sphere_points(self): ], dtype=np.float64) def test_potential_batch(self, body, sphere_points): - V, _, __ = torch_evaluate( + potential, _, __ = torch_evaluate( body["verts_t"], body["faces_t"], 1.0, torch.tensor(sphere_points, dtype=torch.float64), ) - ref_V, _ = _esa_batch(body["esa_poly"], sphere_points) - err = _rel_err_pot(V.detach().numpy(), ref_V) + ref_potential, _ = _cpp_batch(body["cpp_poly"], sphere_points) + err = _rel_err_pot(potential.detach().numpy(), ref_potential) assert err.max() < REL_TOL, \ f"body={body['id']}: max potential rel error = {err.max():.2e}" @@ -232,7 +232,7 @@ def test_acceleration_batch(self, body, sphere_points): body["verts_t"], body["faces_t"], 1.0, torch.tensor(sphere_points, dtype=torch.float64), ) - _, ref_g = _esa_batch(body["esa_poly"], sphere_points) + _, ref_g = _cpp_batch(body["cpp_poly"], sphere_points) err = _rel_err_acc(g.detach().numpy(), ref_g) assert err.max() < REL_TOL, \ f"body={body['id']}: max acceleration rel error = {err.max():.2e}" @@ -244,23 +244,23 @@ class TestGradients: def test_gradient_wrt_vertices_nonzero(self, body): verts = body["verts_t"].clone().requires_grad_(True) q = torch.tensor(body["exterior_points"][:1], dtype=torch.float64) - V, _, __ = torch_evaluate(verts, body["faces_t"], 1.0, q) - V.sum().backward() + potential, _, __ = torch_evaluate(verts, body["faces_t"], 1.0, q) + potential.sum().backward() assert verts.grad is not None and verts.grad.abs().max() > 0 def test_gradient_wrt_scalar_density_nonzero(self, body): density = torch.tensor(1.0, dtype=torch.float64, requires_grad=True) q = torch.tensor(body["exterior_points"][:1], dtype=torch.float64) - V, _, __ = torch_evaluate(body["verts_t"], body["faces_t"], density, q) - V.sum().backward() + potential, _, __ = torch_evaluate(body["verts_t"], body["faces_t"], density, q) + potential.sum().backward() assert density.grad is not None and density.grad.abs() > 0 def test_gradient_consistent_with_finite_differences(self, body): """∂V/∂vertices must agree with central finite differences to 1e-6.""" verts = body["verts_t"].clone().requires_grad_(True) q = torch.tensor(body["grad_point"], dtype=torch.float64) - V, _, __ = torch_evaluate(verts, body["faces_t"], 1.0, q) - V.backward() + potential, _, __ = torch_evaluate(verts, body["faces_t"], 1.0, q) + potential.backward() grad = verts.grad.clone() eps = 1e-5 @@ -268,11 +268,13 @@ def test_gradient_consistent_with_finite_differences(self, body): base = body["verts_t"] for i in range(base.shape[0]): for j in range(3): - vp = base.clone(); vp[i, j] += eps - vm = base.clone(); vm[i, j] -= eps - Vp, _, __ = torch_evaluate(vp, body["faces_t"], 1.0, q) - Vm, _, __ = torch_evaluate(vm, body["faces_t"], 1.0, q) - fd_grad[i, j] = (Vp - Vm) / (2 * eps) + vp = base.clone() + vp[i, j] += eps + vm = base.clone() + vm[i, j] -= eps + potential_plus, _, __ = torch_evaluate(vp, body["faces_t"], 1.0, q) + potential_minus, _, __ = torch_evaluate(vm, body["faces_t"], 1.0, q) + fd_grad[i, j] = (potential_plus - potential_minus) / (2 * eps) rel = (grad - fd_grad).abs() / (fd_grad.abs() + 1e-30) assert rel.max().item() < 1e-6, \ @@ -284,27 +286,27 @@ class TestTensor: def test_tensor_accuracy(self, body): pts = body["exterior_points"] - _, _, T = torch_evaluate( + _, _, tensor = torch_evaluate( body["verts_t"], body["faces_t"], 1.0, torch.tensor(pts, dtype=torch.float64), ) - T_np = T.detach().numpy() + tensor_np = tensor.detach().numpy() for i, p in enumerate(pts): - _, _, esa_T = evaluate(body["esa_poly"], p) - ref = np.array(esa_T) + _, _, cpp_tensor = evaluate(body["cpp_poly"], p) + ref = np.array(cpp_tensor) # Normalise by the problem scale (max component) so near-zero diagonal # components at high-symmetry points don't blow up the relative error. scale = np.abs(ref).max() + 1e-30 - err = np.abs(T_np[i] - ref) / scale + err = np.abs(tensor_np[i] - ref) / scale assert err.max() < REL_TOL, \ f"body={body['id']} point={p}: tensor rel error too large: {err.max():.2e}" def test_tensor_scaling(self, body): q = torch.tensor(body["exterior_points"][:2], dtype=torch.float64) - _, _, T1 = torch_evaluate(body["verts_t"], body["faces_t"], 1.0, q) - _, _, Tr = torch_evaluate(body["verts_t"], body["faces_t"], 2.0, q) - assert torch.allclose(Tr, 2.0 * T1, rtol=1e-12) + _, _, tensor_unit = torch_evaluate(body["verts_t"], body["faces_t"], 1.0, q) + _, _, tensor_scaled = torch_evaluate(body["verts_t"], body["faces_t"], 2.0, q) + assert torch.allclose(tensor_scaled, 2.0 * tensor_unit, rtol=1e-12) class TestDensityScaling: @@ -313,8 +315,8 @@ class TestDensityScaling: @pytest.mark.parametrize("rho", [0.5, 2.0, 1000.0]) def test_linear_scaling(self, body, rho): q = torch.tensor(body["exterior_points"][:2], dtype=torch.float64) - V1, g1, _ = torch_evaluate(body["verts_t"], body["faces_t"], 1.0, q) - Vr, gr, _ = torch_evaluate(body["verts_t"], body["faces_t"], rho, q) - assert torch.allclose(Vr, rho * V1, rtol=1e-12) + potential_unit, g1, _ = torch_evaluate(body["verts_t"], body["faces_t"], 1.0, q) + potential_scaled, gr, _ = torch_evaluate(body["verts_t"], body["faces_t"], rho, q) + assert torch.allclose(potential_scaled, rho * potential_unit, rtol=1e-12) assert torch.allclose(gr, rho * g1, rtol=1e-12) From ab7afc611a8fe292fd6e85e0ac03fb34ebcb17cb Mon Sep 17 00:00:00 2001 From: Pietro Fanti Date: Thu, 18 Jun 2026 11:38:02 +0200 Subject: [PATCH 3/3] fix: make magic attributes accessible again --- python/polyhedral_gravity/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/polyhedral_gravity/__init__.py b/python/polyhedral_gravity/__init__.py index 6eb8646..958b2cc 100644 --- a/python/polyhedral_gravity/__init__.py +++ b/python/polyhedral_gravity/__init__.py @@ -1 +1,2 @@ from ._core import * +from ._core import __version__, __parallelization__, __commit__, __logging__