@@ -92,6 +92,10 @@ tri.circumscribed_circle(simplex)      # → (center, radius)
 tri.volume(simplex)                    # Simplex volume
 tri.volumes()                          # All simplex volumes
 tri.point_in_simplex(point, simplex)   # Containment test
+tri.simplices_containing(point)        # All simplices containing a point, in one call
+                                       # instead of a point_in_simplex loop; pass a known
+                                       # containing simplex via simplex=... to skip the
+                                       # locate step, or restrict with candidates=...
 tri.point_in_circumcircle(pt, simplex) # Circumcircle test
 tri.bowyer_watson(pt_index)            # Direct Bowyer-Watson
 tri.get_opposing_vertices(simplex)     # Facet neighbours' opposite vertices
@@ -111,6 +115,8 @@ from adaptive_triangulation import (
     point_in_simplex,          # Containment test
     volume,                    # Simplex volume
     simplex_volume_in_embedding,  # Volume in embedding space
+    default_loss,              # LearnerND's default loss (embedded simplex volume),
+                               # signature-compatible with loss_per_simplex
     orientation,               # Face orientation
 )
 ```
 
@@ -0,0 +1,118 @@
+"""Measure the LearnerND speedup from `simplices_containing` and `default_loss`.
+
+This monkeypatches adaptive (the PR #493 branch, with the Rust backend
+active) to use the two batched APIs, the way the adaptive-side integration
+would:
+
+- ``tell_pending`` distributes a pending point over the simplices containing
+  it with one ``tri.simplices_containing(point, simplex=hint)`` call instead
+  of a Python loop of per-neighbour ``tri.point_in_simplex`` checks.
+- The default loss is the Rust ``default_loss`` instead of the Python
+  wrapper around ``simplex_volume_in_embedding``.
+
+Requirements:
+    pip install adaptive-triangulation
+    pip install "adaptive @ git+https://github.com/python-adaptive/adaptive@refs/pull/493/head"
+"""
+
+from __future__ import annotations
+
+import time
+from typing import TYPE_CHECKING
+
+import adaptive_triangulation as at
+import numpy as np
+from adaptive.learner import learnerND as lnd_mod
+from adaptive.learner.learnerND import LearnerND
+
+if TYPE_CHECKING:
+    from collections.abc import Callable, Sequence
+
+reference_default_loss = lnd_mod.default_loss
+
+
+def ring_of_fire(xy: Sequence[float]) -> float:
+    x, y = xy
+    a, d = 0.2, 0.5
+    return x + np.exp(-((x**2 + y**2 - d**2) ** 2) / a**4)
+
+
+def sphere_of_fire(xyz: Sequence[float]) -> float:
+    x, y, z = xyz
+    a, d = 0.2, 0.5
+    return x + np.exp(-((x**2 + y**2 + z**2 - d**2) ** 2) / a**4)
+
+
+def _add_pending_to_simplex(self: LearnerND, point: tuple, simplex: tuple) -> tuple:
+    """`_try_adding_pending_point_to_simplex` without the containment check
+    (the caller already knows `simplex` contains `point`)."""
+    if simplex not in self._subtriangulations:
+        vertices = self.tri.get_vertices(simplex)
+        self._subtriangulations[simplex] = self._triangulation_class(vertices)
+    self._pending_to_simplex[point] = simplex
+    return self._subtriangulations[simplex].add_point(point)
+
+
+def tell_pending_batched(self: LearnerND, point: tuple, *, simplex: tuple | None = None) -> None:
+    """tell_pending with the point_in_simplex loop replaced by one
+    `simplices_containing` call (the `simplex` argument becomes the hint)."""
+    point = tuple(point)
+    if not self.inside_bounds(point):
+        return
+    self.pending_points.add(point)
+    if self.tri is None:
+        return
+    containing = self.tri.simplices_containing(point, simplex=simplex)
+    for simpl in containing:
+        _, to_add = _add_pending_to_simplex(self, point, simpl)
+        if to_add is None:
+            continue
+        self._update_subsimplex_losses(simpl, to_add)
+
+
+def run(
+    func: Callable[[Sequence[float]], float],
+    bounds: list[tuple[float, float]],
+    n_points: int,
+    *,
+    rust_loss: bool,
+    batched_tell: bool,
+) -> tuple[float, LearnerND]:
+    lnd_mod.default_loss = at.default_loss if rust_loss else reference_default_loss
+    learner = LearnerND(func, bounds=bounds)
+    if batched_tell:
+        learner.tell_pending = tell_pending_batched.__get__(learner)
+    t0 = time.perf_counter()
+    for _ in range(n_points):
+        points, _ = learner.ask(1)
+        for p in points:
+            learner.tell(p, func(p))
+    elapsed = time.perf_counter() - t0
+    return elapsed, learner
+
+
+CONFIGS = [
+    ("baseline (Rust backend, PR #493)", False, False),
+    ("+ rust default_loss", True, False),
+    ("+ simplices_containing tell_pending", False, True),
+    ("+ both", True, True),
+]
+
+for name, func, bounds, n in [
+    ("2D ring_of_fire", ring_of_fire, [(-1, 1), (-1, 1)], 3000),
+    ("3D sphere_of_fire", sphere_of_fire, [(-1, 1), (-1, 1), (-1, 1)], 1500),
+]:
+    print(f"\n{name}, {n} points:")
+    baseline = None
+    base_points = None
+    for label, rust_loss, batched_tell in CONFIGS:
+        best = min(
+            run(func, bounds, n, rust_loss=rust_loss, batched_tell=batched_tell)[0]
+            for _ in range(3)
+        )
+        _, learner = run(func, bounds, n, rust_loss=rust_loss, batched_tell=batched_tell)
+        pts = sorted(learner.data)
+        if baseline is None:
+            baseline, base_points = best, pts
+        same = "identical points" if pts == base_points else "DIFFERENT POINTS"
+        print(f"  {label:40s} {best:6.2f}s  ({baseline / best:4.2f}x, {same})")
@@ -29,6 +29,7 @@
     VerticesProxy,
     __version__,
     circumsphere,
+    default_loss,
     fast_2d_circumcircle,
     fast_2d_point_in_simplex,
     fast_3d_circumcircle,
@@ -47,6 +48,7 @@
     "VerticesProxy",
     "__version__",
     "circumsphere",
+    "default_loss",
     "fast_2d_circumcircle",
     "fast_2d_point_in_simplex",
     "fast_3d_circumcircle",
 
@@ -16,6 +16,9 @@ from ._rust import (
 from ._rust import (
     circumsphere as circumsphere,
 )
+from ._rust import (
+    default_loss as default_loss,
+)
 from ._rust import (
     fast_2d_circumcircle as fast_2d_circumcircle,
 )
@@ -51,6 +54,7 @@ __all__: list[str] = [
     "VerticesProxy",
     "__version__",
     "circumsphere",
+    "default_loss",
     "fast_2d_circumcircle",
     "fast_2d_point_in_simplex",
     "fast_3d_circumcircle",
 
@@ -94,6 +94,13 @@ class Triangulation:
         simplex: SimplexLike,
         eps: float = 1e-8,
     ) -> bool: ...
+    def simplices_containing(
+        self,
+        point: PointLike,
+        simplex: SimplexLike | None = None,
+        candidates: Iterable[SimplexLike] | None = None,
+        eps: float = 1e-8,
+    ) -> list[Simplex]: ...
     def bowyer_watson(
         self,
         pt_index: int,
@@ -104,6 +111,11 @@ class Triangulation:
     def convex_invariant(self, vertex: int) -> bool: ...
 
 def circumsphere(points: PointsLike) -> tuple[tuple[float, ...], float]: ...
+def default_loss(
+    simplex: PointsLike,
+    values: Sequence[float] | PointsLike,
+    value_scale: float | None = None,
+) -> float: ...
 def fast_2d_circumcircle(points: PointsLike) -> tuple[tuple[float, ...], float]: ...
 def fast_2d_point_in_simplex(point: PointLike, simplex: PointsLike, eps: float = 1e-8) -> bool: ...
 def fast_3d_circumcircle(points: PointsLike) -> tuple[tuple[float, ...], float]: ...
 
@@ -666,6 +666,19 @@ pub fn simplex_volume_in_embedding(vertices: &[Vec<f64>]) -> Result<f64, Geometr
     Ok(vol_square.sqrt())
 }
 
+/// The loss adaptive's `LearnerND` assigns to a simplex by default: the
+/// volume of the simplex embedded in (input + output)-dimensional space,
+/// where each vertex is extended with its function value(s). `simplex` and
+/// `values` are zipped pairwise, like the reference implementation.
+pub fn default_loss(simplex: &[Vec<f64>], values: &[Vec<f64>]) -> Result<f64, GeometryError> {
+    let embedded: Vec<Vec<f64>> = simplex
+        .iter()
+        .zip(values)
+        .map(|(vertex, value)| vertex.iter().chain(value.iter()).copied().collect())
+        .collect();
+    simplex_volume_in_embedding(&embedded)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
 
@@ -8,6 +8,7 @@ pub mod py;
 pub mod tolerances;
 pub mod triangulation;
 
+use numpy::{PyReadonlyArray1, PyReadonlyArray2};
 use pyo3::exceptions::{PyValueError, PyZeroDivisionError};
 use pyo3::prelude::*;
 use pyo3::types::{PyAny, PyModule, PyTuple};
@@ -47,6 +48,7 @@ fn _rust(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(py_fast_norm, m)?)?;
     m.add_function(wrap_pyfunction!(py_volume, m)?)?;
     m.add_function(wrap_pyfunction!(py_simplex_volume_in_embedding, m)?)?;
+    m.add_function(wrap_pyfunction!(py_default_loss, m)?)?;
     m.add_function(wrap_pyfunction!(py_orientation, m)?)?;
 
     Ok(())
@@ -176,6 +178,54 @@ fn py_simplex_volume_in_embedding(vertices: &Bound<'_, PyAny>) -> PyResult<f64>
         .map_err(|err| PyValueError::new_err(err.to_string()))
 }
 
+/// Parse per-vertex function values: one scalar per vertex or one vector per
+/// vertex, with bulk-copy fast paths for 1-D and 2-D f64 numpy arrays (what
+/// `LearnerND._compute_loss` passes after scaling).
+fn parse_values(obj: &Bound<'_, PyAny>) -> PyResult<Vec<Vec<f64>>> {
+    if let Ok(array) = obj.extract::<PyReadonlyArray1<'_, f64>>() {
+        return Ok(array.as_array().iter().map(|&value| vec![value]).collect());
+    }
+    if let Ok(array) = obj.extract::<PyReadonlyArray2<'_, f64>>() {
+        return Ok(array
+            .as_array()
+            .outer_iter()
+            .map(|row| row.to_vec())
+            .collect());
+    }
+    let Ok(iter) = obj.try_iter() else {
+        return Err(PyValueError::new_err(
+            "Expected one scalar or one vector of floats per vertex",
+        ));
+    };
+    let mut values = Vec::new();
+    for item in iter {
+        let item = item?;
+        match item.extract::<f64>() {
+            Ok(value) => values.push(vec![value]),
+            Err(_) => values.push(parse_point(&item)?),
+        }
+    }
+    Ok(values)
+}
+
+/// The default `LearnerND` loss: the volume of the simplex embedded in
+/// (input + output)-dimensional space, where each vertex is extended with
+/// its function value(s). Signature-compatible with adaptive's
+/// `loss_per_simplex` functions; `value_scale` is accepted and unused, like
+/// the reference (`values` arrive already scaled).
+#[pyfunction]
+#[pyo3(name = "default_loss", signature = (simplex, values, value_scale=None))]
+fn py_default_loss(
+    simplex: &Bound<'_, PyAny>,
+    values: &Bound<'_, PyAny>,
+    value_scale: Option<&Bound<'_, PyAny>>,
+) -> PyResult<f64> {
+    let _ = value_scale;
+    let simplex = parse_points_sized(simplex, "Please provide a 2-dimensional list of points")?;
+    let values = parse_values(values)?;
+    geom::default_loss(&simplex, &values).map_err(|err| PyValueError::new_err(err.to_string()))
+}
+
 #[pyfunction]
 #[pyo3(name = "orientation")]
 fn py_orientation(face: &Bound<'_, PyAny>, origin: &Bound<'_, PyAny>) -> PyResult<i32> {
 
@@ -33,6 +33,22 @@ impl From<TriangulationError> for PyErr {
     }
 }
 
+/// Error mapping for point-in-simplex checks: the reference implementation
+/// divides by the simplex volume and solves a linear system, so degenerate
+/// simplices raise `ZeroDivisionError` and singular solves raise
+/// `numpy.linalg.LinAlgError`.
+fn point_in_simplex_error(py: Python<'_>, error: TriangulationError) -> PyErr {
+    match error {
+        TriangulationError::Geometry(GeometryError::DegenerateSimplex) => {
+            PyZeroDivisionError::new_err("division by zero")
+        }
+        TriangulationError::Geometry(GeometryError::SingularMatrix) => {
+            numpy_linalg_error(py, "Singular matrix")
+        }
+        other => other.into(),
+    }
+}
+
 pub(crate) fn parse_point(obj: &Bound<'_, PyAny>) -> PyResult<Vec<f64>> {
     // Bulk-copy f64 numpy arrays instead of iterating Python objects.
     if let Ok(array) = obj.extract::<PyReadonlyArray1<'_, f64>>() {
@@ -875,16 +891,51 @@ impl PyTriangulation {
     ) -> PyResult<bool> {
         let point = parse_point(point)?;
         let simplex = ordered_indices_from_py(simplex, self.core.vertices.len())?;
-        match self.core.point_in_simplex(&point, &simplex, eps) {
-            Ok(result) => Ok(result),
-            Err(TriangulationError::Geometry(GeometryError::DegenerateSimplex)) => {
-                Err(PyZeroDivisionError::new_err("division by zero"))
-            }
-            Err(TriangulationError::Geometry(GeometryError::SingularMatrix)) => {
-                Err(numpy_linalg_error(py, "Singular matrix"))
+        self.core
+            .point_in_simplex(&point, &simplex, eps)
+            .map_err(|error| point_in_simplex_error(py, error))
+    }
+
+    /// All simplices containing `point`, as a sorted list of tuples: among
+    /// `candidates` when given, otherwise found from a simplex containing
+    /// the point (the `simplex` hint when it holds the point, else
+    /// `locate_point`). One call replaces a Python-side loop of
+    /// `point_in_simplex` checks (the hot loop in adaptive's
+    /// `LearnerND.tell_pending`, whose `simplex` argument maps onto the
+    /// hint).
+    #[pyo3(signature = (point, simplex=None, candidates=None, eps=BARYCENTRIC_EPS))]
+    fn simplices_containing(
+        &self,
+        py: Python<'_>,
+        point: &Bound<'_, PyAny>,
+        simplex: Option<&Bound<'_, PyAny>>,
+        candidates: Option<&Bound<'_, PyAny>>,
+        eps: f64,
+    ) -> PyResult<Py<PyAny>> {
+        let point = parse_point(point)?;
+        let hint = match given(simplex) {
+            None => None,
+            Some(value) => Some(canonical_simplex_from_py(value, self.core.vertices.len())?),
+        };
+        let candidates = match given(candidates) {
+            None => None,
+            Some(value) => {
+                let mut canonical = FxHashSet::default();
+                for simplex in parse_signed_simplex_set(value)? {
+                    canonical.insert(canonicalize_simplex(&simplex, self.core.vertices.len())?);
+                }
+                Some(canonical.into_iter().collect())
             }
-            Err(other) => Err(other.into()),
-        }
+        };
+        let containing = self
+            .core
+            .simplices_containing_point(&point, candidates, hint.as_deref(), eps)
+            .map_err(|error| point_in_simplex_error(py, error))?;
+        let items: Vec<Py<PyAny>> = containing
+            .iter()
+            .map(|simplex| simplex_tuple(py, simplex).into())
+            .collect();
+        Ok(PyList::new(py, items)?.into())
     }
 
     #[pyo3(signature = (pt_index, containing_simplex=None, transform=None))]