From e26da9ddce20440b1af3d159f7d035fb869f11bc Mon Sep 17 00:00:00 2001
From: aochuba <aochuba52@gmail.com>
Date: Thu, 18 Jun 2026 23:47:45 +0530
Subject: [PATCH 1/4] Deprecate InverseRTransform and migrate ODE/Poisson
 solvers to forward transforms.  adjusted tests.

---
 src/grid/ode.py                |  46 ++++++++------
 src/grid/poisson.py            |  24 ++++----
 src/grid/rtransform.py         |  16 ++++-
 src/grid/tests/test_ode.py     | 106 ++++++++++++++++-----------------
 src/grid/tests/test_poisson.py |   9 ++-
 5 files changed, 111 insertions(+), 90 deletions(-)

diff --git a/src/grid/ode.py b/src/grid/ode.py
index a5219bd8f..1c055bd8d 100644
--- a/src/grid/ode.py
+++ b/src/grid/ode.py
@@ -109,9 +109,8 @@ def func(x, y):
         # x has shape (1,) and y has shape (K+1,1), output has shape (K+1,1)
         x = np.array([x])
         if transform:
-            # Transform the points back to the original domain.
-            orig_dom = transform.inverse(x)
-            dy_dx = _transform_and_rearrange_to_explicit_ode(orig_dom, y, coeffs, transform, fx)
+            # x here is in the finite domain; pass it directly.
+            dy_dx = _transform_and_rearrange_to_explicit_ode(x, y, coeffs, transform, fx)
         else:
             coeffs_mt = _evaluate_coeffs_on_points(x, coeffs)
             dy_dx = _rearrange_to_explicit_ode(y, coeffs_mt, fx(x))
@@ -120,21 +119,21 @@ def func(x, y):
         return np.vstack((*y[1:, :], dy_dx))
 
     if transform:
-        # first check if the bounds are in the domain
-        if min(x_span) < transform.domain[0] or max(x_span) > transform.domain[1]:
+        # first check if the bounds are in the codomain (infinite domain r ∈ [0,∞))
+        if min(x_span) < transform.codomain[0] or max(x_span) > transform.codomain[1]:
             raise ValueError(
                 f"The x_span {min(x_span), max(x_span)} is not within the transform "
                 f"domain {transform.domain}."
             )
         # Convert the initial value problem to the new derivative space, only transform up to K-1
-        # e.g. the first derivative dV/dx = dV/dr * dr/dx = dV/dr / (dx/dr)
+        # e.g. the first derivative dV/dx = dV/dr * dr/dx = dV/dr * deriv_inverse
         deriv = _derivative_transformation_matrix(
-            [transform.deriv, transform.deriv2, transform.deriv3],
+            [transform.deriv_inverse, transform.deriv2_inverse, transform.deriv3_inverse],
             x_span[0],
             order - 1,  # Only need derivatives up to K-1.
         )
-        # If transform is used, then transform (x_0, x_1) that it integrates up to.
-        x_span = transform.transform(np.array(list(x_span)))
+        # If transform is used, map (r0, r1) to finite domain (x0, x1) via inverse.
+        x_span = transform.inverse(np.array(list(x_span)))
         # Solve for derivatives in original domain by solving A(original derivs) = new derivs
         y_derivs = solve(deriv, np.array(y0[1:]))
         if np.any(np.isinf(y_derivs)):
