diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json
index 6ab23c7..21cf3a0 100644
--- a/.basedpyright/baseline.json
+++ b/.basedpyright/baseline.json
@@ -1207,54 +1207,6 @@
                     "lineCount": 1
                 }
             },
-            {
-                "code": "reportUnknownVariableType",
-                "range": {
-                    "startColumn": 8,
-                    "endColumn": 20,
-                    "lineCount": 1
-                }
-            },
-            {
-                "code": "reportUnknownArgumentType",
-                "range": {
-                    "startColumn": 32,
-                    "endColumn": 13,
-                    "lineCount": 5
-                }
-            },
-            {
-                "code": "reportCallIssue",
-                "range": {
-                    "startColumn": 12,
-                    "endColumn": 24,
-                    "lineCount": 1
-                }
-            },
-            {
-                "code": "reportArgumentType",
-                "range": {
-                    "startColumn": 16,
-                    "endColumn": 23,
-                    "lineCount": 1
-                }
-            },
-            {
-                "code": "reportAny",
-                "range": {
-                    "startColumn": 4,
-                    "endColumn": 13,
-                    "lineCount": 1
-                }
-            },
-            {
-                "code": "reportUnknownArgumentType",
-                "range": {
-                    "startColumn": 23,
-                    "endColumn": 35,
-                    "lineCount": 1
-                }
-            },
             {
                 "code": "reportAny",
                 "range": {
@@ -2481,14 +2433,6 @@
                     "lineCount": 1
                 }
             },
-            {
-                "code": "reportImplicitAbstractClass",
-                "range": {
-                    "startColumn": 6,
-                    "endColumn": 19,
-                    "lineCount": 1
-                }
-            },
             {
                 "code": "reportUnknownParameterType",
                 "range": {
@@ -6273,6 +6217,78 @@
                     "lineCount": 1
                 }
             },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 23,
