Change default [0, 1] to [-1, 1] for Promolecular

Since the bounds [-1, 1] is used in onedgrids.py,
the default bounds is changed to be consistent.

src/grid/protransform.py

@@ -80,6 +80,8 @@ class CubicProTransform(Grid):
     Promolecular Grid Transformation of a Cubic Grid.
 
     Grid is three dimensional and modeled as Tensor Product of Three, one dimensional grids.
+    Theta space is defined to be :math:`[-1, 1]^3`.
+    Real space is defined to be :math:`\mathbb{R}^3.`
 
     Attributes
     ----------
@@ -100,13 +102,13 @@ class CubicProTransform(Grid):
     integrate(trick=False)
         Integral of a real-valued function over Euclidean space.
     jacobian()
-        Jacobian of the transformation from Real space to Theta/Unit cube space.
+        Jacobian of the transformation from Real space to Theta space :math:`[-1, 1]^3`.
     steepest_ascent_theta()
         Direction of steepest-ascent of a function in theta space from gradient in real space.
     transform():
-        Transform Real point to theta/unit-cube point :math:`[0,1]^3`.
+        Transform Real point to theta point :math:`[-1, 1]^3`.
     inverse(bracket=(-10, 10))
-        Transform theta/unit-cube point to Real space :math:`\mathbb{R}^3`.
+        Transform theta point to Real space :math:`\mathbb{R}^3`.
 
     Examples
     --------
@@ -116,10 +118,10 @@ class CubicProTransform(Grid):
     >> coord = np.array([[0., 0., 0.], [2., 2., 2.]])
 
     Define information of the grid and its weights.
+    >> from grid.onedgrid import GaussChebyshev
+
     >> numb_x = 50
-    >> weights = np.array([(1.0 / (numb_x - 2))] * numb_x)  # Simple Riemannian weights.
-    >> grid = np.linspace(0, 1, num=num_pts, endpoint=True)  # Uniform 1D Grid TODO 0,1
-    >> oned = OneDGrid(grid, weights, domain=(0,1)) TODo 0,1
+    >> oned = GaussChebyshev(numb_x)  # Grid is the same in all x, y, z directions.
     >> promol = CubicProTransform([oned, oned, oned], params.c_m, params.e_m, params.coords)
 
     To integrate some function f.
@@ -138,15 +140,14 @@ class CubicProTransform(Grid):
     TODO: Insert Info About Conditional Distribution Method.
     TODO: Add Infor about how boundarys on theta-space are mapped to np.nan.
 
-
     """
     def __init__(self, oned_grids, coeffs, exps, coords):
         if not isinstance(oned_grids, list):
-            pass
+            raise TypeError("oned_grid should be of type list.")
         if not np.all([isinstance(grid, OneDGrid) for grid in oned_grids]):
-            pass
-        if not np.all([grid.domain == (0., 1.) for grid in oned_grids]):
-            pass
+            raise TypeError("Grid in oned_grids should be of type `OneDGrid`.")
+        if not np.all([grid.domain == (-1, 1.) for grid in oned_grids]):
+            raise ValueError("One Dimensional grid domain should be (-1, 1).")
         if not len(oned_grids) == 3:
             raise ValueError("There should be three One-Dimensional grids in `oned_grids`.")
 
@@ -163,9 +164,13 @@ def __init__(self, oned_grids, coeffs, exps, coords):
         self._prointegral = self._promol.integrate_all()
 
         empty_pts = np.empty((np.prod(self._num_pts), self._dim), dtype=np.float64)
-        weights = np.kron(np.kron(oned_grids[0].weights, oned_grids[1].weights),
-                          oned_grids[2].weights)
-        super().__init__(empty_pts, weights * self._prointegral)
+        weights = np.kron(
+            np.kron(oned_grids[0].weights, oned_grids[1].weights),
+            oned_grids[2].weights
+        )
+        # The prointegral is needed because of promolecular integration.
+        # Divide by 8 needed because the grid is in [-1, 1] rather than [0, 1].
+        super().__init__(empty_pts, weights * self._prointegral / 2.0**self._dim)
         # Transform Cubic Grid in Theta-Space to Real-space.
         self._transform(oned_grids)
 
@@ -192,7 +197,7 @@ def promol(self):
 
     def transform(self, real_pt):
         r"""