@@ -236,9 +235,8 @@ def solve_ode_bvp(
     def func(x, y):
         # x has shape (N,) and y has shape (K+1, N), output has shape (K+1, N)
         if transform:
-            # Transform the points back to the original domain.
-            orig_dom = transform.inverse(x)
-            dy_dx = _transform_and_rearrange_to_explicit_ode(orig_dom, y, coeffs, transform, fx)
+            # x here is in the finite domain; pass it directly.
+            dy_dx = _transform_and_rearrange_to_explicit_ode(x, y, coeffs, transform, fx)
         else:
             coeffs_mt = _evaluate_coeffs_on_points(x, coeffs)
             dy_dx = _rearrange_to_explicit_ode(y, coeffs_mt, fx(x))
@@ -262,7 +260,8 @@ def bc(ya, yb):
 
     # Solve the ODE
     if transform:
-        pts_tf = transform.transform(x)
+        # Map the radial points r to finite domain x via the forward transform's inverse.
+        pts_tf = transform.inverse(x)
         res = solve_bvp(func, bc, pts_tf, y=initial_guess_y, tol=tol, max_nodes=max_nodes)
     else:
         res = solve_bvp(func, bc, x, y=initial_guess_y, tol=tol, max_nodes=max_nodes)
@@ -284,7 +283,8 @@ def _transform_solution_to_original_domain(result, tf, no_derivs, order):
 
     # Note this is it's own function becuase it is used twice for solve_ode_ivp and bv.
     def interpolate_wrt_original_var(pt):
-        transf_pts = tf.transform(pt)
+        # Map radial points r to finite domain x via the forward transform's inverse.
+        transf_pts = tf.inverse(pt)
         # Row is which func/deriv and Col is points.
         interpolated = result.sol(transf_pts)
         # If derivatives are not wanted then only return y(x).
@@ -292,11 +292,14 @@ def interpolate_wrt_original_var(pt):
             if interpolated.ndim == 1:
                 return interpolated
             return interpolated[0, :]
-        deriv_funcs = [tf.deriv, tf.deriv2, tf.deriv3]
+        # Use inverse derivative methods from BaseTransform (added in #307).
+        # The ODE is solved in the transformed domain x; converting back to r requires
+        # dx/dr derivatives, i.e. the derivatives of the inverse transformation.
+        deriv_funcs = [tf.deriv_inverse, tf.deriv2_inverse, tf.deriv3_inverse]
         new_interpolate = np.zeros(interpolated.shape)
         new_interpolate[0, :] = interpolated[0, :]
         for i in range(interpolated.shape[1]):
-            # Calculate the jacobian dr/dx of the original domain.
+            # Calculate the jacobian dx/dr of the original domain.
             deriv = _derivative_transformation_matrix(deriv_funcs, pt[i], order - 1)
             new_interpolate[1:, i] = deriv.dot(interpolated[1:, i])
         return new_interpolate
@@ -407,8 +410,12 @@ def _transform_ode_from_rtransform(coeff_a: list | np.ndarray, tf: BaseTransform
         Coefficients :math:`b_j(r)` of the new ODE with respect to transformed variable :math:`r`.
 
     """
-    deriv_func = [tf.deriv, tf.deriv2, tf.deriv3]
-    return _transform_ode_from_derivs(coeff_a, deriv_func, x)
+    # x here is the finite domain points of tf.
+    # We want to transform the ODE from the infinite domain r = tf(x) to the finite domain x.
+    r = tf.transform(x)
+    # The derivatives we need are the derivatives of the inverse transform (dr -> dx), i.e. tf.deriv_inverse
+    deriv_func = [tf.deriv_inverse, tf.deriv2_inverse, tf.deriv3_inverse]
+    return _transform_ode_from_derivs(coeff_a, deriv_func, r)
 
 
 def _transform_and_rearrange_to_explicit_ode(
@@ -442,8 +449,9 @@ def _transform_and_rearrange_to_explicit_ode(
         to explicit form evaluated on all N points.
 
     """
+    orig_dom = tf.transform(x)
     coeff_b = _transform_ode_from_rtransform(coeff_a, tf, x)
-    result = _rearrange_to_explicit_ode(y, coeff_b, fx_func(x))
+    result = _rearrange_to_explicit_ode(y, coeff_b, fx_func(orig_dom))
     return result
 
 
diff --git a/src/grid/poisson.py b/src/grid/poisson.py
index 8b7871ee1..8cfd95a76 100644
--- a/src/grid/poisson.py
+++ b/src/grid/poisson.py
@@ -194,9 +194,11 @@ def solve_poisson_ivp(
         spherical harmonic basis.
     func_vals : ndarray(N,)
         The function values evaluated on all :math:`N` points on the molecular grid.
-    transform : BaseTransform, optional
-        Transformation from infinite domain :math:`r \in [0, \infty)` to another
-        domain that is a finite.
+    transform : BaseTransform
+        Forward transformation from the finite domain :math:`x` to the infinite radial
+        domain :math:`r \in [0, \infty)`, i.e. the same transform used to construct the
+        atomic grid. Inverse-related quantities are computed internally via
+        ``BaseTransform.deriv_inverse`` and related methods.
     r_interval : tuple, optional
         The interval :math:`(b, a)` of :math:`r` for which the ODE solver will start from and end,
         where :math:`b>a`\. The value :math:`b` should be large as it determines the asymptotic
@@ -244,10 +246,10 @@ def _solve_poisson_bvp_atomgrid(
         sph_o_l = generate_real_spherical_harmonics(0, np.array([0.1]), np.array([0.1]))
         boundary = atomgrid.integrate(func_vals) / sph_o_l[0, 0]
 
-    # Check if the domain of transform is in [0, \infty)
-    domain = transform.domain
-    if domain[0] < 0.0:
-        raise ValueError(f"The domain of the transform {domain} should be in [0, infinity).")
+    # Check if the codomain of transform is in [0, \infty)
+    codomain = transform.codomain
+    if codomain[0] < 0.0:
+        raise ValueError(f"The codomain of the transform {codomain} should be in [0, infinity).")
 
     # Get the radial components from expanding func into real spherical harmonics.
     radial_components = atomgrid.radial_component_splines(func_vals)
@@ -346,9 +348,11 @@ def solve_poisson_bvp(
         harmonic basis.
     func_vals : ndarray(N,)
         The function values evaluated on all :math:`N` points on the molecular grid.
-    transform : BaseTransform, optional
-        Transformation from infinite domain :math:`r \in [0, \infty)` to another
-        domain that is a finite.
+    transform : BaseTransform
+        Forward transformation from the finite domain :math:`x` to the infinite radial
+        domain :math:`r \in [0, \infty)`, i.e. the same transform used to construct the
+        atomic grid. Inverse-related quantities are computed internally via
+        ``BaseTransform.deriv_inverse`` and related methods.
     boundary : float, optional
         The boundary value of :math:`g` in the limit of r to infinity.
     include_origin : bool, optional
diff --git a/src/grid/rtransform.py b/src/grid/rtransform.py
index 56daa9fb6..372ccc192 100644
--- a/src/grid/rtransform.py
+++ b/src/grid/rtransform.py
@@ -515,7 +515,14 @@ def inverse(self, r: np.ndarray):
 
 
 class InverseRTransform(BaseTransform):
-    """Inverse transformation class for any general transformation."""
+    """Inverse transformation class for any general transformation.
+
+    .. deprecated::
+        ``InverseRTransform`` is deprecated and will be removed in a future release.
+        The inverse transformation derivatives (``deriv_inverse``, ``deriv2_inverse``,
+        ``deriv3_inverse``) are now available directly on every ``BaseTransform`` subclass.
+        Pass the forward transform and call those methods instead.
+    """
 
     def __init__(self, transform: BaseTransform):
         """Construct InverseRTransform instance.
@@ -526,6 +533,13 @@ def __init__(self, transform: BaseTransform):
             One-dimension transformation instance.
 
         """
+        warnings.warn(
+            "InverseRTransform is deprecated and will be removed in a future release. "
+            "Use the forward transform directly and call deriv_inverse(), deriv2_inverse(), "
+            "or deriv3_inverse() on it instead.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
         if not isinstance(transform, BaseTransform):
             raise TypeError(f"Input need to be a transform instance, got {type(transform)}.")
         self._tfm = transform
diff --git a/src/grid/tests/test_ode.py b/src/grid/tests/test_ode.py
index 11b3f1eb8..43b2d43ef 100644
--- a/src/grid/tests/test_ode.py
+++ b/src/grid/tests/test_ode.py
@@ -41,7 +41,6 @@
     BaseTransform,
     BeckeRTransform,
     IdentityRTransform,
-    InverseRTransform,
     KnowlesRTransform,
     LinearFiniteRTransform,
 )
@@ -80,10 +79,15 @@ def __init__(self, exp=2, extra=0):
         """Initialize power transform instance."""
         self._exp = exp
         self._extra = extra
+        self._domain = (-1.0, 1.0)
+        if exp % 2 == 0:
+            self._codomain = (extra, 1.0 + extra)
+        else:
+            self._codomain = (-1.0 + extra, 1.0 + extra)
 
     @property
     def domain(self):
-        return (-1.0, 1.0)
+        return self._domain
 
     def transform(self, x):
         """Transform given array."""
@@ -112,7 +116,7 @@ def deriv3(self, x):
         [IdentityRTransform(), fx_ones, [-1, 1, 1]],
         [SqTF(), fx_ones, np.random.uniform(-100, 100, (3,))],
         [KnowlesRTransform(0.01, 1, 5), fx_ones, np.random.uniform(-10, 10, (3,))],
-        [BeckeRTransform(0.1, 100.0), fx_ones, np.random.uniform(0, 100, (3,))],
+        [BeckeRTransform(0.1, 2.0), fx_ones, np.random.uniform(0, 100, (3,))],
         [SqTF(1, 3), fx_complicated_example, np.random.uniform(0, 100, (3,))],
         [SqTF(3, 1), fx_complicated_example, [2, 3, 2]],
         [SqTF(), fx_quadratic, np.random.uniform(-100, 100, (3,))],
@@ -122,59 +126,55 @@ def deriv3(self, x):
 )
 def test_transform_and_rearrange_to_explicit_ode_with_simple_boundary(transform, fx, coeffs):
     r"""Test transforming second-order ode with simple boundary conditions."""
-    x = np.arange(0.1, 0.99, 0.01)
+    x = np.linspace(0.1, 0.5, 90)
     transform_pts = transform.transform(x)
 
     def bc(ya, yb):
         # Boundary of y is zero on endpoints
         return np.array([ya[0], yb[0]])
 
-    # Run with transformation
-    def func_with_transform(r, y):
-        # Transform back to original domain x
-        original = transform.inverse(r)
-        # Apply the ode transfomration
-        dy_dx = _transform_and_rearrange_to_explicit_ode(original, y, coeffs, transform, fx)
+    # Run with transformation (integrates on finite domain x)
+    def func_with_transform(x_pts, y):
+        dy_dx = _transform_and_rearrange_to_explicit_ode(x_pts, y, coeffs, transform, fx)
         return np.vstack((*y[1:], dy_dx))
 
     init_guess = np.zeros((2, x.size))
     solution_with_transf = solve_bvp(
-        func_with_transform, bc, transform_pts, init_guess, tol=1e-6, max_nodes=10000
+        func_with_transform, bc, x, init_guess, tol=1e-6, max_nodes=10000
     )
 
-    # Run without transformation
+    # Run without transformation (integrates on infinite domain r)
     def func_without_transform(original_pts, y):
         coeffs_mt = _evaluate_coeffs_on_points(original_pts, coeffs)
         dy_dx = _rearrange_to_explicit_ode(y, coeffs_mt, fx(original_pts))
         return np.vstack((y[1:], dy_dx))
 
-    # Run without transformation
     solution_without_transf = solve_bvp(
         func_without_transform,
         bc,
-        x,
-        solution_with_transf.sol(transform_pts),
+        transform_pts,
+        solution_with_transf.sol(x),
         tol=1e-6,
         max_nodes=100000,
     )
     # Check if the solution at x is the same as the solution (with transform) at r=g(x)
     assert_allclose(
-        solution_with_transf.sol(transform_pts)[0],
-        solution_without_transf.sol(x)[0],
+        solution_with_transf.sol(x)[0],
+        solution_without_transf.sol(transform_pts)[0],
         atol=1e-4,
     )
     # Check if they're similar at random points on interval with lower accuracy tolerance
-    random_pts = np.random.uniform(np.min(x), np.max(x), transform_pts.shape)
+    random_pts = np.random.uniform(np.min(x), np.max(x), x.shape)
     transf_pts = transform.transform(random_pts)
     assert_allclose(
-        solution_with_transf.sol(transf_pts)[0],
-        solution_without_transf.sol(random_pts)[0],
+        solution_with_transf.sol(random_pts)[0],
+        solution_without_transf.sol(transf_pts)[0],
         atol=1e-4,
     )
-    # Test the derivative
+    # Test the derivative: dy/dr = dy/dx / transform.deriv(x)
     assert_allclose(
-        solution_with_transf.sol(transform_pts)[1],
-        solution_without_transf.sol(x)[1] / transform.deriv(x),
+        solution_with_transf.sol(x)[1] / transform.deriv(x),
+        solution_without_transf.sol(transform_pts)[1],
         atol=1e-4,
     )
 
@@ -189,12 +189,6 @@ def func_without_transform(original_pts, y):
             [-1, 1, 1],
             [(0, 0, 3), (1, 0, 3)],
         ],
-        [
-            InverseRTransform(BeckeRTransform(1.0, 5.0)),
-            fx_complicated_example3,
-            [-1, 1, 1],
-            [(0, 0, 3), (1, 0, 3)],
-        ],
         [
             BeckeRTransform(1.0, 5.0),
             fx_complicated_example3,
@@ -212,9 +206,11 @@ def func_without_transform(original_pts, y):
 )
 def test_solve_ode_bvp_with_and_without_transormation(transform, fx, coeffs, bd_cond):
     r"""Test solve_ode with and without transformation with different bd conditions."""
-    x = np.linspace(0.01, 0.999, 20)
+    # Use r-domain points (codomain of the forward transform): transform x → r
+    x_domain = np.linspace(-0.9, 0.9, 20)
+    r_pts = transform.transform(x_domain)  # r ∈ [r_min, ∞)
     sol_with_transform = solve_ode_bvp(
-        x,
+        r_pts,
         fx,
         coeffs,
         bd_cond,
@@ -223,23 +219,23 @@ def test_solve_ode_bvp_with_and_without_transormation(transform, fx, coeffs, bd_
         max_nodes=20000,
         no_derivatives=False,
     )
-    init_guess = sol_with_transform(x)
+    init_guess = sol_with_transform(r_pts)
 
     sol_normal = solve_ode_bvp(
-        x, fx, coeffs, bd_cond, tol=1e-8, max_nodes=20000, initial_guess_y=init_guess
+        r_pts, fx, coeffs, bd_cond, tol=1e-8, max_nodes=20000, initial_guess_y=init_guess
     )
     # Test the function values
-    assert_allclose(sol_with_transform(x)[0], sol_normal(x)[0], atol=1e-5)
+    assert_allclose(sol_with_transform(r_pts)[0], sol_normal(r_pts)[0], atol=1e-5)
     # Test the boundary condition
     for bd in bd_cond:
         bnd, deriv, val = bd
-        assert_allclose(sol_with_transform(x)[deriv][-bnd], val, atol=1e-5)
+        assert_allclose(sol_with_transform(r_pts)[deriv][-bnd], val, atol=1e-5)
 
     if len(coeffs) >= 3:
         # Test the first derivative of y.
-        assert_allclose(sol_with_transform(x)[1], sol_normal(x)[1], atol=1e-3)
+        assert_allclose(sol_with_transform(r_pts)[1], sol_normal(r_pts)[1], atol=1e-3)
         if len(coeffs) >= 4:
-            assert_allclose(sol_with_transform(x)[2], sol_normal(x)[2], atol=1e-3)
+            assert_allclose(sol_with_transform(r_pts)[2], sol_normal(r_pts)[2], atol=1e-3)
 
 
 @pytest.mark.parametrize(
@@ -281,29 +277,31 @@ def test_solve_ode_bvp_with_and_without_transormation(transform, fx, coeffs, bd_
 )
 def test_solve_ode_ivp_with_and_without_transformation(transform, fx, coeffs, ivp):
     r"""Test solve_ode_ivp with and without transformation with different initial guesses."""
+    r_start = transform.transform(0.01)
+    r_end = transform.transform(0.5)
     sol_with_transform = solve_ode_ivp(
-        (0.01, 0.999),
+        (r_start, r_end),
         fx,
         coeffs,
         ivp,
         transform,
     )
 
-    sol_normal = solve_ode_ivp((0.01, 0.999), fx, coeffs, ivp)
+    sol_normal = solve_ode_ivp((r_start, r_end), fx, coeffs, ivp)
 
     # Test the initial value problem
     for i, ivp_val in enumerate(ivp):
-        assert_allclose(sol_normal(np.array([0.01]))[i], ivp_val, atol=1e-5)
-        assert_allclose(sol_with_transform(np.array([0.01]))[i], ivp_val, atol=1e-5)
+        assert_allclose(sol_normal(np.array([r_start]))[i], ivp_val, atol=1e-5)
+        assert_allclose(sol_with_transform(np.array([r_start]))[i], ivp_val, atol=1e-5)
 
     # Test the function values
-    x = np.arange(0.01, 0.999, 0.1)
-    assert_allclose(sol_with_transform(x)[0], sol_normal(x)[0], atol=1e-5, rtol=1e-5)
+    r_pts = transform.transform(np.arange(0.01, 0.5, 0.05))
+    assert_allclose(sol_with_transform(r_pts)[0], sol_normal(r_pts)[0], atol=1e-5, rtol=1e-5)
     if len(coeffs) >= 3:
         # Test the first derivative of y.
-        assert_allclose(sol_with_transform(x)[1], sol_normal(x)[1], atol=1e-3, rtol=1e-5)
+        assert_allclose(sol_with_transform(r_pts)[1], sol_normal(r_pts)[1], atol=1e-3, rtol=1e-5)
         if len(coeffs) >= 4:
-            assert_allclose(sol_with_transform(x)[2], sol_normal(x)[2], atol=1e-3, rtol=1e-5)
+            assert_allclose(sol_with_transform(r_pts)[2], sol_normal(r_pts)[2], atol=1e-3, rtol=1e-5)
 
 
 @pytest.mark.parametrize(
@@ -429,14 +427,14 @@ def test_transform_coeff_with_x_and_r():
     """Test coefficient transform between x and r."""
     coeff = np.array([2, 3, 4])
     ltf = LinearFiniteRTransform(1, 10)  # (-1, 1) -> (r0, rmax)
-    inv_tf = InverseRTransform(ltf)  # (r0, rmax) -> (-1, 1)
     x = np.linspace(-1, 1, 20)
     r = ltf.transform(x)
     assert r[0] == 1
     assert r[-1] == 10
-    # Transform ODE from [1, 10) to (-1, 1)
-    coeff_transform = _transform_ode_from_rtransform(coeff, inv_tf, x)
-    derivs_fun = [inv_tf.deriv, inv_tf.deriv2, inv_tf.deriv3]
+    # Transform ODE from r ∈ [1, 10] to x ∈ [-1, 1] using the forward transform and x points
+    coeff_transform = _transform_ode_from_rtransform(coeff, ltf, x)
+    # Equivalent: evaluate inverse derivatives at r-domain points
+    derivs_fun = [ltf.deriv_inverse, ltf.deriv2_inverse, ltf.deriv3_inverse]
     coeff_transform_all_pts = _transform_ode_from_derivs(coeff, derivs_fun, r)
     assert_allclose(coeff_transform, coeff_transform_all_pts)
 
@@ -444,10 +442,9 @@ def test_transform_coeff_with_x_and_r():
 def test_transformation_of_ode_with_identity_transform():
     """Test transformation of ODE with identity transform."""
     # Checks that the identity transform x -> x results in the same answer.
-    # Obtain identity trasnform and derivatives.
+    # Use inverse derivative methods from BaseTransform directly.
     itf = IdentityRTransform()
-    inv_tf = InverseRTransform(itf)
-    derivs_fun = [inv_tf.deriv, inv_tf.deriv2, inv_tf.deriv3]
+    derivs_fun = [itf.deriv_inverse, itf.deriv2_inverse, itf.deriv3_inverse]
     # d^2y / dx^2 = 1
     coeff = np.array([0, 0, 1])
     x = np.linspace(0, 1, 10)
@@ -462,9 +459,8 @@ def test_transformation_of_ode_with_linear_transform():
     x = GaussLaguerre(10).points
     # Obtain linear transformation with rmin = 1 and rmax = 10.
     ltf = LinearFiniteRTransform(1, 10)
-    # The inverse is x_i = \frac{r_i - r_{min} - R} {r_i - r_{min} + R}
-    inv_ltf = InverseRTransform(ltf)
-    derivs_fun = [inv_ltf.deriv, inv_ltf.deriv2, inv_ltf.deriv3]
+    # Use inverse derivative methods from BaseTransform directly.
+    derivs_fun = [ltf.deriv_inverse, ltf.deriv2_inverse, ltf.deriv3_inverse]
     # Test with 2y + 3y` + 4y``
     coeff = np.array([2, 3, 4])
     coeff_b = _transform_ode_from_derivs(coeff, derivs_fun, x)
diff --git a/src/grid/tests/test_poisson.py b/src/grid/tests/test_poisson.py
index c8e9de749..46fac5b00 100644
--- a/src/grid/tests/test_poisson.py
+++ b/src/grid/tests/test_poisson.py
@@ -31,7 +31,6 @@
 from grid.rtransform import (
     BeckeRTransform,
     IdentityRTransform,
-    InverseRTransform,
     LinearFiniteRTransform,
 )
 from grid.utils import convert_cart_to_sph, generate_real_spherical_harmonics
@@ -221,7 +220,7 @@ def test_poisson_bvp_on_unit_charge_distribution(oned, tf, remove_large_pts, cen
     potential = solve_poisson_bvp(
         molgrids,
         charge_distribution(molgrids.points, centers=centers),
-        InverseRTransform(tf),
+        tf,
         remove_large_pts=remove_large_pts,
         include_origin=True,
     )
@@ -265,7 +264,7 @@ def test_poisson_bvp_gives_the_correct_laplacian(func, centers):
     potential = solve_poisson_bvp(
         molgrids,
         func(molgrids.points, centers),
-        InverseRTransform(btf),
+        btf,
         include_origin=True,
         remove_large_pts=10.0,
     )
@@ -320,7 +319,7 @@ def test_poisson_ivp_gives_the_correct_laplacian(func, centers):
     potential = solve_poisson_ivp(
         molgrids,
         func(molgrids.points, centers),
-        InverseRTransform(btf),
+        btf,
         r_interval=(np.max(radial.points), np.min(radial.points)),
     )
 
@@ -377,7 +376,7 @@ def test_poisson_ivp_on_unit_charge_distribution(centers):
     potential = solve_poisson_ivp(
         molgrids,
         charge_distribution(molgrids.points, centers=centers),
-        InverseRTransform(btf),
+        btf,
         r_interval=(1000.0, 1e-3),
     )
 

From 2a33a3ef4a1520e05a9028ef73195152fb0a1a0c Mon Sep 17 00:00:00 2001
From: aochuba <aochuba52@gmail.com>
Date: Fri, 19 Jun 2026 01:18:22 +0530
Subject: [PATCH 2/4]  Unpack InverseRTransform in poisson solvers for backward
 compatibility

---
 src/grid/poisson.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/grid/poisson.py b/src/grid/poisson.py
index 8cfd95a76..db1cae9b4 100644
--- a/src/grid/poisson.py
+++ b/src/grid/poisson.py
@@ -171,6 +171,12 @@ def solve_poisson_ivp(
     r_interval: tuple = (1000, 1e-5),
     ode_params: dict | type(None) = None,
 ):
+    from grid.rtransform import InverseRTransform
+
+    if isinstance(transform, InverseRTransform):
+        transform = transform._tfm
+
+
     r"""
     Return interpolation of the solution to the Poisson equation solved as an initial value problem.
 
@@ -325,6 +331,11 @@ def solve_poisson_bvp(
     remove_large_pts: float = 1e6,
     ode_params: dict | type(None) = None,
 ):
+    from grid.rtransform import InverseRTransform
+
+    if isinstance(transform, InverseRTransform):
+        transform = transform._tfm
+
     r"""
     Return interpolation of the solution to the Poisson equation solved as a boundary value problem.
 

From a05b1d686c531e50260b1883c6be824c9958af08 Mon Sep 17 00:00:00 2001
From: aochuba <aochuba52@gmail.com>
Date: Mon, 22 Jun 2026 19:00:13 +0530
Subject: [PATCH 3/4]  Move InverseRTransform unwrap block after function
 docstring

---
 src/grid/poisson.py | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/src/grid/poisson.py b/src/grid/poisson.py
index db1cae9b4..32b524c0d 100644
--- a/src/grid/poisson.py
+++ b/src/grid/poisson.py
@@ -171,12 +171,6 @@ def solve_poisson_ivp(
     r_interval: tuple = (1000, 1e-5),
     ode_params: dict | type(None) = None,
 ):
-    from grid.rtransform import InverseRTransform
-
-    if isinstance(transform, InverseRTransform):
-        transform = transform._tfm
-
-
     r"""
     Return interpolation of the solution to the Poisson equation solved as an initial value problem.
 
@@ -219,6 +213,11 @@ def solve_poisson_ivp(
         The solution to Poisson equaiton/potential :math:`g : \mathbb{R}^3 \rightarrow \mathbb{R}`\.
 
     """
+    from grid.rtransform import InverseRTransform
+
+    if isinstance(transform, InverseRTransform):
+        transform = transform._tfm
+
     return _interpolate_molgrid_helper(
         molgrid,
         func_vals,
@@ -331,11 +330,6 @@ def solve_poisson_bvp(
     remove_large_pts: float = 1e6,
     ode_params: dict | type(None) = None,
 ):
-    from grid.rtransform import InverseRTransform
-
-    if isinstance(transform, InverseRTransform):
-        transform = transform._tfm
-
     r"""
     Return interpolation of the solution to the Poisson equation solved as a boundary value problem.
 
@@ -386,6 +380,11 @@ def solve_poisson_bvp(
            polyatomic molecules. The Journal of chemical physics, 89(5), 2993-2997.
 
     """
+    from grid.rtransform import InverseRTransform
+
+    if isinstance(transform, InverseRTransform):
+        transform = transform._tfm
+
     return _interpolate_molgrid_helper(
         molgrid,
         func_vals,

From 8081039052849efecfdd38c5005f9c71afa32dee Mon Sep 17 00:00:00 2001
From: aochuba <aochuba52@gmail.com>
Date: Mon, 22 Jun 2026 19:16:43 +0530
Subject: [PATCH 4/4] misleading error message

---
 src/grid/ode.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/grid/ode.py b/src/grid/ode.py
index 1c055bd8d..adcd8d7ed 100644
--- a/src/grid/ode.py
+++ b/src/grid/ode.py
@@ -123,7 +123,7 @@ def func(x, y):
         if min(x_span) < transform.codomain[0] or max(x_span) > transform.codomain[1]:
             raise ValueError(
                 f"The x_span {min(x_span), max(x_span)} is not within the transform "
-                f"domain {transform.domain}."
+                f"codomain {transform.codomain}."
             )
         # Convert the initial value problem to the new derivative space, only transform up to K-1
         # e.g. the first derivative dV/dx = dV/dr * dr/dx = dV/dr * deriv_inverse