+                    "endColumn": 82,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 27,
+                    "endColumn": 54,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 34,
+                    "endColumn": 46,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 27,
+                    "endColumn": 54,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 34,
+                    "endColumn": 46,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 27,
+                    "endColumn": 54,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 34,
+                    "endColumn": 46,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 33,
+                    "endColumn": 47,
+                    "lineCount": 1
+                }
+            },
+            {
+                "code": "reportAny",
+                "range": {
+                    "startColumn": 40,
+                    "endColumn": 46,
+                    "lineCount": 1
+                }
+            },
             {
                 "code": "reportAny",
                 "range": {
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 3bde834..2be23f6 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -32,7 +32,7 @@ jobs:
                 python-version: '3.x'
         -   name: "Main Script"
             run: |
-                EXTRA_INSTALL="matplotlib pytest scipy scipy-stubs basedpyright"
+                EXTRA_INSTALL="matplotlib pytest scipy scipy-stubs basedpyright optype"
                 curl -L -O https://tiker.net/ci-support-v0
                 . ./ci-support-v0
                 build_py_project_in_venv
diff --git a/doc/conf.py b/doc/conf.py
index 9dbc49a..45311f4 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -23,6 +23,11 @@
 
 autodoc_member_order = "bysource"
 
+nitpick_ignore_regex = [
+    # NOTE: optype does not have Sphinx compatible documentation
+    ["py:class", r"onp.*"],
+]
+
 sphinxconfig_missing_reference_aliases = {
     # numpy
     "NDArray": "obj:numpy.typing.NDArray",
diff --git a/modepy/__init__.py b/modepy/__init__.py
index 6f22b14..7da8902 100644
--- a/modepy/__init__.py
+++ b/modepy/__init__.py
@@ -45,6 +45,7 @@
 from modepy.modes import (
     Basis,
     BasisNotOrthonormal,
+    SimplexBasis,
     TensorProductBasis,
     basis_for_space,
     grad_jacobi,
@@ -133,6 +134,7 @@
     "QuadratureRuleUnavailable",
     "Shape",
     "Simplex",
+    "SimplexBasis",
     "TensorProductBasis",
     "TensorProductQuadrature",
     "TensorProductShape",
diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py
index b278d13..ecc1c58 100644
--- a/modepy/modal_decay.py
+++ b/modepy/modal_decay.py
@@ -23,14 +23,18 @@
 THE SOFTWARE.
 """
 
-
+from functools import singledispatch
 from typing import TYPE_CHECKING
 
 import numpy as np
 import numpy.linalg as la
+from typing_extensions import deprecated
+
+from modepy.modes import Basis, SimplexBasis, TensorProductBasis
 
 
 if TYPE_CHECKING:
+    import optype.numpy as onp
     from numpy.typing import NDArray
 
     from modepy.typing import ArrayF
@@ -49,25 +53,82 @@
 
 .. versionadded:: 2013.2
 
-.. autofunction:: simplex_interp_error_coefficient_estimator_matrix
+.. autofunction:: degree_vector
+.. autofunction:: interp_error_coefficient_estimator_matrix
 .. autofunction:: fit_modal_decay
 .. autofunction:: estimate_relative_expansion_residual
 """
 
 
-def simplex_interp_error_coefficient_estimator_matrix(
-        unit_nodes: ArrayF,
-        order: int,
-        n_tail_orders: int) -> ArrayF:
+@singledispatch
+def degree_vector(basis: Basis) -> onp.Array1D[np.integer]:
+    """Return a vector of integers indicating the 'overall order' of each basis
+    function in the basis, in order. For simplex discretizations, consider
+    the total degree. For tensor product discretizations, the max degree.
+
+    .. versionadded:: 2026.1.1
+    """
+    raise TypeError(f"unsupported basis type: {type(basis)}")
+
+
+@degree_vector.register(SimplexBasis)
+def degree_vector_simplex(basis: SimplexBasis) -> onp.Array1D[np.integer]:
+    return np.array([
+                        sum(mid_tuple)
+                        for mid_tuple in basis.mode_ids
+                    ])
+
+
+@degree_vector.register(TensorProductBasis)
+def degree_vector_tp(basis: TensorProductBasis) -> onp.Array1D[np.integer]:
+    degree_vectors = [degree_vector(b) for b in basis.bases]
+    return np.array([
+                        max(
+                            ov[mode_id]
+                            for mode_id, ov
+                            in zip(mid_tuple, degree_vectors, strict=True)
+                        )
+                        for mid_tuple in basis._mode_index_tuples  # pyright: ignore[reportPrivateUsage]
+                    ])
+
+
+def interp_error_coefficient_estimator_matrix(
+            unit_nodes: ArrayF,
+            n_tail_orders: int,
+            basis: Basis,
+        ) -> ArrayF:
     """Return a matrix :math:`C` that, when multiplied by a vector of nodal values,
-    yields the coeffiicients belonging to the basis functions of the *n_tail_orders*
+    yields the coefficients belonging to the basis functions of the *n_tail_orders*
     highest orders.
 
     The 2-norm of the resulting coefficients can be used as an estimate of the
     interpolation error.
 
-    .. versionadded:: 2018.1
+    .. versionadded:: 2026.1.1
     """
+    from modepy.matrices import vandermonde
+    vdm = vandermonde(basis.functions, unit_nodes)
+    vdm_inv = la.inv(vdm)
+
+    degree_vec = degree_vector(basis)
+
+    max_order = np.max(degree_vec)
+    return vdm_inv[degree_vec > max_order-n_tail_orders]
+
+
+@deprecated("Use interp_error_coefficient_estimator_matrix instead")
+def simplex_interp_error_coefficient_estimator_matrix(
+            unit_nodes: ArrayF,
+            order: int,
+            n_tail_orders: int,
+        ) -> ArrayF:
+    from warnings import warn
+    warn("simplex_interp_error_coefficient_estimator_matrix is deprecated, "
+         "use interp_error_coefficient_estimator_matrix instead. "
+         "This will stop working in 2027.x.",
+         DeprecationWarning, stacklevel=2,
+     )
+
     dim, _nunit_nodes = unit_nodes.shape
 
     from modepy.shapes import Simplex
@@ -77,24 +138,8 @@ def simplex_interp_error_coefficient_estimator_matrix(
 
     from modepy.modes import orthonormal_basis_for_space
     basis = orthonormal_basis_for_space(space, shape)
-
-    from modepy.matrices import vandermonde
-    vdm = vandermonde(basis.functions, unit_nodes)
-    vdm_inv = la.inv(vdm)
-
-    if dim > 1:
-        order_vector = np.array([
-            # NOTE: basis.mode_ids are declared as `Tuple[Hashable, ...]`, which
-            # cannot be summed
-            sum(mode_id) for mode_id in basis.mode_ids
-            ])
-    else:
-        order_vector = np.array(basis.mode_ids)
-
-    max_order = np.max(order_vector)
-    assert max_order == order
-
-    return vdm_inv[order_vector > max_order-n_tail_orders]
+    return interp_error_coefficient_estimator_matrix(
+                        unit_nodes, n_tail_orders, basis)
 
 
 def make_mode_number_vector(
diff --git a/modepy/modes.py b/modepy/modes.py
index 263ce47..9374a63 100644
--- a/modepy/modes.py
+++ b/modepy/modes.py
@@ -83,6 +83,10 @@
 
 .. autofunction:: symbolicize_function
 
+Base classes
+------------
+.. autoclass:: SimplexBasis
+
 Tensor product adapter
 ----------------------
 
@@ -873,7 +877,7 @@ def _grad_pkdo_1d(order: tuple[int], r: ArrayF) -> tuple[ArrayF]:
     return (grad_jacobi(0, 0, i, r0),)
 
 
-class _SimplexBasis(Basis):
+class SimplexBasis(Basis, ABC):
     def __init__(self, dim, order):
         self._dim = dim
         self._order = order
@@ -889,7 +893,7 @@ def mode_ids(self):
         return tuple(gnitsam(self._order, self._dim))
 
 
-class _SimplexONB(_SimplexBasis):
+class _SimplexONB(SimplexBasis):
     is_orthonormal = True
 
     def orthonormality_weight(self):
@@ -920,7 +924,7 @@ def gradients(self):
             raise NotImplementedError(f"gradient in {self._dim} dimensions")
 
 
-class _SimplexMonomialBasis(_SimplexBasis):
+class _SimplexMonomialBasis(SimplexBasis):
     def orthonormality_weight(self) -> float:
         raise BasisNotOrthonormal
 
diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py
index 2d8d781..22fdc7e 100644
--- a/modepy/test/test_tools.py
+++ b/modepy/test/test_tools.py
@@ -157,6 +157,99 @@ def estimate_resid(inner_n: int) -> ArrayF:
 # }}}
 
 
+# {{{ test_interp_error_coefficient_estimator_matrix
+
+@pytest.mark.parametrize("shape", [
+    mp.Simplex(1),
+    mp.Simplex(2),
+    mp.Simplex(3),
+    mp.Hypercube(1),
+    mp.Hypercube(2),
+    mp.Hypercube(3),
+])
+def test_interp_error_coefficient_estimator_matrix(shape: mp.Shape) -> None:
+    """Check that interp_error_coefficient_estimator_matrix gives estimates that
+    track the actual interpolation error for functions of varying smoothness.
+
+    For each shape the test:
+    - Builds the estimator matrix C at polynomial order 7.
+    - For four test functions (C^1, C^2, C^3, analytic) computes both the
+      estimated error ``||C @ f_nodes||_2`` and the actual RMS interpolation
+      error sampled at many random points inside the reference element.
+    - Asserts that estimated and actual errors are within a factor of 100 of
+      each other.
+    - Asserts that estimated errors are ordered correctly by smoothness:
+      C^1 > C^2 > C^3 > analytic.
+
+    Test functions use ``|x - 0.3|^alpha`` with fractional exponents alpha,
+    which are C^{floor(alpha)} but not C^{ceil(alpha)}, evaluated on the first
+    reference coordinate.  Fractional powers avoid the piecewise-polynomial
+    coincidences that cause accidental zero modal coefficients.
+    """
+    from modepy.modal_decay import interp_error_coefficient_estimator_matrix
+
+    order = 7
+    space = mp.space_for_shape(shape, order)
+    nodes = mp.edge_clustered_nodes_for_space(space, shape)
+    basis = mp.orthonormal_basis_for_space(space, shape)
+
+    # Build the estimator matrix (one tail order = highest-degree modes only)
+    err_est_matrix = interp_error_coefficient_estimator_matrix(nodes, 1, basis)
+
+    # Dense random sample for measuring actual interpolation error.
+    # Using random sampling instead of quadrature avoids availability and
+    # negative-weight issues with high-order rules for 3D simplices.
+    n_fine = 2000
+    rng = np.random.default_rng(seed=42)
+    fine_nodes = mp.random_nodes_for_shape(shape, n_fine, rng)
+    resamp = mp.resampling_matrix(basis.functions, fine_nodes, nodes)
+
+    def _errors(f: Callable[[ArrayF], ArrayF]) -> tuple[float, float]:
+        f_nodes = f(nodes)
+        estimated = float(la.norm(err_est_matrix @ f_nodes))
+        # RMS error at random fine nodes (normalised by sqrt(n_fine) so
+        # it approximates the L2 error per unit volume)
+        actual = float(la.norm(f(fine_nodes) - resamp @ f_nodes) / np.sqrt(n_fine))
+        return estimated, actual
+
+    # |x - 0.3|^alpha with fractional alpha: C^k = C^{floor(alpha)}, not C^{k+1}.
+    # Kink at x = 0.3, which lies inside every reference element used here.
+    test_cases: list[tuple[str, Callable[[ArrayF], ArrayF]]] = [
+        ("c1", lambda pts: np.abs(pts[0] - 0.3) ** 1.5),
+        ("c2", lambda pts: np.abs(pts[0] - 0.3) ** 2.5),
+        ("c3", lambda pts: np.abs(pts[0] - 0.3) ** 3.5),
+        ("analytic", lambda pts: np.exp(pts[0])),
+    ]
+
+    estimated_errors: dict[str, float] = {}
+
+    for name, f in test_cases:
+        est, act = _errors(f)
+        estimated_errors[name] = est
+        logger.info("shape=%s %s: estimated=%.5e actual=%.5e",
+                    shape, name, est, act)
+
+        # Estimated and actual errors should be within a factor of 100
+        if act > 1e-14:
+            ratio = est / act
+            assert 0.05 < ratio < 100, (
+                f"shape={shape} {name}: estimated/actual={ratio:.3e} "
+                f"out of [0.05, 100]"
+            )
+        else:
+            assert est < 1e-10, (
+                f"shape={shape} {name}: actual error tiny ({act:.2e}) "
+                f"but estimated is {est:.2e}"
+            )
+
+    # Smoother functions should yield strictly smaller estimated errors
+    assert estimated_errors["c1"] > estimated_errors["c2"], shape
+    assert estimated_errors["c2"] > estimated_errors["c3"], shape
+    assert estimated_errors["c3"] > estimated_errors["analytic"], shape
+
+# }}}
+
+
 # {{{ test_resampling_matrix
 
 @pytest.mark.parametrize("dims", [1, 2, 3])
diff --git a/pyproject.toml b/pyproject.toml
index 709f8ce..c633fd0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,6 +35,7 @@ dependencies = [
     "pytools",
     # for Never
     "typing_extensions>=4.1",
+    "optype",
 ]
 
 [project.optional-dependencies]