diff --git a/docs/src/references/api/python/misc.rst b/docs/src/references/api/python/misc.rst index 6288741a8..8a08fca83 100644 --- a/docs/src/references/api/python/misc.rst +++ b/docs/src/references/api/python/misc.rst @@ -13,5 +13,3 @@ Miscelaneous .. autoclass:: rascaline.Profiler :members: :undoc-members: - -.. autofunction:: rascaline.generate_splines diff --git a/docs/src/references/api/python/utils/index.rst b/docs/src/references/api/python/utils/index.rst index cdc1a4783..abe7e11a0 100644 --- a/docs/src/references/api/python/utils/index.rst +++ b/docs/src/references/api/python/utils/index.rst @@ -8,3 +8,4 @@ Utility functions and classes that extend the core usage of rascaline. :maxdepth: 1 power-spectrum + splines diff --git a/docs/src/references/api/python/utils/splines.rst b/docs/src/references/api/python/utils/splines.rst new file mode 100644 index 000000000..7e93ba9db --- /dev/null +++ b/docs/src/references/api/python/utils/splines.rst @@ -0,0 +1 @@ +.. automodule:: rascaline.utils.splines diff --git a/python/rascaline-torch/rascaline/torch/__init__.py b/python/rascaline-torch/rascaline/torch/__init__.py index c1bb8ec45..177d0c688 100644 --- a/python/rascaline-torch/rascaline/torch/__init__.py +++ b/python/rascaline-torch/rascaline/torch/__init__.py @@ -37,7 +37,6 @@ __all__ = [ "AtomicComposition", - "CalculatorBase", "LodeSphericalExpansion", "NeighborList", "SoapPowerSpectrum", diff --git a/python/rascaline/examples/splined-radial-integral.py b/python/rascaline/examples/splined-radial-integral.py index d3b5b0f02..24432d7d4 100644 --- a/python/rascaline/examples/splined-radial-integral.py +++ b/python/rascaline/examples/splined-radial-integral.py @@ -1,4 +1,6 @@ """ +.. _example-splines: + Splined radial integrals ======================== @@ -20,7 +22,8 @@ from scipy.special import jv from scipy.special import spherical_jn as j_l -from rascaline import SphericalExpansion, generate_splines +from rascaline import SphericalExpansion +from rascaline.utils import RadialIntegralFromFunction # %% @@ -114,16 +117,16 @@ def laplacian_eigenstate_basis_derivative(n, el, r): # %% # # The radial basis functions and their derivatives can be input into a spline -# generator. This will output the positions of the spline points, the +# generator class. This will output the positions of the spline points, the # values of the basis functions evaluated at the spline points, and the # corresponding derivatives. -spline_points = generate_splines( - laplacian_eigenstate_basis, - laplacian_eigenstate_basis_derivative, - max_radial, - max_angular, - cutoff, - requested_accuracy=1e-5, +spliner = RadialIntegralFromFunction( + radial_integral=laplacian_eigenstate_basis, + radial_integral_derivative=laplacian_eigenstate_basis_derivative, + spline_cutoff=cutoff, + max_radial=max_radial, + max_angular=max_angular, + accuracy=1e-5, ) # %% @@ -136,11 +139,7 @@ def laplacian_eigenstate_basis_derivative(n, el, r): "max_radial": max_radial, "max_angular": max_angular, "center_atom_weight": 0.0, - "radial_basis": { - "TabulatedRadialIntegral": { - "points": spline_points, - } - }, + "radial_basis": spliner.compute(), "atomic_gaussian_width": 1.0, # ignored "cutoff_function": {"Step": {}}, } diff --git a/python/rascaline/rascaline/__init__.py b/python/rascaline/rascaline/__init__.py index ea26f9d82..c11deb7a8 100644 --- a/python/rascaline/rascaline/__init__.py +++ b/python/rascaline/rascaline/__init__.py @@ -15,7 +15,6 @@ ) from .log import set_logging_callback # noqa from .profiling import Profiler # noqa -from .splines import generate_splines # noqa from .status import RascalError # noqa from .systems import IntoSystem, SystemBase # noqa from .version import __version__ # noqa @@ -23,7 +22,6 @@ __all__ = [ "AtomicComposition", - "CalculatorBase", "LodeSphericalExpansion", "NeighborList", "SoapPowerSpectrum", diff --git a/python/rascaline/rascaline/splines.py b/python/rascaline/rascaline/splines.py deleted file mode 100644 index 22cd54d68..000000000 --- a/python/rascaline/rascaline/splines.py +++ /dev/null @@ -1,219 +0,0 @@ -from typing import Callable, Optional - -import numpy as np - - -def generate_splines( - radial_basis: Callable[[int, int, np.ndarray], np.ndarray], - radial_basis_derivatives: Callable[[int, int, np.ndarray], np.ndarray], - max_radial: int, - max_angular: int, - cutoff_radius: float, - n_spline_points: Optional[int] = None, - requested_accuracy: float = 1e-8, -): - """Spline generator for tabulated radial integrals. - - Besides some self-explanatory parameters, this function takes as inputs two - functions, namely radial_basis and radial_basis_derivatives. These must be able to - calculate the radial basis functions by taking n, l, and r as their inputs, where n - and l are integers and r is a numpy 1-D array that contains the spline points at - which the radial basis function (or its derivative) needs to be evaluated. These - functions should return a numpy 1-D array containing the values of the radial basis - function (or its derivative) corresponding to the specified n and l, and evaluated - at all points in the r array. If specified, n_spline_points determines how many - spline points will be used for each splined radial basis function. Alternatively, - the user can specify a requested accuracy. Spline points will be added until either - the relative error or the absolute error fall below the requested accuracy on - average across all radial basis functions. - """ - - def value_evaluator_3D(positions): - values = [] - for el in range(max_angular + 1): - for n in range(max_radial): - value = radial_basis(n, el, positions) - values.append(value) - values = np.array(values).T - values = values.reshape(len(positions), max_angular + 1, max_radial) - return values - - def derivative_evaluator_3D(positions): - derivatives = [] - for el in range(max_angular + 1): - for n in range(max_radial): - derivative = radial_basis_derivatives(n, el, positions) - derivatives.append(derivative) - derivatives = np.array(derivatives).T - derivatives = derivatives.reshape(len(positions), max_angular + 1, max_radial) - return derivatives - - if n_spline_points is not None: # if user specifies the number of spline points - positions = np.linspace(0.0, cutoff_radius, n_spline_points) # spline positions - values = value_evaluator_3D(positions) - derivatives = derivative_evaluator_3D(positions) - else: - dynamic_spliner = DynamicSpliner( - 0.0, - cutoff_radius, - value_evaluator_3D, - derivative_evaluator_3D, - requested_accuracy, - ) - positions, values, derivatives = dynamic_spliner.spline() - - # Convert positions, values, derivatives into the appropriate json formats: - spline_points = [] - for position, value, derivative in zip(positions, values, derivatives): - spline_points.append( - { - "position": position, - "values": { - # this is the data representation used by ndarray through serde - "v": 1, - "dim": value.shape, - "data": value.flatten().tolist(), - }, - "derivatives": { - "v": 1, - "dim": derivative.shape, - "data": derivative.flatten().tolist(), - }, - } - ) - - return spline_points - - -class DynamicSpliner: - def __init__( - self, - start: float, - stop: float, - values_fn: Callable[[np.ndarray], np.ndarray], - derivatives_fn: Callable[[np.ndarray], np.ndarray], - requested_accuracy: float = 1e-8, - ) -> None: - """Dynamic spline generator. - - This class can be used to spline any set of functions defined within the - start-stop interval. Cubic Hermite splines - (https://en.wikipedia.org/wiki/Cubic_Hermite_spline) are used. The same spline - points will be used for all functions, and more will be added until either the - relative error or the absolute error fall below the requested accuracy on - average across all functions. The functions are specified via values_fn and - derivatives_fn. These must be able to take a numpy 1D array of positions as - their input, and they must output a numpy array where the first dimension - corresponds to the input positions, while other dimensions are arbitrary and can - correspond to any way in which the target functions can be classified. The - splines can be obtained via the spline method. - """ - - self.start = start - self.stop = stop - self.values_fn = values_fn - self.derivatives_fn = derivatives_fn - self.requested_accuracy = requested_accuracy - - # initialize spline with 11 points - positions = np.linspace(start, stop, 11) - self.spline_positions = positions - self.spline_values = values_fn(positions) - self.spline_derivatives = derivatives_fn(positions) - - self.number_of_custom_axes = len(self.spline_values.shape) - 1 - - def spline(self): - """Calculates and outputs the splines. - - The outputs of this function are, respectively: - - A numpy 1D array containing the spline positions. These are equally spaced in - the start-stop interval. - - A numpy ndarray containing the values of the splined functions at the spline - positions. The first dimension corresponds to the spline positions, while all - subsequent dimensions are consistent with the values_fn and - `get_function_derivative` provided during initialization of the class. - - A numpy ndarray containing the derivatives of the splined functions at the - spline positions, with the same structure as that of the ndarray of values. - """ - - while True: - n_intermediate_positions = len(self.spline_positions) - 1 - - if n_intermediate_positions >= 50000: - raise ValueError( - "Maximum number of spline points reached. \ - There might be a problem with the functions to be splined" - ) - - half_step = (self.spline_positions[1] - self.spline_positions[0]) / 2 - intermediate_positions = np.linspace( - self.start + half_step, self.stop - half_step, n_intermediate_positions - ) - - estimated_values = self._compute_from_spline(intermediate_positions) - new_values = self.values_fn(intermediate_positions) - - mean_absolute_error = np.mean(np.abs(estimated_values - new_values)) - with np.errstate(divide="ignore"): # Ignore divide-by-zero warnings - mean_relative_error = np.mean( - np.abs((estimated_values - new_values) / new_values) - ) - - if ( - mean_absolute_error < self.requested_accuracy - or mean_relative_error < self.requested_accuracy - ): - break - - new_derivatives = self.derivatives_fn(intermediate_positions) - - concatenated_positions = np.concatenate( - [self.spline_positions, intermediate_positions], axis=0 - ) - concatenated_values = np.concatenate( - [self.spline_values, new_values], axis=0 - ) - concatenated_derivatives = np.concatenate( - [self.spline_derivatives, new_derivatives], axis=0 - ) - - sort_indices = np.argsort(concatenated_positions, axis=0) - - self.spline_positions = concatenated_positions[sort_indices] - self.spline_values = concatenated_values[sort_indices] - self.spline_derivatives = concatenated_derivatives[sort_indices] - - return self.spline_positions, self.spline_values, self.spline_derivatives - - def _compute_from_spline(self, positions): - x = positions - delta_x = self.spline_positions[1] - self.spline_positions[0] - n = (np.floor(x / delta_x)).astype(np.int32) - - t = (x - n * delta_x) / delta_x - t_2 = t**2 - t_3 = t**3 - - h00 = 2.0 * t_3 - 3.0 * t_2 + 1.0 - h10 = t_3 - 2.0 * t_2 + t - h01 = -2.0 * t_3 + 3.0 * t_2 - h11 = t_3 - t_2 - - p_k = self.spline_values[n] - p_k_1 = self.spline_values[n + 1] - - m_k = self.spline_derivatives[n] - m_k_1 = self.spline_derivatives[n + 1] - - new_shape = (-1,) + (1,) * self.number_of_custom_axes - h00 = h00.reshape(new_shape) - h10 = h10.reshape(new_shape) - h01 = h01.reshape(new_shape) - h11 = h11.reshape(new_shape) - - interpolated_values = ( - h00 * p_k + h10 * delta_x * m_k + h01 * p_k_1 + h11 * delta_x * m_k_1 - ) - - return interpolated_values diff --git a/python/rascaline/rascaline/utils/__init__.py b/python/rascaline/rascaline/utils/__init__.py index 1d5ab3774..88d7e9e90 100644 --- a/python/rascaline/rascaline/utils/__init__.py +++ b/python/rascaline/rascaline/utils/__init__.py @@ -2,10 +2,8 @@ import metatensor -from .power_spectrum import PowerSpectrum - - -__all__ = ["PowerSpectrum"] +from .power_spectrum import PowerSpectrum # noqa +from .splines import RadialIntegralFromFunction, RadialIntegralSplinerBase # noqa # path that can be used with cmake to access the rascaline library and headers diff --git a/python/rascaline/rascaline/utils/splines.py b/python/rascaline/rascaline/utils/splines.py new file mode 100644 index 000000000..247cbdfc6 --- /dev/null +++ b/python/rascaline/rascaline/utils/splines.py @@ -0,0 +1,403 @@ +""" +Splined radial integrals +======================== + +Classes for generating spline points which can be used as input for the for the radial +integral hyper paramater option. For an complete example for the LE how to use the +classes see :ref:`example-splines`. + +.. autoclass:: rascaline.utils.RadialIntegralSplinerBase + :members: + :show-inheritance: + +.. autoclass:: rascaline.utils.RadialIntegralFromFunction + :members: + :show-inheritance: +""" + +from abc import ABC, abstractmethod +from typing import Callable, Dict, Optional, Union + +import numpy as np + + +class RadialIntegralSplinerBase(ABC): + """Base class for splining arbitrary radial integrals. + + If ``_radial_integral_derivative`` is not implemented in a child class it will + computed based on finite differences. + + :parameter max_angular: number of radial components + :parameter max_radial: number of angular components + :parameter spline_cutoff: cutoff radius for the spline interpolation. This is also + the maximal value that can be interpolated. + :parameter accuracy: accuracy of the numerical integration and the splining. + Accuracy is reached when either the mean absolute error or the mean relative + error gets below the ``accuracy`` threshold. + """ + + def __init__( + self, + max_radial: int, + max_angular: int, + spline_cutoff: float, + accuracy: float, + ): + self.max_radial = max_radial + self.max_angular = max_angular + self.spline_cutoff = spline_cutoff + self.accuracy = accuracy + + @abstractmethod + def _radial_integral(self, n: int, ell: int, positions: np.ndarray) -> np.ndarray: + """Method calculating the radial integral.""" + ... + + @property + def _center_contribution(self) -> Union[None, np.ndarray]: + r"""Contribution of the central atom required for LODE calculations.""" + + return None + + def _radial_integral_derivative( + self, n: int, ell: int, positions: np.ndarray + ) -> np.ndarray: + """Method calculating the derivatice of the radial integral.""" + displacement = 1e-6 + mean_abs_positions = np.abs(positions).mean() + + if mean_abs_positions <= 1.0: + raise ValueError( + "Numerically derivative of the radial integral can not be performed " + "since positions are too small. Mean of the absolute positions is " + f"{mean_abs_positions:.1e} but should be at least 1." + ) + + radial_integral_pos = self._radial_integral( + n, ell, positions + displacement / 2 + ) + radial_integral_neg = self._radial_integral( + n, ell, positions - displacement / 2 + ) + + return (radial_integral_pos - radial_integral_neg) / displacement + + def _value_evaluator_3D( + self, + positions: np.ndarray, + derivative: bool, + ): + values = np.zeros([len(positions), self.max_angular + 1, self.max_radial]) + for ell in range(self.max_angular + 1): + for n in range(self.max_radial): + if derivative: + values[:, ell, n] = self._radial_integral_derivative( + n, ell, positions + ) + else: + values[:, ell, n] = self._radial_integral(n, ell, positions) + + return values + + def compute( + self, + n_spline_points: Optional[int] = None, + ) -> Dict: + """Compute the spline for rascaline's tabulated radial integrals. + + :parameter n_spline_points: Use fixed number of spline points instead of find + the number based on the provided ``accuracy``. + :returns dict: dictionary for the input as the ``radial_basis`` parameter of a + rascaline calculator. + """ + + def value_evaluator_3D(positions): + return self._value_evaluator_3D(positions, derivative=False) + + def derivative_evaluator_3D(positions): + return self._value_evaluator_3D(positions, derivative=True) + + if n_spline_points is not None: + positions = np.linspace(0, self.spline_cutoff, n_spline_points) + values = value_evaluator_3D(positions) + derivatives = derivative_evaluator_3D(positions) + else: + dynamic_spliner = DynamicSpliner( + 0, + self.spline_cutoff, + value_evaluator_3D, + derivative_evaluator_3D, + self.accuracy, + ) + positions, values, derivatives = dynamic_spliner.spline() + + # Convert positions, values, derivatives into the appropriate json formats: + spline_points = [] + for position, value, derivative in zip(positions, values, derivatives): + spline_points.append( + { + "position": position, + "values": { + "v": 1, + "dim": value.shape, + "data": value.flatten().tolist(), + }, + "derivatives": { + "v": 1, + "dim": derivative.shape, + "data": derivative.flatten().tolist(), + }, + } + ) + + parameters = {"points": spline_points} + + center_contribution = self._center_contribution + if center_contribution is not None: + parameters["center_contribution"] = center_contribution + + return {"TabulatedRadialIntegral": parameters} + + +class DynamicSpliner: + def __init__( + self, + start: float, + stop: float, + values_fn: Callable[[np.ndarray], np.ndarray], + derivatives_fn: Callable[[np.ndarray], np.ndarray], + accuracy: float = 1e-8, + ) -> None: + """Dynamic spline generator. + + This class can be used to spline any set of functions defined within the + start-stop interval. Cubic Hermite splines + (https://en.wikipedia.org/wiki/Cubic_Hermite_spline) are used. The same spline + points will be used for all functions, and more will be added until either the + relative error or the absolute error fall below the requested accuracy on + average across all functions. The functions are specified via values_fn and + derivatives_fn. These must be able to take a numpy 1D array of positions as + their input, and they must output a numpy array where the first dimension + corresponds to the input positions, while other dimensions are arbitrary and can + correspond to any way in which the target functions can be classified. The + splines can be obtained via the spline method. + """ + + self.start = start + self.stop = stop + self.values_fn = values_fn + self.derivatives_fn = derivatives_fn + self.requested_accuracy = accuracy + + # initialize spline with 11 points + positions = np.linspace(start, stop, 11) + self.spline_positions = positions + self.spline_values = values_fn(positions) + self.spline_derivatives = derivatives_fn(positions) + + self.number_of_custom_axes = len(self.spline_values.shape) - 1 + + def spline(self): + """Calculates and outputs the splines. + + The outputs of this function are, respectively: - A numpy 1D array containing + the spline positions. These are equally spaced in the start-stop interval. + - A numpy ndarray containing the values of the splined functions at the spline + positions. The first dimension corresponds to the spline positions, while all + subsequent dimensions are consistent with the values_fn and + `get_function_derivative` provided during initialization of the class. + - A numpy ndarray containing the derivatives of the splined functions at the + spline positions, with the same structure as that of the ndarray of values. + """ + + while True: + n_intermediate_positions = len(self.spline_positions) - 1 + + if n_intermediate_positions >= 50000: + raise ValueError( + "Maximum number of spline points reached. \ + There might be a problem with the functions to be splined" + ) + + half_step = (self.spline_positions[1] - self.spline_positions[0]) / 2 + intermediate_positions = np.linspace( + self.start + half_step, self.stop - half_step, n_intermediate_positions + ) + + estimated_values = self._compute_from_spline(intermediate_positions) + new_values = self.values_fn(intermediate_positions) + + mean_absolute_error = np.mean(np.abs(estimated_values - new_values)) + with np.errstate(divide="ignore"): # Ignore divide-by-zero warnings + mean_relative_error = np.mean( + np.abs((estimated_values - new_values) / new_values) + ) + + if ( + mean_absolute_error < self.requested_accuracy + or mean_relative_error < self.requested_accuracy + ): + break + + new_derivatives = self.derivatives_fn(intermediate_positions) + + concatenated_positions = np.concatenate( + [self.spline_positions, intermediate_positions], axis=0 + ) + concatenated_values = np.concatenate( + [self.spline_values, new_values], axis=0 + ) + concatenated_derivatives = np.concatenate( + [self.spline_derivatives, new_derivatives], axis=0 + ) + + sort_indices = np.argsort(concatenated_positions, axis=0) + + self.spline_positions = concatenated_positions[sort_indices] + self.spline_values = concatenated_values[sort_indices] + self.spline_derivatives = concatenated_derivatives[sort_indices] + + return self.spline_positions, self.spline_values, self.spline_derivatives + + def _compute_from_spline(self, positions): + x = positions + delta_x = self.spline_positions[1] - self.spline_positions[0] + n = (np.floor(x / delta_x)).astype(np.int32) + + t = (x - n * delta_x) / delta_x + t_2 = t**2 + t_3 = t**3 + + h00 = 2.0 * t_3 - 3.0 * t_2 + 1.0 + h10 = t_3 - 2.0 * t_2 + t + h01 = -2.0 * t_3 + 3.0 * t_2 + h11 = t_3 - t_2 + + p_k = self.spline_values[n] + p_k_1 = self.spline_values[n + 1] + + m_k = self.spline_derivatives[n] + m_k_1 = self.spline_derivatives[n + 1] + + new_shape = (-1,) + (1,) * self.number_of_custom_axes + h00 = h00.reshape(new_shape) + h10 = h10.reshape(new_shape) + h01 = h01.reshape(new_shape) + h11 = h11.reshape(new_shape) + + interpolated_values = ( + h00 * p_k + h10 * delta_x * m_k + h01 * p_k_1 + h11 * delta_x * m_k_1 + ) + + return interpolated_values + + +class RadialIntegralFromFunction(RadialIntegralSplinerBase): + r"""Compute the radial integral spline points based on a provided function. + + :parameter radial_integral: Function to compute the radial integral. Function must + take ``n``, ``l``, and ``positions`` as inputs, where ``n`` and ``l`` are + integers and ``positions`` is a numpy 1-D array that contains the spline points + at which the radial integral will be evaluated. The function must return a numpy + 1-D array containing the values of the radial integral. + :parameter spline_cutoff: cutoff radius for the spline interpolation. This is also + the maximal value that can be interpolated. + :parameter max_radial: number of angular componentss + :parameter max_angular: number of radial components + :parameter radial_integral_derivative: The derivative of the radial integral taking + the same paramaters as ``radial_integral``. If it is :py:obj:`None` (default), + finite differences are used to calculate the derivative of the radial integral. + It is recommended to provide this parameter if possible. Derivatives from finite + differences can cause problems when evaluating at the edges of the domain (i.e., + at ``0`` and ``spline_cutoff``) because the function might not be defined + outside of the domain. + :parameter accuracy: accuracy of the numerical integration and the splining. + Accuracy is reached when either the mean absolute error or the mean relative + error gets below the ``accuracy`` threshold. + :parameter center_contribution: Contribution of the central atom required for LODE + calculations. The ``center_contribution`` is defined as + + .. math:: + c_n = \sqrt{4π}\int_0^\infty dr r^2 R_n(r) g(r) + + where :math:`g(r)` is the radially symmetric density function, `R_n(r)` the + radial basis function and :math:`n` the current radial channel. This should be + implemented in the child classes. + + Example + ------- + First define a ``radial_integral`` function + + >>> def radial_integral(n, ell, r): + ... return np.sin(r) + ... + + and provide this as input to the spline generator + + >>> spliner = RadialIntegralFromFunction( + ... radial_integral=radial_integral, + ... max_radial=12, + ... max_angular=9, + ... spline_cutoff=8.0, + ... ) + + Finally, we can use the ``spliner`` directly in the ``radial_integral`` section of a + calculator + + >>> from rascaline import SoapPowerSpectrum + >>> calculator = SoapPowerSpectrum( + ... cutoff=8.0, + ... max_radial=12, + ... max_angular=9, + ... center_atom_weight=1.0, + ... radial_basis=spliner.compute(), + ... atomic_gaussian_width=1.0, # ignored + ... cutoff_function={"Step": {}}, + ... ) + + The ``atomic_gaussian_width`` paramater is required by the calculator but will be + will be ignored during the feature computation. + + A more in depth example using a "rectangular" Laplacian eigenstate basis + is provided in the :ref:`example section`. + """ + + def __init__( + self, + radial_integral: Callable[[int, int, np.ndarray], np.ndarray], + spline_cutoff: float, + max_radial: int, + max_angular: int, + radial_integral_derivative: Optional[ + Callable[[int, int, np.ndarray], np.ndarray] + ] = None, + center_contribution: Optional[np.ndarray] = None, + accuracy: float = 1e-8, + ): + self.radial_integral_function = radial_integral + self.radial_integral_derivative_funcion = radial_integral_derivative + self.center_contribution = center_contribution + + super().__init__( + max_radial=max_radial, + max_angular=max_angular, + spline_cutoff=spline_cutoff, + accuracy=accuracy, + ) + + def _radial_integral(self, n: int, ell: int, positions: np.ndarray) -> np.ndarray: + return self.radial_integral_function(n, ell, positions) + + @property + def _center_contribution(self) -> Union[None, np.ndarray]: + # Test that ``len(self.center_contribution) == max_radial`` is performed by the + # calculator. + return self.center_contribution + + def _radial_integral_derivative( + self, n: int, ell: int, positions: np.ndarray + ) -> np.ndarray: + if self.radial_integral_derivative_funcion is None: + return super()._radial_integral_derivative(n, ell, positions) + else: + return self.radial_integral_derivative_funcion(n, ell, positions) diff --git a/python/rascaline/tests/splines.py b/python/rascaline/tests/splines.py deleted file mode 100644 index 02b53a6d1..000000000 --- a/python/rascaline/tests/splines.py +++ /dev/null @@ -1,99 +0,0 @@ -import numpy as np - -import rascaline - - -def sine(n, ell, r): - return np.sin(r) - - -def cosine(n, ell, r): - return np.cos(r) - - -def test_splines_with_n_spline_points(): - cutoff_radius = 8.0 - - spline_points = rascaline.generate_splines( - sine, - cosine, - max_radial=12, - max_angular=9, - cutoff_radius=cutoff_radius, - n_spline_points=1234, - ) - - # check that the first spline point is at 0 - assert spline_points[0]["position"] == 0.0 - - # check that the last spline point is at the cutoff radius - assert spline_points[-1]["position"] == 8.0 - - # ensure correct length for values representation - assert len(spline_points[52]["values"]["data"]) == (9 + 1) * 12 - - # ensure correct length for derivatives representation - assert len(spline_points[23]["derivatives"]["data"]) == (9 + 1) * 12 - - # check values at r = 0.0 - assert np.allclose( - np.array(spline_points[0]["values"]["data"]), np.zeros((9 + 1) * 12) - ) - - # check derivatives at r = 0.0 - assert np.allclose( - np.array(spline_points[0]["derivatives"]["data"]), np.ones((9 + 1) * 12) - ) - - n_spline_points = len(spline_points) - random_spline_point = 123 - random_x = random_spline_point * cutoff_radius / (n_spline_points - 1) - - # check value of a random spline point - assert np.allclose( - np.array(spline_points[random_spline_point]["values"]["data"]), - np.sin(random_x) * np.ones((9 + 1) * 12), - ) - - -def test_splines_with_accuracy(): - cutoff_radius = 8.0 - spline_points = rascaline.generate_splines( - sine, - cosine, - max_radial=12, - max_angular=9, - cutoff_radius=cutoff_radius, - ) - - # check that the first spline point is at 0 - assert spline_points[0]["position"] == 0.0 - - # check that the last spline point is at the cutoff radius - assert spline_points[-1]["position"] == 8.0 - - # ensure correct length for values representation - assert len(spline_points[52]["values"]["data"]) == (9 + 1) * 12 - - # ensure correct length for derivatives representation - assert len(spline_points[23]["derivatives"]["data"]) == (9 + 1) * 12 - - # check values at r = 0.0 - assert np.allclose( - np.array(spline_points[0]["values"]["data"]), np.zeros((9 + 1) * 12) - ) - - # check derivatives at r = 0.0 - assert np.allclose( - np.array(spline_points[0]["derivatives"]["data"]), np.ones((9 + 1) * 12) - ) - - n_spline_points = len(spline_points) - random_spline_point = 123 - random_x = random_spline_point * cutoff_radius / (n_spline_points - 1) - - # check value of a random spline point - assert np.allclose( - np.array(spline_points[random_spline_point]["values"]["data"]), - np.sin(random_x) * np.ones((9 + 1) * 12), - ) diff --git a/python/rascaline/tests/utils/splines.py b/python/rascaline/tests/utils/splines.py new file mode 100644 index 000000000..a08cc54e3 --- /dev/null +++ b/python/rascaline/tests/utils/splines.py @@ -0,0 +1,104 @@ +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_equal + +from rascaline.utils import RadialIntegralFromFunction + + +def sine(n: int, ell: int, positions: np.ndarray) -> np.ndarray: + return np.sin(positions) + + +def cosine(n: int, ell: int, positions: np.ndarray) -> np.ndarray: + return np.cos(positions) + + +@pytest.mark.parametrize("n_spline_points", [None, 1234]) +def test_splines_with_n_spline_points(n_spline_points): + spline_cutoff = 8.0 + + spliner = RadialIntegralFromFunction( + radial_integral=sine, + max_radial=12, + max_angular=9, + spline_cutoff=spline_cutoff, + radial_integral_derivative=cosine, + ) + + radial_integral = spliner.compute(n_spline_points=n_spline_points)[ + "TabulatedRadialIntegral" + ] + + # check central contribution is not added + with pytest.raises(KeyError): + radial_integral["center_contribution"] + + spline_points = radial_integral["points"] + + # check that the first spline point is at 0 + assert spline_points[0]["position"] == 0.0 + + # check that the last spline point is at the cutoff radius + assert spline_points[-1]["position"] == 8.0 + + # ensure correct length for values representation + assert len(spline_points[52]["values"]["data"]) == (9 + 1) * 12 + + # ensure correct length for derivatives representation + assert len(spline_points[23]["derivatives"]["data"]) == (9 + 1) * 12 + + # check values at r = 0.0 + assert np.allclose( + np.array(spline_points[0]["values"]["data"]), np.zeros((9 + 1) * 12) + ) + + # check derivatives at r = 0.0 + assert np.allclose( + np.array(spline_points[0]["derivatives"]["data"]), np.ones((9 + 1) * 12) + ) + + n_spline_points = len(spline_points) + random_spline_point = 123 + random_x = random_spline_point * spline_cutoff / (n_spline_points - 1) + + # check value of a random spline point + assert np.allclose( + np.array(spline_points[random_spline_point]["values"]["data"]), + np.sin(random_x) * np.ones((9 + 1) * 12), + ) + + +def test_splines_numerical_derivative(): + kwargs = { + "radial_integral": sine, + "max_radial": 12, + "max_angular": 9, + "spline_cutoff": 8.0, + } + + spliner = RadialIntegralFromFunction(**kwargs, radial_integral_derivative=cosine) + spliner_numerical = RadialIntegralFromFunction(**kwargs) + + spline_points = spliner.compute()["TabulatedRadialIntegral"]["points"] + spline_points_numerical = spliner_numerical.compute()["TabulatedRadialIntegral"][ + "points" + ] + + for s, s_num in zip(spline_points, spline_points_numerical): + assert_equal(s["values"]["data"], s_num["values"]["data"]) + assert_allclose( + s["derivatives"]["data"], s_num["derivatives"]["data"], rtol=1e-7 + ) + + +def test_splines_numerical_derivative_error(): + kwargs = { + "radial_integral": sine, + "max_radial": 12, + "max_angular": 9, + "spline_cutoff": 1e-3, + } + + match = "Numerically derivative of the radial integral can not be performed" + with pytest.raises(ValueError, match=match): + RadialIntegralFromFunction(**kwargs).compute()