-        Transform a real point in three-dimensional Reals to theta/unit cube.
+        Transform a real point in three-dimensional Reals to theta space.
 
         Parameters
         ----------
@@ -202,7 +207,7 @@ def transform(self, real_pt):
         Returns
         -------
         theta_pt : np.ndarray(3)
-            Point in :math:`[0, 1]^3`.
+            Point in :math:`[-1, 1]^3`.
 
         """
         return np.array(
@@ -212,12 +217,12 @@ def transform(self, real_pt):
 
     def inverse(self, theta_pt, bracket=(-10, 10)):
         r"""
-        Transform a theta/unit-cube point to three-dimensional Real space.
+        Transform a theta space point to three-dimensional Real space.
 
         Parameters
         ----------
         theta_pt : np.ndarray(3)
-            Point in :math:`[0, 1]^3`
+            Point in :math:`[-1, 1]^3`
         bracket : (float, float), optional
             Interval where root is suspected to be in Reals.
             Used for "brentq" root-finding method. Default is (-10, 10).
@@ -301,7 +306,7 @@ def integrate(self, *value_arrays, trick=False, tol=1e-10):
 
     def jacobian(self, real_pt):
         r"""
-        Jacobian of the transformation from real space to unit-cube/theta space.
+        Jacobian of the transformation from real space to theta space.
 
         Precisely, it is the lower-triangular matrix
         .. math::
@@ -373,7 +378,7 @@ def jacobian(self, real_pt):
                     )
                     jacobian[i_var, j_deriv] /= transf_den ** 2.0
 
-        return jacobian
+        return 2.0 * jacobian
 
     def derivative(self, real_pt, real_derivative):
         r"""
@@ -401,7 +406,7 @@ def derivative(self, real_pt, real_derivative):
 
     def steepest_ascent_theta(self, real_pt, real_grad):
         r"""
-        Steepest ascent direction of a function in theta/unit-cube space.
+        Steepest ascent direction of a function in theta space.
 
         Steepest ascent is the gradient ie direction of maximum change of a function.
         This guarantees moving in direction of steepest ascent in real-space
@@ -417,34 +422,34 @@ def steepest_ascent_theta(self, real_pt, real_grad):
         Returns
         -------
         theta_grad : np.ndarray(3)
-            Gradient of a function in theta/unit-cube space.
+            Gradient of a function in theta space.
 
         """
         jacobian = self.jacobian(real_pt)
         return jacobian.dot(real_grad)
 
     def _transform(self, oned_grids):
-        # Indices (i, j, k) start from bottom, left-most corner of the unit cube.
+        # Indices (i, j, k) start from bottom, left-most corner of the [-1, 1]^3 cube.
         counter = 0
         for ix in range(self.num_pts[0]):
             cart_pt = [None, None, None]
-            unit_x = oned_grids[0].points[ix]
+            theta_x = oned_grids[0].points[ix]
 
             brack_x = self._get_bracket((ix,), 0)
-            cart_pt[0] = _inverse_coordinate(unit_x, 0, cart_pt, self.promol, brack_x)
+            cart_pt[0] = _inverse_coordinate(theta_x, 0, cart_pt, self.promol, brack_x)
 
             for iy in range(self.num_pts[1]):
-                unit_y = oned_grids[1].points[iy]
+                theta_y = oned_grids[1].points[iy]
 
                 brack_y = self._get_bracket((ix, iy), 1)
-                cart_pt[1] = _inverse_coordinate(unit_y, 1, cart_pt, self.promol, brack_y)
+                cart_pt[1] = _inverse_coordinate(theta_y, 1, cart_pt, self.promol, brack_y)
 
                 for iz in range(self.num_pts[2]):
