diff --git a/pyproject.toml b/pyproject.toml index 9ba1edde..84df5147 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ dependencies = [ "pytest>=8.0.0", "scipy>=1.4", "importlib_resources", + "dataclasses", "sympy", ] dynamic = ["version"] diff --git a/src/grid/cubic.py b/src/grid/cubic.py index 1c2998f6..d5f77ccf 100644 --- a/src/grid/cubic.py +++ b/src/grid/cubic.py @@ -101,8 +101,10 @@ def get_points_along_axes(self): x = self.points[coords_x, 0] return x, y - def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, method="cubic"): - r"""Interpolate function and its derivatives on cubic grid. + def interpolate( + self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, method="cubic", grid_pts=None + ): + r"""Interpolate function value at a given point. Only implemented in three-dimensions. @@ -130,10 +132,14 @@ def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, met If zero, then the function in z-direction is interpolated. If greater than zero, then the "nu_z"th-order derivative in the z-direction is interpolated. - method : str, optional + method: str, optional The method of interpolation to perform. Supported are "cubic" (most accurate but computationally expensive), "linear", or "nearest" (least accurate but cheap computationally). The last two methods use SciPy's RegularGridInterpolator function. + grid_pts: list[OneDGrids], optional + If provided, then uses `grid_pts` rather than the points of the HyperRectangle class + `self.points` to construct interpolation. Useful when doing a promolecular + transformation. Returns ------- @@ -141,6 +147,18 @@ def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, met The interpolation of a function (or of it's derivatives) at a :math:`M` point. """ + # Needed because CubicProTransform is a subclass of this method and has its own + # interpolate function. Since interpolate references itself, it chooses + # CubicProTransform rather than _HyperRectangleGrid class. + return self._interpolate( + points, values, use_log, nu_x, nu_y, nu_z, method, grid_pts + ) + + def _interpolate( + self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, method="cubic", grid_pts=None + ): + r"""Core of the Interpolate Algorithm.""" + if method not in ["cubic", "linear", "nearest"]: raise ValueError( f"Argument method should be either cubic, linear, or nearest , got {method}" @@ -154,6 +172,10 @@ def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, met f"Number of function values {values.shape[0]} does not match number of " f"grid points {np.prod(self.shape)}." ) + if grid_pts is not None and not isinstance(grid_pts, np.ndarray): + raise TypeError( + f"The grid points {type(grid_pts)} should have type None or numpy array." + ) if use_log: values = np.log(values) @@ -165,6 +187,10 @@ def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, met interpolate = RegularGridInterpolator((x, y, z), values, method=method) return interpolate(points) + # If grid_pts isn't specified, then use the grid stored as the class attribute. + if grid_pts is None: + grid_pts = self.points + # Interpolate the Z-Axis. def z_spline(z, x_index, y_index, nu_z=nu_z): # x_index, y_index is assumed to be in the grid while z is not assumed. @@ -173,7 +199,7 @@ def z_spline(z, x_index, y_index, nu_z=nu_z): small_index = self.coordinates_to_index((x_index, y_index, 1)) large_index = self.coordinates_to_index((x_index, y_index, self.shape[2] - 2)) val = CubicSpline( - self.points[small_index:large_index, 2], + grid_pts[small_index:large_index, 2], values[small_index:large_index], )(z, nu_z) return val @@ -183,7 +209,7 @@ def y_splines(y, x_index, z, nu_y=nu_y): # The `1` and `self.num_puts[1] - 2` is needed because I don't want the boundary. # Assumes x_index is in the grid while y, z may not be. val = CubicSpline( - self.points[np.arange(1, self.shape[1] - 2) * self.shape[2], 1], + grid_pts[np.arange(1, self.shape[1] - 2) * self.shape[2], 1], [z_spline(z, x_index, y_index, nu_z) for y_index in range(1, self.shape[1] - 2)], )(y, nu_y) # Trying to vectorize over z-axis and y-axis, this computes the interpolation for every @@ -195,7 +221,7 @@ def y_splines(y, x_index, z, nu_y=nu_y): # Interpolate the point (x, y, z) from a list of interpolated points on x,y-axis. def x_spline(x, y, z, nu_x): val = CubicSpline( - self.points[np.arange(1, self.shape[0] - 2) * self.shape[1] * self.shape[2], 0], + grid_pts[np.arange(1, self.shape[0] - 2) * self.shape[1] * self.shape[2], 0], [y_splines(y, x_index, z, nu_y) for x_index in range(1, self.shape[0] - 2)], )(x, nu_x) # Trying to vectorize over x-axis, this computes the interpolation for every @@ -207,7 +233,9 @@ def x_spline(x, y, z, nu_x): if use_log: # All derivatives require the interpolation of f at (x,y,z) interpolated = np.exp( - self.interpolate(points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0) + self._interpolate( + points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, grid_pts=grid_pts + ) ) # Only consider taking the derivative in only one direction one_var_deriv = sum([nu_x == 0, nu_y == 0, nu_z == 0]) == 2 @@ -218,21 +246,28 @@ def x_spline(x, y, z, nu_x): elif one_var_deriv: # Taking the k-th derivative wrt to only one variable (x, y, z) # Interpolate d^k ln(f) d"deriv_var" for all k from 1 to "deriv_var" + # Each entry of `derivs` is the interpolation of the derivative eval on points. if nu_x > 0: derivs = [ - self.interpolate(points, values, use_log=False, nu_x=i, nu_y=0, nu_z=0) + self._interpolate( + points, values, use_log=False, nu_x=i, nu_y=0, nu_z=0, grid_pts=grid_pts + ) for i in range(1, nu_x + 1) ] deriv_var = nu_x elif nu_y > 0: derivs = [ - self.interpolate(points, values, use_log=False, nu_x=0, nu_y=i, nu_z=0) + self._interpolate( + points, values, use_log=False, nu_x=0, nu_y=i, nu_z=0, grid_pts=grid_pts + ) for i in range(1, nu_y + 1) ] deriv_var = nu_y else: derivs = [ - self.interpolate(points, values, use_log=False, nu_x=0, nu_y=0, nu_z=i) + self._interpolate( + points, values, use_log=False, nu_x=0, nu_y=0, nu_z=i, grid_pts=grid_pts + ) for i in range(1, nu_z + 1) ] deriv_var = nu_z diff --git a/src/grid/protransform.py b/src/grid/protransform.py new file mode 100644 index 00000000..6ac11f8f --- /dev/null +++ b/src/grid/protransform.py @@ -0,0 +1,1031 @@ +# GRID is a numerical integration module for quantum chemistry. +# +# Copyright (C) 2011-2019 The GRID Development Team +# +# This file is part of GRID. +# +# GRID is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# GRID is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# -- +r"""Promolecular Grid Transformation.""" + +from dataclasses import astuple, dataclass, field + +from grid.basegrid import Grid, OneDGrid +from grid.cubic import _HyperRectangleGrid + +import numpy as np + +from scipy.interpolate import CubicSpline +from scipy.linalg import solve_triangular +from scipy.optimize import root_scalar +from scipy.special import erf + +__all__ = ["CubicProTransform"] + + +class CubicProTransform(_HyperRectangleGrid): + r""" + Promolecular Grid Transformation of a Cubic Grid in :math:`[-1, 1]^3`. + + 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 + ---------- + shape : (int, int, int) + The number of points, including both of the end/boundary points, in x, y, and z direction. + prointegral : float + The integration value of the promolecular density over Euclidean space. + promol : namedTuple + Data about the Promolecular density. + points : np.ndarray(N, 3) + The transformed points in real space. + weights : np.ndarray(N,) + The weights multiplied by `prointegral`. + + Methods + ------- + integrate(trick=False) + Integral of a real-valued function over Euclidean space. Can use promolecular trick. + jacobian() + Jacobian of the transformation from Real space to Theta space :math:`[-1, 1]^3`. + hessian() + Hessian 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 space :math:`[-1, 1]^3`. + inverse(bracket=(-10, 10)) + Transform Theta point to Real space :math:`\mathbb{R}^3`. + interpolate(use_log=False, nu=0) + Interpolate a function (or its logarithm) at a real point. Can interpolate its derivative. + + Examples + -------- + Define information of the Promolecular Density. + >> c = np.array([[5.], [10.]]) + >> e = np.array([[2.], [3.]]) + >> 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 + This is a grid in :math:`[-1, 1]`. + >> oned = GaussChebyshev(numb_x) + One dimensional 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. + >> def f(pt): + >> return np.exp(-0.1 * np.linalg.norm(pt, axis=1)**2.) + >> func_values = f(promol.points) + >> print("The integral is %.4f" % promol.integrate(func_values, trick=False) + + References + ---------- + .. [1] J. I. Rodríguez, D. C. Thompson, P. W. Ayers, and A. M. Koster, "Numerical integration + of exchange-correlation energies and potentials using transformed sparse grids." + + Notes + ----- + Let :math:`\rho^o(x, y, z) = \sum_{i=1}^M \sum_{j=1}^D e^{}` be the Promolecular density of a \ + linear combination of Gaussian functions. + + The conditional distribution transformation from :math:`\mathbb{R}^3` to :math:`[-1, 1]^3` + transfers the (x, y, z) coordinates in :math:`\mathbb{R}^3` to a set of coordinates, + denoted as :math:`(\theta_x, \theta_y, \theta_z)`, in :math:`[-1,1]^3` that are "bunched" + up where :math:`\rho^o` is large. + + Precisely it is, + + .. math:: + \begin{eqnarray} + \theta_x(x) :&= + -1 + 2 \frac{\int_{-\infty}^x \int \int \rho^o(x, y, z)dx dy dz } + {\int \int \int \rho^o(x, y, z)dxdydz}\\ + \theta_y(x, y) :&= + -1 + 2 \frac{\int_{-\infty}^y \int \rho^o(x, y, z)dy dz } + {\int \int \rho^o(x, y, z)dydz} \\ + \theta_z(x, y, z) :&= + -1 + 2 \frac{\int_{-\infty}^z \rho^o(x, y, z)dz } + {\int \rho^o(x, y, z)dz}\\ + \end{eqnarray} + + Integration of a integrable function :math:`f : \mathbb{R}^3 \rightarrow \mathbb{R}` can be + done as follows in theta space: + + .. math:: + \int \int \int f(x, y, z)dxdy dz \approx + \frac{1}{8} N \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 \frac{f(\theta_x, \theta_y, \theta_z)} + {\rho^o(\theta_x, \theta_y, \theta_z)} d\theta_x d\theta_y d\theta_z, + + \text{where } N = \int \int \int \rho^o(x, y, z) dx dy dz. + + Note that this class always assumed the boundary of [-1, 1]^3 is always included. + + """ + + def __init__(self, oned_grids, coeffs, exps, coords, cut_off=1e-8): + r""" + Construct CubicProTransform object. + + Parameters + ---------- + oned_grids: List[OneDGrid] + List of three one-dimensional grid representing the grids along x-axis. + coeffs: List[List[float]] + Coefficients of the promolecular transformation over :math:`M` centers. + exps: List[List[float]] + Exponents of the promolecular transformation over :math:`M` centers. + coords: ndarray(M, 3) + The coordinates of the promolecular expansion. + cut_off: float + If the distance between a point in theta-space to the boundary is less than the + cut_off, then the point is considered to be part of the boundary. + + """ + if not isinstance(oned_grids, list): + raise TypeError("oned_grid should be of type list.") + if not np.all([isinstance(grid, OneDGrid) for grid in oned_grids]): + 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`.") + self._l_bnd = -1.0 + self._u_bnd = 1.0 + self._shape = tuple([grid.size for grid in oned_grids]) + dimension = len(oned_grids) + + # pad coefficients and exponents with zeros to have the same size, easier to use numpy. + coeffs, exps = _pad_coeffs_exps_with_zeros(coeffs, exps) + self._promol = _PromolParams(coeffs, exps, coords, dimension) + self._prointegral = self._promol.integrate_all() + + weights = np.kron( + np.kron(oned_grids[0].weights, oned_grids[1].weights), oned_grids[2].weights + ) + # Transform Cubic Grid in Theta-Space to Real-space. + points = self._transform(oned_grids, cut_off) + # 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__(points, weights * self._prointegral / 2.0**dimension, self._shape) + + @property + def l_bnd(self): + r"""float: Lower bound in theta-space. Any point in theta-space that contains this point + gets mapped to infinity.""" + return self._l_bnd + + @property + def u_bnd(self): + r"""float: Upper bound in theta-space. Any point in theta-space that contains this point + gets mapped to infinity.""" + return self._u_bnd + + @property + def prointegral(self): + r"""Return integration of Promolecular density.""" + return self._prointegral + + @property + def promol(self): + r"""Return `PromolParams` data class.""" + return self._promol + + def transform(self, real_pt): + r""" + Transform a real point in three-dimensional Reals to theta space. + + Parameters + ---------- + real_pt : np.ndarray(3) + Point in :math:`\mathbb{R}^3` + + Returns + ------- + theta_pt : np.ndarray(3) + Point in :math:`[-1, 1]^3`. + + """ + return np.array( + [_transform_coordinate(real_pt, i, self.promol) for i in range(0, self.promol.dim)] + ) + + def inverse(self, theta_pt, bracket=(-10, 10)): + r""" + Transform a theta space point to three-dimensional Real space. + + Parameters + ---------- + theta_pt : np.ndarray(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). + + Returns + ------- + real_pt : np.ndarray(3) + Point in :math:`\mathbb{R}^3` + + Notes + ----- + - If a point is far away from the promolecular density, then it will be mapped + to `np.nan`. + + """ + real_pt = [] + for i in range(0, self.promol.dim): + scalar = _inverse_coordinate(theta_pt[i], i, real_pt[:i], self.promol, bracket) + real_pt.append(scalar) + return np.array(real_pt) + + def integrate(self, *value_arrays, trick=False, tol=1e-10): + r""" + Integrate any real-valued function on Euclidean space. + + Assumes the function decays faster than the promolecular density. + + Parameters + ---------- + *value_arrays : (np.ndarray(N, dtype=float),) + One or multiple value array to integrate. + trick : bool, optional + If true, uses the promolecular trick. + tol : float, optional + Integrand is set to zero whenever promolecular density is less than tolerance. + Default value is 1e-10. + + Returns + ------- + float : + Return the integration of the function. + + Raises + ------ + TypeError + Input integrand is not of type np.ndarray. + ValueError + Input integrand array is given or not of proper shape. + + Notes + ----- + - Formula for the integration of a integrable function + :math:`f : \mathbb{R}^3 \rightarrow \mathbb{R}` is done as follows: + + .. math:: + \int \int \int f(x, y, z)dxdy dz \approx + \frac{1}{8} N \int_{-1}^1 \int_{-1}^1 \int_{-1}^1 \frac{f(\theta_x, \theta_y, \theta_z)} + {\rho^o(\theta_x, \theta_y, \theta_z)} d\theta_x d\theta_y d\theta_z, + + \text{where } N = \int \int \int \rho^o(x, y, z) dx dy dz. + + - This method assumes function f decays faster than the promolecular density. + + """ + promolecular = self.promol.promolecular(self.points) + # Integrand is set to zero when promolecular is less than certain value and, + # When on the boundary (hence when promolecular is nan). + cond = (promolecular <= tol) | (np.isinf(promolecular)) + promolecular = np.ma.masked_where(cond, promolecular, copy=False) + + integrands = [] + for arr in value_arrays: + # This is needed as it gives incorrect results when arr.dtype isn't object. + assert arr.dtype != object, "Array dtype should not be object." + # This may be refactored to fit in the general promolecular trick in `grid`. + # Masked array is needed since division by promolecular contains nan. + if trick: + integrand = np.ma.divide(arr - promolecular, promolecular) + else: + integrand = np.ma.divide(arr, promolecular) + # Function/Integrand evaluated at points on the boundary is set to zero. + np.ma.fix_invalid(integrand, copy=False, fill_value=0) + integrands.append(integrand) + + if trick: + return self.prointegral + super().integrate(*integrands) + return super().integrate(*integrands) + + def derivative(self, real_pt, real_derivative): + r""" + Directional derivative in theta space. + + Parameters + ---------- + real_pt : np.ndarray(3) + Point in :math:`\mathbb{R}^3`. + real_derivative : np.ndarray(3) + Derivative of a function in real space with respect to x, y, z coordinates. + + Returns + ------- + theta_derivative : np.ndarray(3) + Derivative of a function in theta space with respect to theta coordinates. + + Notes + ----- + This does not preserve the direction of steepest-ascent/gradient. + + See Also + -------- + steepest_ascent_theta : Steepest-ascent direction. + + """ + jacobian = self.jacobian(real_pt) + return solve_triangular(jacobian.T, real_derivative) + + def steepest_ascent_theta(self, real_pt, real_grad): + r""" + 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 + corresponds to moving in the direction of the gradient in theta-space. + + Parameters + ---------- + real_pt : np.ndarray(3) + Point in :math:`\mathbb{R}^3` + real_grad : np.ndarray(3) + Gradient of a function in real space. + + Returns + ------- + theta_grad : np.ndarray(3) + Gradient of a function in theta space. + + """ + jacobian = self.jacobian(real_pt) + return jacobian.dot(real_grad) + + def interpolate( + self, + points, + values, + oned_grids, + use_log=False, + nu=0, + ): + r""" + Interpolate function at a point. + + Parameters + ---------- + real_pt : np.ndarray(3,) + Point in :math:`\mathbb{R}^3` that needs to be interpolated. + values : np.ndarray(N,) + Function values at each point of the grid `points`. + oned_grids = list(OneDGrids) + List containing three one-dimensional grid corresponding to x, y, z direction + use_log : bool + If true, then logarithm is applied before interpolating to the function values, + including `func_values`. + nu : int + If zero, then the function is interpolated. + If one, then the derivative is interpolated. + + Returns + ------- + float : + If nu is 0: Returns the interpolated of a function at a real point. + If nu is 1: Returns the interpolated derivative of a function at a real point. + + """ + if nu not in (0, 1): + raise ValueError("The parameter nu %d is either zero or one " % nu) + + grid_pts = ( + np.vstack( + np.meshgrid( + oned_grids[0].points, + oned_grids[1].points, + oned_grids[2].points, + indexing="ij", + ) + ) + .reshape(3, -1) + .T + ) + theta_points = np.array([self.transform(x) for x in points], dtype=float) + + if nu == 0: + interpolation = super()._interpolate( + theta_points, values, use_log, 0, 0, 0, grid_pts=grid_pts + ) + if nu == 1: + # Interpolate in the theta-space and take the derivatives + interpolate_x = super()._interpolate( + theta_points, values, use_log, 1, 0, 0, grid_pts=grid_pts + ) + interpolate_y = super()._interpolate( + theta_points, values, use_log, 0, 1, 0, grid_pts=grid_pts + ) + interpolate_z = super()._interpolate( + theta_points, values, use_log, 0, 0, 1, grid_pts=grid_pts + ) + interpolation = np.vstack((interpolate_x.T, interpolate_y.T, interpolate_z.T)).T + # Transform the derivatives back to real-space. + interpolation = np.array( + [self.jacobian(points[i]).dot(interpolation[i]) for i in range(len(interpolation))] + ) + + return interpolation + + def jacobian(self, real_pt): + r""" + Jacobian of the transformation from real space to theta space. + + Precisely, it is the lower-triangular matrix + + .. math:: + \begin{bmatrix} + \frac{\partial \theta_x}{\partial X} & 0 & 0 \\ + \frac{\partial \theta_y}{\partial X} & \frac{\partial \theta_y}{\partial Y} & 0 \\ + \frac{\partial \theta_z}{\partial X} & \frac{\partial \theta_z}{\partial Y} & + \frac{\partial \theta_z}{\partial Z} + \end{bmatrix}. + + Parameters + ---------- + real_pt : np.ndarray(3,) + Point in :math:`\mathbb{R}^3`. + + Returns + ------- + np.ndarray(3, 3) : + Jacobian of transformation. + + """ + jacobian = np.zeros((3, 3), dtype=np.float64) + + c_m, e_m, coords, pi_over_exps, dim = astuple(self.promol) + + # Distance to centers/nuclei`s and Prefactors. + diff_coords = real_pt - coords + diff_squared = diff_coords**2.0 + # If i_var is zero, then distance is just all zeros. + for i_var in range(0, self.promol.dim): + + # Basic-Level arrays for integration and derivatives. + ( + distance, + single_gauss, + integrate_till_pt_x, + transf_num, + transf_den, + ) = self.promol.helper_for_derivatives(diff_squared, diff_coords, i_var) + + for j_deriv in range(0, i_var + 1): + if i_var == j_deriv: + # Derivative eliminates `integrate_till_pt_x`, and adds a Gaussian. + inner_term = single_gauss * np.exp(-e_m * diff_squared[:, i_var][:, np.newaxis]) + jacobian[i_var, i_var] = np.sum(inner_term) / transf_den + elif j_deriv < i_var: + # Derivative of inside of Gaussian. + deriv_inside = self.promol.derivative_gaussian(diff_coords, j_deriv) + deriv_num = np.sum(single_gauss * integrate_till_pt_x * deriv_inside) + deriv_den = np.sum(single_gauss * deriv_inside * pi_over_exps) + # Quotient Rule + jacobian[i_var, j_deriv] = deriv_num * transf_den - transf_num * deriv_den + jacobian[i_var, j_deriv] /= transf_den**2.0 + + return 2.0 * jacobian + + def hessian(self, real_pt): + r""" + Hessian of the transformation. + + The Hessian :math:`H` is a three-dimensional array with (i, j, k)th entry: + + .. math:: + H_{i, j, k} = \frac{\partial^2 \theta_i(x_0, \cdots, x_{i-1}}{\partial x_i \partial x_j} + + \text{where } (x_0, x_1, x_2) := (x, y, z). + + Parameters + ---------- + real_pt : np.ndarray(3,) + Real point in :math:`\mathbb{R}^3`. + + Returns + ------- + hessian : np.ndarray(3, 3, 3) + The (i, j, k)th entry is the partial derivative of the ith transformation function + with respect to the jth, kth coordinate. e.g. when i = 0, then hessian entry at + (i, j, k) is zero unless j = k = 0. + + """ + hessian = np.zeros((self.ndim, self.ndim, self.ndim), dtype=np.float64) + + c_m, e_m, coords, pi_over_exps, dim = astuple(self.promol) + + # Distance to centers/nuclei`s and Prefactors. + diff_coords = real_pt - coords + diff_squared = diff_coords**2.0 + + # i_var is the transformation to theta-space. + # j_deriv is the first partial derivative wrt x, y, z. + # k_deriv is the second partial derivative wrt x, y, z. + for i_var in range(0, 3): + + # Basic-Level arrays for integration and derivatives. + ( + distance, + single_gauss, + integrate_till_pt_x, + transf_num, + transf_den, + ) = self.promol.helper_for_derivatives(diff_squared, diff_coords, i_var) + + for j_deriv in range(0, i_var + 1): + for k_deriv in range(0, i_var + 1): + # num is the numerator of transformation function. + # den is the denominator of transformation function. + # dnum_dk is the derivative of numerator wrt to k_deriv. + # dnum_dkdj is the derivative of num wrt to j_deriv then k_deriv. + # The derivative will store the result and pass it to the Hessian. + derivative = 0.0 + + if i_var == j_deriv: + gauss_extra = single_gauss * np.exp( + -e_m * diff_squared[:, j_deriv][:, np.newaxis] + ) + if j_deriv == k_deriv: + # Diagonal derivative e.g. d(theta_X)(dx dx) + gauss_extra *= self.promol.derivative_gaussian(diff_coords, j_deriv) + derivative = np.sum(gauss_extra) / transf_den + else: + # Partial derivative of diagonal derivative e.g. d^2(theta_y)(dy dx). + deriv_inside = self.promol.derivative_gaussian(diff_coords, k_deriv) + dnum_dkdj = np.sum(gauss_extra * deriv_inside) + dden_dk = np.sum(single_gauss * deriv_inside * pi_over_exps) + # Numerator is different from `transf_num` since Gaussian is added. + dnum_dj = np.sum(gauss_extra) + # Quotient Rule + derivative = dnum_dkdj * transf_den - dnum_dj * dden_dk + derivative /= transf_den**2.0 + + # Here, quotient rule all the way down. + elif j_deriv < i_var: + if k_deriv == i_var: + gauss_extra = single_gauss * np.exp( + -e_m * diff_squared[:, k_deriv][:, np.newaxis] + ) + deriv_inside = self.promol.derivative_gaussian(diff_coords, j_deriv) + ddnum_djdi = np.sum(gauss_extra * deriv_inside) + dden_dj = np.sum(single_gauss * deriv_inside * pi_over_exps) + # Quotient Rule + dnum_dj = np.sum(gauss_extra) + derivative = ddnum_djdi * transf_den - dnum_dj * dden_dj + derivative /= transf_den**2.0 + + elif k_deriv == j_deriv: + # Double Quotient Rule. + # See wikipedia "Quotient Rules Higher order formulas". + deriv_inside = self.promol.derivative_gaussian(diff_coords, k_deriv) + dnum_dj = np.sum(single_gauss * integrate_till_pt_x * deriv_inside) + dden_dj = np.sum(single_gauss * pi_over_exps * deriv_inside) + + prod_rule = deriv_inside**2.0 - 2.0 * e_m + sec_deriv_num = np.sum(single_gauss * integrate_till_pt_x * prod_rule) + sec_deriv_den = np.sum(single_gauss * pi_over_exps * prod_rule) + + output = sec_deriv_num * transf_den - dnum_dj * dden_dj + output /= transf_den**2.0 + quot = transf_den * (dnum_dj * dden_dj + transf_num * sec_deriv_den) + quot -= 2.0 * transf_num * dden_dj * dden_dj + derivative = output - quot / transf_den**3.0 + + elif k_deriv != j_deriv: + # K is i_Sec_diff and i is i_diff + deriv_inside = self.promol.derivative_gaussian(diff_coords, j_deriv) + deriv_inside_sec = self.promol.derivative_gaussian(diff_coords, k_deriv) + gauss_and_inte_x = single_gauss * integrate_till_pt_x + gauss_and_inte = single_gauss * pi_over_exps + + dnum_di = np.sum(gauss_and_inte_x * deriv_inside) + dden_di = np.sum(gauss_and_inte * deriv_inside) + + dnum_dk = np.sum(gauss_and_inte_x * deriv_inside_sec) + dden_dk = np.sum(gauss_and_inte * deriv_inside_sec) + + ddnum_dkdk = np.sum(gauss_and_inte_x * deriv_inside * deriv_inside_sec) + ddden_dkdk = np.sum(gauss_and_inte * deriv_inside * deriv_inside_sec) + + output = ddnum_dkdk / transf_den + output -= dnum_di * dden_dk / transf_den**2.0 + product = dnum_dk * dden_di + transf_num * ddden_dkdk + derivative = output + derivative -= product / transf_den**2.0 + derivative += 2.0 * transf_num * dden_di * dden_dk / transf_den**3.0 + + # The 2.0 is needed because we're in [-1, 1] rather than [0, 1]. + hessian[i_var, j_deriv, k_deriv] = 2.0 * derivative + + return hessian + + def _transform(self, oned_grids, cut_off=1e-8): + # Transform the entire grid. + # Indices (i, j, k) start from bottom, left-most corner of the [-1, 1]^3 cube. + def _is_boundary_pt(theta_pt): + if np.abs(theta_pt - self.l_bnd) < cut_off or np.abs(theta_pt - self.u_bnd) < cut_off: + return True + return False + + counter = 0 + points = np.empty((np.prod(self.shape), len(oned_grids)), dtype=np.float64) + for ix in range(self.shape[0]): + cart_pt = [None, None, None] + theta_x = oned_grids[0].points[ix] + is_boundary = _is_boundary_pt(theta_x) + brack_x = self._get_bracket((ix,), 0, points, is_boundary) + cart_pt[0] = _inverse_coordinate(theta_x, 0, cart_pt, self.promol, brack_x) + + for iy in range(self.shape[1]): + theta_y = oned_grids[1].points[iy] + is_boundary = _is_boundary_pt(theta_y) + brack_y = self._get_bracket((ix, iy), 1, points, is_boundary) + cart_pt[1] = _inverse_coordinate(theta_y, 1, cart_pt, self.promol, brack_y) + + for iz in range(self.shape[2]): + theta_z = oned_grids[2].points[iz] + is_boundary = _is_boundary_pt(theta_z) + brack_z = self._get_bracket((ix, iy, iz), 2, points, is_boundary) + cart_pt[2] = _inverse_coordinate(theta_z, 2, cart_pt, self.promol, brack_z) + points[counter] = cart_pt.copy() + counter += 1 + return points + + def _get_bracket(self, indices, i_var, prev_points, is_boundary): + r""" + Obtain brackets for root-finder based on the coordinate of the point. + + Parameters + ---------- + indices : tuple(int, int, int) + The indices of a point, where (0, 0, 0) is the bottom, left-most, down point + of the cube. + i_var : int + Index of point being transformed. + prev_points: ndarray(M, 3) + Points that are transformed and empty points that haven't been transformed yet. + The points that are transformed is used to obtain brackets for this current + point that hasn't been transformed yet. + is_boundary: bool + If the current indices is a boundary point, then returns inf, as it maps to infinity. + + Returns + ------- + (float, float) : + The bracket for the root-finder solver. + + """ + # If it is a boundary point, then return nan. Done by indices. + if is_boundary: + return np.inf, np.inf + # If it is a new point, with no nearby point, get a large initial guess. + else: + min = np.min(self.promol.coords[:, i_var]) - 20.0 + max = np.max(self.promol.coords[:, i_var]) + 20.0 + return min, max + # If the previous point has been converted, use that as a initial guess. + # if i_var == 0: + # print(" shape = %s" % str(self.shape)) + # index = (indices[0] - 1) * self.shape[1] * self.shape[2] + # print(" index = %d" % index) + # elif i_var == 1: + # index = indices[0] * self.shape[1] * self.shape[2] + self.shape[2] * (indices[1] - 1) + # elif i_var == 2: + # index = ( + # indices[0] * self.shape[1] * self.shape[2] + # + self.shape[2] * indices[1] + # + indices[2] + # - 1 + # ) + # print(" prev_points = %s" % str(prev_points)) + # print(" shape prev_points = %s" % str(prev_points.shape)) + # print(f" returned {prev_points[index, i_var]}") + + # return prev_points[index, i_var], prev_points[index, i_var] + 10.0 + + +@dataclass +class _PromolParams: + r""" + Private class for Promolecular Density information. + + Contains helper-functions for Promolecular Transformation. + They are coded as pipe-lines for this special purpose and + the reason why "diff_coords" is chosen as a attribute rather + than a generic "[x, y, z]" point. + + """ + + c_m: np.ndarray # Coefficients of Promolecular. + e_m: np.ndarray # Exponents of Promolecular. + coords: np.ndarray # Centers/Coordinates of Each Gaussian. + pi_over_exponents: np.ndarray = field(init=False) + dim: int = 3 + + def __post_init__(self): + r"""Initialize pi_over_exponents.""" + # Rather than computing this repeatedly. It is fixed. + with np.errstate(divide="ignore"): + self.pi_over_exponents = np.sqrt(np.pi / self.e_m) + self.pi_over_exponents[self.e_m == 0.0] = 0.0 + + def integrate_all(self): + r"""Integration of Gaussian over Entire Real space ie :math:`\mathbb{R}^D`.""" + return np.sum(self.c_m * self.pi_over_exponents**self.dim) + + def derivative_gaussian(self, diff_coords, j_deriv): + r"""Return derivative of single Gaussian but without exponential.""" + return -self.e_m * 2.0 * diff_coords[:, j_deriv][:, np.newaxis] + + def integration_gaussian_till_point(self, diff_coords, i_var, with_factor=False): + r"""Integration of Gaussian wrt to `i_var` variable till a point (inside diff_coords).""" + coord_ivar = diff_coords[:, i_var][:, np.newaxis] + integration = (erf(np.sqrt(self.e_m) * coord_ivar) + 1.0) / 2.0 + if with_factor: + # Included the (pi / exponents), this becomes the actual integral. + # Not including the (pi / exponents) increasing computation slightly faster. + return integration * self.pi_over_exponents + return integration + + def single_gaussians(self, distance): + r"""Return matrix with entries a single gaussian evaluated at the float distance.""" + return self.c_m * np.exp(-self.e_m * distance) + + def promolecular(self, points): + r""" + Evaluate the promolecular density over a grid. + + Parameters + ---------- + points : np.ndarray(N, D) + Points in :math:`\mathbb{R}^D`. + + Returns + ------- + np.ndarray(N,) : + Promolecular density evaluated at the points. + + """ + # M is the number of centers/atoms. + # D is the number of dimensions, usually 3. + # K is maximum number of gaussian functions over all M atoms. + cm, em, coords = self.c_m, self.e_m, self.coords + # Shape (N, M, D), then Summing gives (N, M, 1) + distance = np.sum((points - coords[:, np.newaxis]) ** 2.0, axis=2, keepdims=True) + # At each center, multiply Each Distance of a Coordinate, with its exponents. + exponen = np.exp(-np.einsum("MND, MK-> MNK", distance, em)) + # At each center, multiply the exponential with its coefficients. + gaussian = np.einsum("MNK, MK -> MNK", exponen, cm) + # At each point, sum for each center, then sum all centers together. + return np.einsum("MNK -> N", gaussian, dtype=np.float64) + + def helper_for_derivatives(self, diff_squared, diff_coords, i_var): + r""" + Return Arrays for computing the derivative of transformation functions wrt x, y, z. + + Parameters + ---------- + diff_squared : np.ndarray + The squared of difference of position to the center of the Promoleculars. + diff_coords : np.ndarray + The difference of position to the center of the Promoleculars. This is the square + root of `diff_squared`. + i_var : int + Index of one of x, y, z. + + Returns + ------- + distance : np.ndarray + The squared distance from the position to the center of the Promoleculars. + single_gauss : np.ndarray + Array with entries of a single Gaussian e^(-a distance) with factor (pi / a). + integrate_till_pt_x : np.ndarray + Integration of a Gaussian from -inf to x, ie + (pi / a)^0.5 * (erf(a^0.5 (x - center of Gaussian) + 1) / 2 + transf_num : float + The numerator of the transformation. Mostly used for quotient rule. + transf_den : float + The denominator of the transformation. Mostly used for quotient rule. + + """ + distance = np.sum(diff_squared[:, :i_var], axis=1)[:, np.newaxis] + + # Gaussian Integrals Over Entire Space For Numerator and Denomator. + single_gauss = self.single_gaussians(distance) + single_gauss *= self.pi_over_exponents ** (self.dim - i_var - 1) + + # Get integral of Gaussian till a point. + integrate_till_pt_x = self.integration_gaussian_till_point( + diff_coords, i_var, with_factor=True + ) + # Numerator and Denominator of Original Transformation. + transf_num = np.sum(single_gauss * integrate_till_pt_x) + transf_den = np.sum(single_gauss * self.pi_over_exponents) + return distance, single_gauss, integrate_till_pt_x, transf_num, transf_den + + +def _transform_coordinate(real_pt, i_var, promol): + r""" + Transform the `i_var` coordinate of a real point to [-1, 1] using promolecular density. + + Parameters + ---------- + real_pt : np.ndarray(D,) + Real point being transformed. + i_var : int + Index that is being tranformed. Less than D. + promol : _PromolParams + Promolecular Data Class. + + Returns + ------- + theta_pt : float + The transformed point in :math:`[-1, 1]`. + + """ + _, _, coords, pi_over_exps, dim = astuple(promol) + + # Distance to centers/nuclei`s and Prefactors. + diff_coords = real_pt[: i_var + 1] - coords[:, : i_var + 1] + diff_squared = diff_coords**2.0 + distance = np.sum(diff_squared[:, :i_var], axis=1)[:, np.newaxis] + # If i_var is zero, then distance is just all zeros. + + # Single Gaussians Including Integration of Exponential over `(dim - i_var)` variables. + single_gauss = promol.single_gaussians(distance) * pi_over_exps ** (dim - i_var) + + # Get the integral of Gaussian till a point excluding a prefactor. + # prefactor (pi / exponents) is included in `gaussian_integrals`. + cdf_gauss = promol.integration_gaussian_till_point(diff_coords, i_var, with_factor=False) + + # Final Result. + transf_num = np.sum(single_gauss * cdf_gauss) + transf_den = np.sum(single_gauss) + transform_value = transf_num / transf_den + + # -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): + r""" + Equation to solve for the root to find inverse coordinate from theta space to Real space. + + Parameters + ---------- + init_guess : float + Initial guess of Real point that transforms to `theta_pt`. + prev_trans_pts : list[`i_var` - 1] + The previous points in real-space that were already transformed. + theta_pt : float + The point in [-1, 1] being transformed to the Real space. + i_var : int + Index of variable being transformed. + promol : _PromolParams + Promolecular Data Class. + + Returns + ------- + float : + The difference between `theta_pt` and the transformed point based on + `init_guess` and `prev_trans_pts`. + + """ + all_points = np.append(prev_trans_pts, init_guess) + transf_pt = _transform_coordinate(all_points, i_var, promol) + return theta_pt - transf_pt + + +def _inverse_coordinate(theta_pt, i_var, transformed, promol, bracket=(-10, 10)): + r""" + Transform a point in [-1, 1] to the real space corresponding to `i_var` variable. + + Parameters + ---------- + theta_pt : float + Point in [-1, 1]. + i_var : int + Index that is being tranformed. Less than D. + transformed : list[`i_var` - 1] + The set of previous points before index `i_var` that were transformed to real space. + promol : _PromolParams + Promolecular Data Class. + bracket : (float, float) + Interval where root is suspected to be in Reals. Used for "brentq" root-finding method. + Default is (-10, 10). + + Returns + ------- + real_pt : float + Return the transformed real point. + + Raises + ------ + AssertionError : If the root did not converge, or brackets did not have opposite sign. + RuntimeError : If dynamic bracketing reached maximum iteration. + + Notes + ----- + - If the theta point is on the boundary or it is itself a nan, then it get's mapped to nan. + Further, if nan is in `transformed[:i_var]` then this function will return nan. + - If Brackets do not have the opposite sign, will change the brackets by adding/subtracting + the value 10 to the lower or upper bound that is closest to zero. + + """ + # 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 - -1.0) < 1e-10: + return np.inf + if np.abs(theta_pt - 1.0) < 1e-10: + return np.inf + # This condition is added for transformation of the entire grid. + # The [:i_var] is needed because of the way I've set-up transforming points in _transform. + # Likewise for the bracket, see the function `get_bracket`. + if np.inf in bracket or np.inf in transformed[:i_var]: + return np.inf + + def _dynamic_bracketing(l_bnd, u_bnd, maxiter=5000): + r"""Dynamically changes the lower (or upper bound) to have different sign values.""" + bounds = [l_bnd, u_bnd] + + def is_same_sign(x, y): + return (x >= 0 and y >= 0) or (x < 0 and y < 0) + + f_l_bnd = _root_equation(l_bnd, *args) + f_u_bnd = _root_equation(u_bnd, *args) + + # Get Index of the one that is closest to zero, the one that needs to change. + f_bnds = np.abs([f_l_bnd, f_u_bnd]) + idx = f_bnds.argmin() + # Check if they have the same sign. + same_sign = is_same_sign(*bounds) + counter = 0 + while same_sign and counter < maxiter: + # Add 10 to the upper bound or subtract 10 to the lower bound to the one that + # is closest to zero. This is done based on the sign. + bounds[idx] = np.sign(idx - 0.5) * 10 + bracket[idx] + # Update info for next iteration. + if idx == 0: + f_l_bnd = _root_equation(bracket[0], *args) + else: + f_u_bnd = _root_equation(bracket[1], *args) + same_sign = is_same_sign(f_l_bnd, f_u_bnd) + counter += 1 + + if counter == maxiter: + raise RuntimeError( + f"Couldn't find correct bounds {bracket} for the root-solver " + f"to solve for the inverse." + ) + return tuple(bounds) + + # Set up Arguments for root_equation with dynamic bracketing. + args = (transformed[:i_var], theta_pt, i_var, promol) + bracket = _dynamic_bracketing(bracket[0], bracket[1]) + root_result = root_scalar( + _root_equation, + args=args, + method="brentq", + bracket=[bracket[0], bracket[1]], + maxiter=50, + xtol=2e-15, + ) + assert root_result.converged + return root_result.root + + +def _pad_coeffs_exps_with_zeros(coeffs, exps): + r"""Pad Promolecular coefficients and exponents with zero. Results in same size array.""" + max_numb_of_gauss = max(len(c) for c in coeffs) + coeffs = np.array( + [ + np.pad(a, (0, max_numb_of_gauss - len(a)), "constant", constant_values=0.0) + for a in coeffs + ], + dtype=np.float64, + ) + exps = np.array( + [np.pad(a, (0, max_numb_of_gauss - len(a)), "constant", constant_values=0.0) for a in exps], + dtype=np.float64, + ) + return coeffs, exps diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py index 42f3f30b..64820c20 100644 --- a/src/grid/tests/test_cubic.py +++ b/src/grid/tests/test_cubic.py @@ -188,7 +188,7 @@ def test_point_and_weights_are_correct(self): assert_allclose(actual_weight, cubic.weights[index]) index += 1 - def test_interpolation_of_gaussian_vertorized(self): + def test_interpolation_of_gaussian_vectorized(self): r"""Test interpolation of a Gaussian function with vectorization.""" oned = MidPoint(50) cubic = Tensor1DGrids(oned, oned, oned) diff --git a/src/grid/tests/test_protransform.py b/src/grid/tests/test_protransform.py new file mode 100644 index 00000000..9a003798 --- /dev/null +++ b/src/grid/tests/test_protransform.py @@ -0,0 +1,744 @@ +# GRID is a numerical integration module for quantum chemistry. +# +# Copyright (C) 2011-2019 The GRID Development Team +# +# This file is part of GRID. +# +# GRID is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# GRID is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# -- +from grid.basegrid import OneDGrid +from grid.onedgrid import GaussChebyshevLobatto +from grid.protransform import ( + CubicProTransform, + _PromolParams, + _pad_coeffs_exps_with_zeros, + _transform_coordinate, +) + +import numpy as np +from numpy.testing import assert_allclose + +import pytest + +from scipy.optimize import approx_fprime +from scipy.special import erf + + +class TestTwoGaussianDiffCenters: + r"""Test a Sum of Two Gaussian function against analytic formulas and numerical procedures.""" + + def setUp(self, ss=0.1, return_obj=False, add_boundary=True): + r"""Set up a two parameter Gaussian function.""" + c = np.array([[5.0], [10.0]]) + e = np.array([[2.0], [3.0]]) + coord = np.array([[1.0, 2.0, 3.0], [2.0, 2.0, 2.0]]) + params = _PromolParams(c, e, coord, dim=3) + if return_obj: + num_pts = int(2 / ss) + 1 + weights = np.array([2.0 / (num_pts - 2)] * num_pts) + if add_boundary: + l_bnd = -1.0 + end_point = True + else: + l_bnd = -0.99 + end_point = False + oned = OneDGrid( + np.linspace(l_bnd, 1, num=num_pts, endpoint=end_point), + weights, + domain=(-1, 1), + ) + + obj = CubicProTransform( + [oned, oned, oned], + params.c_m, + params.e_m, + params.coords, + ) + return params, obj + return params + + def promolecular(self, x, y, z, params): + r"""Hard-Code the promolecular density.""" + # Promolecular in UniformProTransform class uses einsum, this tests it against that. + # Also could be used for integration tests. + cm, em, coords = params.c_m, params.e_m, params.coords + promol = 0.0 + for i, coeffs in enumerate(cm): + xc, yc, zc = coords[i] + for j, coeff in enumerate(coeffs): + distance = (x - xc) ** 2.0 + (y - yc) ** 2.0 + (z - zc) ** 2.0 + promol += np.sum(coeff * np.exp(-em[i, j] * distance)) + return promol + + def test_promolecular_density(self): + r"""Test Promolecular Density function against hard-coded.""" + grid = np.array( + [ + [-1.0, 2.0, 3.0], + [5.0, 10.0, -1.0], + [0.0, 0.0, 0.0], + [3.0, 2.0, 1.0], + [10.0, 20.0, 30.0], + [0.0, 10.0, 0.2], + ] + ) + params = self.setUp() + + true_ans = [] + for pt in grid: + x, y, z = pt + true_ans.append(self.promolecular(x, y, z, params)) + + desired = params.promolecular(grid) + assert np.all(np.abs(np.array(true_ans) - desired) < 1e-8) + + @pytest.mark.parametrize( + "pts, add_boundary", [[np.arange(-5.0, 5.0, 0.5), False], [np.arange(-5.0, 5.0, 0.5), True]] + ) + def test_transforming_x_against_formula(self, pts, add_boundary): + r"""Test transformming the X-transformation against analytic formula.""" + for pt in pts: + true_ans = _transform_coordinate([pt], 0, self.setUp(add_boundary=add_boundary)) + + def formula_transforming_x(x): + r"""Return closed form formula for transforming x coordinate.""" + first_factor = (5.0 * np.pi**1.5 / (4 * 2**0.5)) * (erf(2**0.5 * (x - 1)) + 1.0) + + sec_fac = ((10.0 * np.pi**1.5) / (6.0 * 3**0.5)) * (erf(3.0**0.5 * (x - 2)) + 1.0) + + 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 + + assert np.abs(true_ans - formula_transforming_x(pt)) < 1e-8 + + @pytest.mark.parametrize( + "pts_xy, add_boundary", + [ + [np.random.uniform(-10.0, 10.0, size=(100, 2)), False], + [np.random.uniform(-10.0, 10.0, size=(100, 2)), True], + ], + ) + def test_transforming_y_against_formula(self, pts_xy, add_boundary): + r"""Test transforming the Y-transformation against analytic formula.""" + for x, y in pts_xy: + true_ans = _transform_coordinate([x, y], 1, self.setUp(add_boundary=add_boundary)) + + def formula_transforming_y(x, y): + r"""Return closed form formula for transforming y coordinate.""" + fac1 = 5.0 * np.sqrt(np.pi / 2.0) * np.exp(-2.0 * (x - 1) ** 2) + fac1 *= np.sqrt(np.pi) * (erf(2.0**0.5 * (y - 2)) + 1.0) / (2.0 * np.sqrt(2.0)) + fac2 = 10.0 * np.sqrt(np.pi / 3.0) * np.exp(-3.0 * (x - 2.0) ** 2.0) + fac2 *= np.sqrt(np.pi) * (erf(3.0**0.5 * (y - 2)) + 1.0) / (2.0 * np.sqrt(3.0)) + num = fac1 + fac2 + + 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 -1.0 + 2.0 * (num / den) + + assert np.abs(true_ans - formula_transforming_y(x, y)) < 1e-8 + + @pytest.mark.parametrize( + "pts, add_boundary", + [ + [np.random.uniform(-10.0, 10.0, size=(100, 3)), False], + [np.random.uniform(-10.0, 10.0, size=(100, 3)), True], + ], + ) + def test_transforming_z_against_formula(self, pts, add_boundary): + r"""Test transforming the Z-transformation against analytic formula.""" + params, obj = self.setUp(ss=0.5, return_obj=True, add_boundary=add_boundary) + + def formula_transforming_z(x, y, z): + r"""Return closed form formula for transforming z coordinate.""" + a1, a2, a3 = (x - 1.0), (y - 2.0), (z - 3.0) + erfx = erf(2.0**0.5 * a3) + 1.0 + fac1 = 5.0 * np.exp(-2.0 * (a1**2.0 + a2**2.0)) * erfx * np.pi**0.5 / (2.0 * 2.0**0.5) + b1, b2, b3 = (x - 2.0), (y - 2.0), (z - 2.0) + erfy = erf(3.0**0.5 * b3) + 1.0 + fac2 = 10.0 * np.exp(-3.0 * (b1**2.0 + b2**2.0)) * erfy * np.pi**0.5 / (2.0 * 3.0**0.5) + 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 -1.0 + 2.0 * (fac1 + fac2) / den + + for x, y, z in pts: + true_ans = formula_transforming_z(x, y, z) + # Test function + actual = _transform_coordinate([x, y, z], 2, params) + assert np.abs(true_ans - actual) < 1e-8 + + # Test Method + actual = obj.transform(np.array([x, y, z]))[2] + assert np.abs(true_ans - actual) < 1e-8 + + def test_transforming_simple_grid(self): + r"""Test transforming a grid that only contains one non-boundary point.""" + ss = 1.0 + params, obj = self.setUp(ss, return_obj=True) + num_pt = int(2.0 / ss) + 1 # number of points in one-direction. + assert obj.points.shape == (num_pt**3, 3) + non_boundary_pt_index = num_pt**2 + num_pt + 1 + real_pt = obj.points[non_boundary_pt_index] + # Test that this point is not the boundary. + assert real_pt[0] != np.inf + assert real_pt[1] != np.inf + assert real_pt[2] != np.inf + # 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) < 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.isinf(all_nans), axis=1)) + + 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. + # So this isn't really the inverse, in the mathematical sense. + param, obj = self.setUp(0.5, return_obj=True) + + pt = np.array([1, 2, 3], dtype=np.float64) + transf = obj.transform(pt) + reverse = obj.inverse(transf) + assert np.all(np.abs(reverse - pt) < 1e-10) + + def test_integrating_itself(self): + 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: + promol.append(self.promolecular(pt[0], pt[1], pt[2], params)) + promol = np.array(promol, dtype=np.float64) + desired = obj.prointegral + actual = obj.integrate(promol) + assert np.abs(actual - desired) < 1e-8 + + actual = obj.integrate(promol, trick=True) + assert np.abs(actual - desired) < 1e-8 + + @pytest.mark.parametrize("pts", [np.arange(-5.0, 5.0, 0.75)]) + def test_derivative_tranformation_x_finite_difference(self, pts): + r"""Test the derivative of X-transformation against finite-difference.""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for pt in pts: + pt = np.array([pt, 2.0, 3.0]) + + actual = obj.jacobian(pt) + + def tranformation_x(pt): + return _transform_coordinate(pt, 0, params) + + grad = approx_fprime([pt[0]], tranformation_x, 1e-6) + assert np.abs(grad - actual[0, 0]) < 1e-4 + + @pytest.mark.parametrize("pts_xy", [np.random.uniform(-5.0, 5.0, size=(100, 2))]) + def test_derivative_tranformation_y_finite_difference(self, pts_xy): + r"""Test the derivative of Y-transformation against finite-difference.""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for x, y in pts_xy: + actual = obj.jacobian(np.array([x, y, 3.0])) + + def tranformation_y(pt): + return _transform_coordinate([x, pt[0]], 1, params) + + grad = approx_fprime([y], tranformation_y, 1e-8) + assert np.abs(grad - actual[1, 1]) < 1e-5 + + def transformation_y_wrt_x(pt): + return _transform_coordinate([pt[0], y], 1, params) + + h = 1e-8 + deriv = np.imag(transformation_y_wrt_x([complex(x, h)])) / h + assert np.abs(deriv - actual[1, 0]) < 1e-4 + + @pytest.mark.parametrize("pts", [np.random.uniform(-3.0, 3.0, size=(100, 3))]) + def test_derivative_tranformation_z_finite_difference(self, pts): + r"""Test the derivative of Z-transformation against finite-difference.""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for x, y, z in pts: + actual = obj.jacobian(np.array([x, y, z])) + + def tranformation_z(pt): + return _transform_coordinate([x, y, pt[0]], 2, params) + + grad = approx_fprime([z], tranformation_z, 1e-8) + assert np.abs(grad - actual[2, 2]) < 1e-5 + + def transformation_z_wrt_y(pt): + return _transform_coordinate([x, pt[0], z], 2, params) + + deriv = approx_fprime([y], transformation_z_wrt_y, 1e-8) + assert np.abs(deriv - actual[2, 1]) < 1e-4 + + def transformation_z_wrt_x(pt): + a = _transform_coordinate([pt[0], y, z], 2, params) + return a + + h = 1e-8 + deriv = np.imag(transformation_z_wrt_x([complex(x, h)])) / h + + assert np.abs(deriv - actual[2, 0]) < 1e-4 + + @pytest.mark.parametrize("pts", [np.random.uniform(-3.0, 3.0, size=(100, 3))]) + def test_steepest_ascent_direction_with_numerics(self, pts): + r"""Test steepest-ascent direction match in real and theta space. + + The function to test is x^2 + y^2 + z^2. + """ + params, obj = self.setUp(ss=0.2, return_obj=True) + + def grad(pt): + # Gradient of x^2 + y^2 + z^2. + return np.array([2.0 * pt[0], 2.0 * pt[1], 2.0 * pt[2]]) + + for x, y, z in pts: + # Take a step in real-space. + pt = np.array([x, y, z]) + grad_pt = grad(pt) + step = 1e-8 + pt_step = pt + grad_pt * step + + # Convert the steps in theta space and calculate finite-difference gradient. + transf = obj.transform(pt) + transf_step = obj.transform(pt_step) + grad_finite = (transf_step - transf) / step + + # Test the actual actual. + actual = obj.steepest_ascent_theta(pt, grad_pt) + assert np.all(np.abs(actual - grad_finite) < 1e-4) + + def test_second_derivative_of_theta_x_dx(self): + r"""Test the second derivative of d(theta_x)/d(x) .""" + params, obj = self.setUp(ss=0.2, return_obj=True) + x, y, z = 0.0, 0.0, 0.0 + actual = obj.hessian(np.array([x, y, z])) + + # test second derivative of theta_x wrt to dx dx. + def tranformation_x(pt): + return obj.jacobian([pt[0], y, z])[0, 0] + + dthetax_dx_dx = approx_fprime([x], tranformation_x, 1e-8) + assert np.abs(dthetax_dx_dx - actual[0, 0, 0]) < 1e-5 + # Since theta_x only depends on x, then all other elements + # are zero. + for i in range(0, 3): + for j in range(0, 3): + if i != j: + assert np.abs(actual[0, i, j]) < 1e-10 + + @pytest.mark.parametrize("pts", [np.random.uniform(-3.0, 3.0, size=(100, 3))]) + def test_second_derivative_of_theta_y(self, pts): + r"""Test the second derivative of d(theta_y).""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for x, y, z in pts: + actual = obj.hessian(np.array([x, y, z])) + + # test second derivative of theta_y wrt to dy dy. + def dtheta_y_dy(pt): + return obj.jacobian([x, pt[0], z])[1, 1] + + dthetay_dy_dy = approx_fprime([y], dtheta_y_dy, 1e-8) + assert np.abs(dthetay_dy_dy - actual[1, 1, 1]) < 1e-5 + + # test second derivative of d(theta_y)/d(y) wrt dx. + def dtheta_y_dy(pt): + return obj.jacobian([pt[0], y, z])[1, 1] + + dthetay_dy_dx = approx_fprime([x], dtheta_y_dy, 1e-8) + assert np.abs(dthetay_dy_dx - actual[1, 1, 0]) < 1e-5 + + # test second derivative of d(theta_y)/(dx) wrt dy + def dtheta_y_dx(pt): + return obj.jacobian([x, pt[0], z])[1, 0] + + dtheta_y_dx_dy = approx_fprime([y], dtheta_y_dx, 1e-8) + assert np.abs(dtheta_y_dx_dy - actual[1, 0, 1]) < 1e-4 + + # test second derivative of d(theta_y)/(dx) wrt dx + def dtheta_y_dx(pt): + return obj.jacobian([pt[0], y, z])[1, 0] + + dtheta_y_dx_dx = approx_fprime([x], dtheta_y_dx, 1e-8) + assert np.abs(dtheta_y_dx_dx - actual[1, 0, 0]) < 1e-5 + + # Test other elements are all zeros. + for i in range(0, 3): + assert np.abs(actual[1, 2, i]) < 1e-10 + assert np.abs(actual[1, 0, 2]) < 1e-10 + assert np.abs(actual[1, 1, 2]) < 1e-10 + + @pytest.mark.parametrize("pts", [np.random.uniform(-3.0, 3.0, size=(100, 3))]) + def test_second_derivative_of_theta_z(self, pts): + r"""Test the second derivative of d(theta_z).""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for x, y, z in pts: + actual = obj.hessian(np.array([x, y, z])) + + # test second derivative of theta_z wrt to dz dz. + def dtheta_z_dz(pt): + return obj.jacobian([x, y, pt[0]])[2, 2] + + dthetaz_dz_dz = approx_fprime([z], dtheta_z_dz, 1e-8) + assert np.abs(dthetaz_dz_dz - actual[2, 2, 2]) < 1e-5 + + # test second derivative of theta_z wrt to dz dy. + def dtheta_z_dz(pt): + return obj.jacobian([x, pt[0], z])[2, 2] + + dthetaz_dz_dy = approx_fprime([y], dtheta_z_dz, 1e-8) + assert np.abs(dthetaz_dz_dy - actual[2, 2, 1]) < 1e-5 + + # test second derivative of theta_z wrt to dz dx. + def dtheta_z_dz(pt): + return obj.jacobian([pt[0], y, z])[2, 2] + + dthetaz_dz_dx = approx_fprime([x], dtheta_z_dz, 1e-8) + assert np.abs(dthetaz_dz_dx - actual[2, 2, 0]) < 1e-5 + + # test second derivative of theta_z wrt to dy dz. + def dtheta_z_dy(pt): + return obj.jacobian([x, y, pt[0]])[2, 1] + + dthetaz_dy_dz = approx_fprime([z], dtheta_z_dy, 1e-8) + assert np.abs(dthetaz_dy_dz - actual[2, 1, 2]) < 1e-5 + + # test second derivative of theta_z wrt to dy dy. + def dtheta_z_dy(pt): + return obj.jacobian([x, pt[0], z])[2, 1] + + dthetaz_dy_dy = approx_fprime([y], dtheta_z_dy, 1e-8) + assert np.abs(dthetaz_dy_dy - actual[2, 1, 1]) < 1e-5 + + # test second derivative of theta_z wrt to dy dx. + def dtheta_z_dy(pt): + return obj.jacobian([pt[0], y, z])[2, 1] + + dthetaz_dy_dx = approx_fprime([x], dtheta_z_dy, 1e-8) + assert np.abs(dthetaz_dy_dx - actual[2, 1, 0]) < 1e-5 + + # test second derivative of theta_z wrt to dx dz. + def dtheta_z_dx(pt): + return obj.jacobian([x, y, pt[0]])[2, 0] + + dthetaz_dx_dz = approx_fprime([z], dtheta_z_dx, 1e-8) + assert np.abs(dthetaz_dx_dz - actual[2, 0, 2]) < 1e-5 + + # test second derivative of theta_z wrt to dx dy. + def dtheta_z_dx(pt): + return obj.jacobian([x, pt[0], z])[2, 0] + + dthetaz_dx_dy = approx_fprime([y], dtheta_z_dx, 1e-8) + assert np.abs(dthetaz_dx_dy - actual[2, 0, 1]) < 1e-5 + + # test second derivative of theta_z wrt to dx dx. + def dtheta_z_dx(pt): + return obj.jacobian([pt[0], y, z])[2, 0] + + dthetaz_dx_dx = approx_fprime([x], dtheta_z_dx, 1e-8) + assert np.abs(dthetaz_dx_dx - actual[2, 0, 0]) < 1e-5 + + +class TestOneGaussianAgainstNumerics: + r"""Tests With Numerical Integration of a One Gaussian function.""" + + def setUp(self, ss=0.1, return_obj=False): + r"""Return a one Gaussian example and its UniformProTransform object.""" + c = np.array([[5.0]]) + e = np.array([[2.0]]) + coord = np.array([[1.0, 2.0, 3.0]]) + params = _PromolParams(c, e, coord, dim=3) + if return_obj: + num_pts = int(1 / ss) + 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_x, oned_x], params.c_m, params.e_m, params.coords) + return params, obj + return params + + @pytest.mark.parametrize("pts", [np.arange(-5.0, 5.0, 0.5)]) + def test_transforming_x_against_numerics(self, pts): + r"""Test transforming X against numerical algorithms.""" + for pt in pts: + true_ans = _transform_coordinate([pt], 0, self.setUp()) + + def promolecular_in_x(grid, every_grid): + r"""Construct the formula of promolecular for integration.""" + promol_x = 5.0 * np.exp(-2.0 * (grid - 1.0) ** 2.0) + promol_x_all = 5.0 * np.exp(-2.0 * (every_grid - 1.0) ** 2.0) + return promol_x, promol_x_all + + 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 = -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("pts_xy", [np.random.uniform(-10.0, 10.0, size=(100, 2))]) + def test_transforming_y_against_numerics(self, pts_xy): + r"""Test transformation y against numerical algorithms.""" + + def promolecular_in_y(grid, every_grid): + r"""Construct the formula of promolecular for integration.""" + promol_y = 5.0 * np.exp(-2.0 * (grid - 2.0) ** 2.0) + promol_y_all = 5.0 * np.exp(-2.0 * (every_grid - 2.0) ** 2.0) + return promol_y_all, promol_y + + for x, y in pts_xy: + true_ans = _transform_coordinate([x, y], 1, self.setUp()) + + 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 = -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("pts", [np.random.uniform(-10.0, 10.0, size=(100, 3))]) + def test_transforming_z_against_numerics(self, pts): + r"""Test transforming Z against numerical algorithms.""" + + def promolecular_in_z(grid, every_grid): + r"""Construct the formula of promolecular for integration.""" + promol_z = 5.0 * np.exp(-2.0 * (grid - 3.0) ** 2.0) + promol_z_all = 5.0 * np.exp(-2.0 * (every_grid - 3.0) ** 2.0) + return promol_z_all, promol_z + + for x, y, z in pts: + 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 = -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-4 + + @pytest.mark.parametrize("pts", [np.random.uniform(-3.0, 3.0, size=(100, 3))]) + def test_jacobian(self, pts): + r"""Test that the jacobian of the transformation.""" + params, obj = self.setUp(ss=0.2, return_obj=True) + for x, y, z in pts: + actual = obj.jacobian(np.array([x, y, z])) + + # assert lower-triangular component is zero. + assert np.abs(actual[1, 0]) < 1e-5 + assert np.abs(actual[2, 0]) < 1e-5 + assert np.abs(actual[2, 1]) < 1e-5 + + # test derivative wrt to x + def tranformation_x(pt): + return _transform_coordinate([pt[0], y, z], 0, params) + + grad = approx_fprime([x], tranformation_x, 1e-8) + assert np.abs(grad - actual[0, 0]) < 1e-5 + + # test derivative wrt to y + def tranformation_y(pt): + return _transform_coordinate([x, pt[0]], 1, params) + + grad = approx_fprime([y], tranformation_y, 1e-8) + assert np.abs(grad - actual[1, 1]) < 1e-5 + + # Test derivative wrt to z + def tranformation_z(pt): + return _transform_coordinate([x, y, pt[0]], 2, params) + + grad = approx_fprime([z], tranformation_z, 1e-8) + assert np.abs(grad - actual[2, 2]) < 1e-5 + + +class TestInterpolation: + r"""Test Interpolation Methods and Cubic Grid Utility Functions.""" + + def setUp(self, ss=0.1): + r"""Set up a simple one Gaussian function with uniform cubic grid.""" + c = np.array([[5.0]]) + e = np.array([[2.0]]) + coord = np.array([[1.0, 2.0, 3.0]]) + params = _PromolParams(c, e, coord, dim=3) + + num_pts = int(2 / ss) + 1 + weights = np.array([2.0 / (num_pts - 2)] * num_pts) + oned_x = OneDGrid( + np.linspace(-1.0, 1.0, num=num_pts, endpoint=True), weights, domain=(-1, 1) + ) + + obj = CubicProTransform([oned_x, oned_x, oned_x], params.c_m, params.e_m, params.coords) + return params, obj, [oned_x, oned_x, oned_x] + + def test_interpolate_cubic_function(self): + r"""Interpolate a cubic function.""" + param, obj, oned_grids = self.setUp(ss=0.08) + + # Function to interpolate. + def func(pts): + return (pts[:, 0] - 1.0) ** 3.0 + (pts[:, 1] - 2.0) ** 3.0 + (pts[:, 2] - 3.0) ** 3.0 + + # Test over a grid that is very close to the nucleus. + grid = np.vstack( + ( + np.random.uniform(1.0, 1.5, (100,)).T, + np.random.uniform(1.5, 2.5, (100,)).T, + np.random.uniform(2.5, 3.5, (100,)).T, + ) + ).T + actuals = obj.interpolate(grid, func(obj.points), oned_grids, use_log=False) + desired = func(grid) + assert_allclose(desired, actuals, atol=1e-3, rtol=1e-4) + + def test_interpolate_gaussian(self): + r"""Interpolate a Gaussian function with use_log set to True.""" + param, obj, oned_grids = self.setUp(ss=0.08) + + # Function to interpolate. + def func(pts, alpha=2.0): + return np.exp(-alpha * np.linalg.norm(pts - param.coords[0], axis=1) ** 2.0) + + # Test over a grid. Pytest isn't used for effiency reasons. + # TODO: the grid points need to be close to center to achieve good accuracy. + real_grid = np.vstack( + ( + np.random.uniform(0.75, 1.25, (100,)).T, + np.random.uniform(0.75, 2.25, (100,)).T, + np.random.uniform(3.75, 3.25, (100,)).T, + ) + ).T + actuals = obj.interpolate(real_grid, func(obj.points), oned_grids, use_log=True) + desired = func(real_grid) + assert_allclose(desired, actuals, atol=1e-1) + + def test_interpolate_derivative_cubic_function(self): + r"""Interpolate the derivative of some simple function.""" + param, obj, oned_grids = self.setUp(ss=0.08) + + # Function to interpolate. + def func(pts): + return (pts[:, 0] - 1.0) * (pts[:, 1] - 2.0) * (pts[:, 2] - 3.0) + + def derivative(pts): + return np.vstack( + [ + (pts[:, 1] - 2.0) * (pts[:, 2] - 3.0), + (pts[:, 0] - 1.0) * (pts[:, 2] - 3.0), + (pts[:, 0] - 1.0) * (pts[:, 1] - 2.0), + ] + ).T + + # Test over a grid. Pytest isn't used for effiency reasons. + grid = np.vstack( + ( + np.random.uniform(1.0, 1.5, (100,)).T, + np.random.uniform(1.5, 2.5, (100,)).T, + np.random.uniform(2.5, 3.5, (100,)).T, + ) + ).T + actual = obj.interpolate(grid, func(obj.points), oned_grids, nu=1) + desired = derivative(grid) + assert_allclose(actual, desired, atol=1e-2) + + def test_interpolate_derivative_cubic_function2(self): + r"""Interpolate the derivative of some simple function.""" + param, obj, oned_grids = self.setUp(ss=0.08) + + # Function to interpolate. + def func(pts): + return (pts[:, 0] - 1.0) ** 2.0 * (pts[:, 1] - 2.0) ** 2.0 * (pts[:, 2] - 3.0) ** 2.0 + + def derivative(pts): + x, y, z = pts[:, 0], pts[:, 1], pts[:, 2] + return np.vstack( + [ + 2.0 * (x - 1.0) * (y - 2.0) ** 2.0 * (z - 3.0) ** 2.0, + 2.0 * (x - 1.0) ** 2.0 * (y - 2.0) * (z - 3.0) ** 2.0, + 2.0 * (x - 1.0) ** 2.0 * (y - 2.0) ** 2.0 * (z - 3.0), + ] + ).T + + # Test over a grid. Pytest isn't used for effiency reasons. + grid = np.vstack( + ( + np.random.uniform(1.0, 1.5, (100,)).T, + np.random.uniform(1.5, 2.5, (100,)).T, + np.random.uniform(2.5, 3.5, (100,)).T, + ) + ).T + actual = obj.interpolate(grid, func(obj.points), oned_grids, nu=1, use_log=False) + desired = derivative(grid) + assert_allclose(actual, desired, atol=1e-2) + # Test with logarithm set to True + actual = obj.interpolate(grid, func(obj.points), oned_grids, nu=1, use_log=True) + assert_allclose(actual, desired, atol=1e-2) + + +class TestIntegration: + r""" + Only one integration test as of this moment. + + Choose to make it, it's own seperate class since many choices of oned grid is possible. + + """ + + def setUp_one_gaussian(self, ss=0.03): + r"""Return a one Gaussian example and its UniformProTransform object.""" + c = np.array([[5.0]]) + e = np.array([[2.0]]) + coord = np.array([[1.0, 2.0, 3.0]]) + params = _PromolParams(c, e, coord, dim=3) + num_pts = int(1 / ss) + 1 + oned_x = GaussChebyshevLobatto(num_pts) + obj = CubicProTransform([oned_x, oned_x, oned_x], params.c_m, params.e_m, params.coords) + return params, obj + + def test_integration_perturbed_gaussian_with_promolecular_trick(self): + r"""Test integration of a slightly perturbed function.""" + # Only Measured against one decimal place and very similar exponent. + _, obj = self.setUp_one_gaussian(ss=0.03) + + # Gaussian exponent is slightly perturbed from 2. + exponent = 2.01 + + def gaussian(grid): + return 5.0 * np.exp( + -exponent * np.sum((grid - np.array([1.0, 2.0, 3.0])) ** 2.0, axis=1) + ) + + func_vals = gaussian(obj.points) + desired = 5.0 * np.sqrt(np.pi / exponent) ** 3.0 + actual = obj.integrate(func_vals, trick=True) + + assert np.abs(actual - desired) < 1e-3 + + +def test_padding_arrays(): + r"""Test different array sizes are correctly padded.""" + coeff = np.array([[1.0, 2.0], [1.0, 2.0, 3.0, 4.0], [5.0]], dtype=object) + exps = np.array([[4.0, 5.0], [5.0, 6.0, 7.0, 8.0], [9.0]], dtype=object) + coeff_pad, exps_pad = _pad_coeffs_exps_with_zeros(coeff, exps) + coeff_desired = np.array( + [[1.0, 2.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0], [5.0, 0.0, 0.0, 0.0]], dtype=object + ) + np.testing.assert_array_equal(coeff_desired, coeff_pad) + exp_desired = np.array( + [[4.0, 5.0, 0.0, 0.0], [5.0, 6.0, 7.0, 8.0], [9.0, 0.0, 0.0, 0.0]], dtype=object + ) + np.testing.assert_array_equal(exp_desired, exps_pad)