-                    unit_z = oned_grids[2].points[iz]
+                    theta_z = oned_grids[2].points[iz]
 
                     brack_z = self._get_bracket((ix, iy, iz), 2)
                     cart_pt[2] = _inverse_coordinate(
-                        unit_z, 2, cart_pt, self.promol, brack_z
+                        theta_z, 2, cart_pt, self.promol, brack_z
                     )
 
                     self.points[counter] = cart_pt.copy()
@@ -496,7 +501,7 @@ def _get_bracket(self, indices, i_var):
 
 def _transform_coordinate(real_pt, i_var, promol):
     r"""
-    Transform the `i_var` coordinate of a real point to [0, 1] using promolecular density.
+    Transform the `i_var` coordinate of a real point to [-1, 1] using promolecular density.
 
     Parameters
     ----------
@@ -509,8 +514,8 @@ def _transform_coordinate(real_pt, i_var, promol):
 
     Returns
     -------
-    unit_pt : float
-        The transformed point in [0,1]^D.
+    theta_pt : float
+        The transformed point in :math:`[-1, 1]`.
 
     """
     c_m, e_m, coords, pi_over_exps, dim = astuple(promol)
@@ -535,7 +540,8 @@ def _transform_coordinate(real_pt, i_var, promol):
     transf_den = np.sum(coeff_num)
     transform_value = transf_num / transf_den
 
-    return transform_value
+    # -1. + 2. is needed to transform to [-1, 1], rather than [0, 1].
+    return -1.0 + 2.0 * transform_value
 
 
 def _root_equation(init_guess, prev_trans_pts, theta_pt, i_var, promol):
@@ -549,7 +555,7 @@ def _root_equation(init_guess, prev_trans_pts, theta_pt, i_var, promol):
     prev_trans_pts : list[`i_var` - 1]
         The previous points in real-space that were already transformed.
     theta_pt : float
-        The point in [0, 1] being transformed to the Real space.
+        The point in [-1, 1] being transformed to the Real space.
     i_var : int
         Index of variable being transformed.
     promol : _PromolParams
@@ -569,12 +575,12 @@ def _root_equation(init_guess, prev_trans_pts, theta_pt, i_var, promol):
 
 def _inverse_coordinate(theta_pt, i_var, transformed, promol, bracket=(-10, 10)):
     r"""
-    Transform a point in [0, 1] to the real space corresponding to `i_var` variable.
+    Transform a point in [-1, 1] to the real space corresponding to `i_var` variable.
 
     Parameters
     ----------
     theta_pt : float
-        Point in [0, 1].
+        Point in [-1, 1].
     i_var : int
         Index that is being tranformed. Less than D.
     transformed : list[`i_var` - 1]
@@ -637,7 +643,7 @@ def is_same_sign(x, y):
 
     # Check's if this is a boundary points which is mapped to np.nan
     # These two conditions are added for individual point transformation.
-    if np.abs(theta_pt - 0.0) < 1e-10:
+    if np.abs(theta_pt - -1.0) < 1e-10:
         return np.nan
     if np.abs(theta_pt - 1.0) < 1e-10:
         return np.nan

src/grid/tests/test_protransform.py

@@ -58,11 +58,13 @@ def setUp(self, ss=0.1, return_obj=False):
         params = _PromolParams(c, e, coord, pi_over_exponents=np.sqrt(np.pi / e), dim=3)
         if return_obj:
             num_pts = int(1 / ss) + 1
-            weights = np.array([(1.0 / (num_pts - 2))] * num_pts)
-            oned = OneDGrid(np.linspace(0, 1, num=num_pts, endpoint=True), weights, domain=(0,1))
+            weights = np.array([(2.0 / (num_pts - 2))] * num_pts)
+            oned = OneDGrid(
+                np.linspace(-1, 1, num=num_pts, endpoint=True), weights, domain=(-1, 1),
+            )
 
             obj = CubicProTransform(
-                [oned, oned, oned], params.c_m, params.e_m, params.coords
+                [oned, oned, oned], params.c_m, params.e_m, params.coords,
             )
             return params, obj
         return params
@@ -116,9 +118,10 @@ def formula_transforming_x(x):
                 erf(3.0 ** 0.5 * (x - 2)) + 1.0
             )
 
-            return (first_factor + sec_fac) / (
+            ans = (first_factor + sec_fac) / (
                 5.0 * (np.pi / 2) ** 1.5 + 10.0 * (np.pi / 3.0) ** 1.5
             )
+            return -1.0 + 2.0 * ans
 
         true_ans = _transform_coordinate([pt], 0, self.setUp())
         assert np.abs(true_ans - formula_transforming_x(pt)) < 1e-8
@@ -147,7 +150,7 @@ def formula_transforming_y(x, y):
             dac1 = 5.0 * (np.pi / 2.0) * np.exp(-2.0 * (x - 1.0) ** 2.0)
             dac2 = 10.0 * (np.pi / 3.0) * np.exp(-3.0 * (x - 2.0) ** 2.0)
             den = dac1 + dac2
-            return num / den
+            return -1.0 + 2.0 * (num / den)
 
         true_ans = _transform_coordinate([x, y], 1, self.setUp())
         assert np.abs(true_ans - formula_transforming_y(x, y)) < 1e-8
@@ -180,7 +183,7 @@ def formula_transforming_z(x, y, z):
             )
             den = 5.0 * (np.pi / 2.0) ** 0.5 * np.exp(-2.0 * (a1 ** 2.0 + a2 ** 2.0))
             den += 10.0 * (np.pi / 3.0) ** 0.5 * np.exp(-3.0 * (b1 ** 2.0 + b2 ** 2.0))
-            return (fac1 + fac2) / den
+            return -1.0 + 2.0 * (fac1 + fac2) / den
 
         params, obj = self.setUp(ss=0.5, return_obj=True)
         true_ans = formula_transforming_z(x, y, z)
@@ -207,14 +210,11 @@ def test_transforming_simple_grid(self):
         # Test that converting the point back to unit cube gives [0.5, 0.5, 0.5].
         for i_var in range(0, 3):
             transf = _transform_coordinate(real_pt, i_var, obj.promol)
-            assert np.abs(transf - 0.5) < 1e-5
+            assert np.abs(transf - 0.) < 1e-5
         # Test that all other points are indeed boundary points.
         all_nans = np.delete(obj.points, non_boundary_pt_index, axis=0)
         assert np.all(np.any(np.isnan(all_nans), axis=1))
 
-    # @pytest.mark.parametrize("x", [-2, -2, 0, 2.2])
-    # @pytest.mark.parametrize("y", [-3, 2., -3.2321])
-    # @pytest.mark.parametrize("z", [-2., 1.5, 2.343432])
     def test_transforming_with_inverse_transformation_is_identity(self):
         r"""Test transforming with inverse transformation is identity."""
         # Note that for points far away from the promolecular gets mapped to nan.
@@ -227,7 +227,7 @@ def test_transforming_with_inverse_transformation_is_identity(self):
         assert np.all(np.abs(reverse - pt) < 1e-10)
 
     def test_integrating_itself(self):
-        r"""Test integrating the very same promolecular density."""
+        r"""Test integrating the very same promolecular density used for transformation."""
         params, obj = self.setUp(ss=0.2, return_obj=True)
         promol = []
         for pt in obj.points:
@@ -345,13 +345,13 @@ def setUp(self, ss=0.1, return_obj=False):
         params = _PromolParams(c, e, coord, dim=3, pi_over_exponents=np.sqrt(np.pi / e))
         if return_obj:
             num_pts = int(1 / ss) + 1
-            weights = np.array([(1.0 / (num_pts - 2))] * num_pts)
-            oned_x = OneDGrid(np.linspace(0, 1, num=num_pts, endpoint=True), weights, domain=(0, 1))
-            oned_y = OneDGrid(np.linspace(0, 1, num=num_pts, endpoint=True), weights, domain=(0, 1))
-            oned_z = OneDGrid(np.linspace(0, 1, num=num_pts, endpoint=True), weights, domain=(0, 1))
+            weights = np.array([(2.0 / (num_pts - 2))] * num_pts)
+            oned_x = OneDGrid(
+                np.linspace(-1, 1, num=num_pts, endpoint=True), weights, domain=(-1, 1)
+            )
 
             obj = CubicProTransform(
-                [oned_x, oned_y, oned_z], params.c_m, params.e_m, params.coords
+                [oned_x, oned_x, oned_x], params.c_m, params.e_m, params.coords
             )
             return params, obj
         return params
@@ -367,12 +367,12 @@ def promolecular_in_x(grid, every_grid):
             return promol_x, promol_x_all
 
         true_ans = _transform_coordinate([pt], 0, self.setUp())
-        grid = np.arange(-10.0, pt, 0.00001)  # Integration till a x point
-        every_grid = np.arange(-10.0, 10.0, 0.00001)  # Full Integration
+        grid = np.arange(-4.0, pt, 0.000005)  # Integration till a x point
+        every_grid = np.arange(-5.0, 10.0, 0.00001)  # Full Integration
         promol_x, promol_x_all = promolecular_in_x(grid, every_grid)
 
         # Integration over y and z cancel out from numerator and denominator.
-        actual = np.trapz(promol_x, grid) / np.trapz(promol_x_all, every_grid)
+        actual = -1.0 + 2.0 * np.trapz(promol_x, grid) / np.trapz(promol_x_all, every_grid)
         assert np.abs(true_ans - actual) < 1e-5
 
     @pytest.mark.parametrize("x", [-10.0, -2.0, 0.0, 2.2, 1.23])
@@ -387,13 +387,13 @@ def promolecular_in_y(grid, every_grid):
             return promol_y_all, promol_y
 
         true_ans = _transform_coordinate([x, y], 1, self.setUp())
-        grid = np.arange(-10.0, y, 0.00001)  # Integration till a x point
-        every_grid = np.arange(-10.0, 10.0, 0.00001)  # Full Integration
+        grid = np.arange(-5.0, y, 0.000001)  # Integration till a x point
+        every_grid = np.arange(-5.0, 10.0, 0.00001)  # Full Integration
         promol_y_all, promol_y = promolecular_in_y(grid, every_grid)
 
         # Integration over z cancel out from numerator and denominator.
         # Further, gaussian at a point does too.
-        actual = np.trapz(promol_y, grid) / np.trapz(promol_y_all, every_grid)
+        actual = -1.0 + 2.0 * np.trapz(promol_y, grid) / np.trapz(promol_y_all, every_grid)
         assert np.abs(true_ans - actual) < 1e-5
 
     @pytest.mark.parametrize("x", [-10.0, -2.0, 0.0, 2.2])
@@ -408,11 +408,11 @@ def promolecular_in_z(grid, every_grid):
             promol_z_all = 5.0 * np.exp(-2.0 * (every_grid - 3.0) ** 2.0)
             return promol_z_all, promol_z
 
-        grid = np.arange(-10.0, z, 0.00001)  # Integration till a x point
-        every_grid = np.arange(-10.0, 10.0, 0.00001)  # Full Integration
+        grid = np.arange(-5.0, z, 0.00001)  # Integration till a x point
+        every_grid = np.arange(-5.0, 10.0, 0.00001)  # Full Integration
         promol_z_all, promol_z = promolecular_in_z(grid, every_grid)
 
-        actual = np.trapz(promol_z, grid) / np.trapz(promol_z_all, every_grid)
+        actual = -1.0 + 2.0 * np.trapz(promol_z, grid) / np.trapz(promol_z_all, every_grid)
         true_ans = _transform_coordinate([x, y, z], 2, self.setUp())
         assert np.abs(true_ans - actual) < 1e-5