From d451f41d148e2641c7969a35d66694a34e6b10bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 14 Sep 2023 10:18:37 -0400 Subject: [PATCH 01/10] Add class docstring --- src/grid/ngrid.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 src/grid/ngrid.py diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py new file mode 100644 index 000000000..83f1a388e --- /dev/null +++ b/src/grid/ngrid.py @@ -0,0 +1,20 @@ +# 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 +# -- +"""Grid for N particle functions.""" From e2afc4d4fe048bc6627a80457730bac497a0f76a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Fri, 15 Sep 2023 15:10:56 -0400 Subject: [PATCH 02/10] Add N-particle grid --- src/grid/ngrid.py | 117 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py index 83f1a388e..dbc79eac0 100644 --- a/src/grid/ngrid.py +++ b/src/grid/ngrid.py @@ -18,3 +18,120 @@ # along with this program; if not, see # -- """Grid for N particle functions.""" + +import numpy as np +from grid.basegrid import Grid +import itertools + + +class Ngrid(Grid): + r""" + Grid class for integration of N particle functions. + """ + + def __init__(self, grid_list=None, n=None, **kwargs): + r""" + Initialize n particle grid. + + Parameters + ---------- + grid_list : list of Grid + List of grids for each particle. + n : int + Number of particles. + """ + # check that grid_list is defined + if grid_list is None: + raise ValueError("The list must be specified") + # TODO check that grid_list is a list of grids + + self.grid_list = grid_list + self.n = n + + def integrate(self, callable, **call_kwargs): + r""" + Integrate callable on the N particle grid. + + Parameters + ---------- + callable : callable + Callable to integrate. It must take a list of N three dimensional tuples as argument + (one for each particle) and return a float (e.g. a function of the form + f((x1,y1,z1), (x2,y2,z2)) -> float). + call_kwargs : dict + Keyword arguments that will be passed to callable. + + Returns + ------- + float + Integral of callable. + """ + if len(self.grid_list) == 0: + raise ValueError("The list must contain at least one grid") + + if len(self.grid_list) == 1 and self.n > 1: + return self._n_integrate(self.grid_list * self.n, callable, **call_kwargs) + else: + return self._n_integrate(self.grid_list, callable, **call_kwargs) + + def _n_integrate(self, grid_list, callable, **call_kwargs): + r""" + Integrate callable on the space spanned domain union of the grids in grid_list. + + Parameters + ---------- + grid_list : list of Grid + List of grids for each particle. + callable : callable + Callable to integrate. It must take a list of N three dimensional tuples as argument + (one for each particle) and return a float (e.g. a function of the form + f((x1,y1,z1), (x2,y2,z2)) -> float). + kwcallargs : dict + Keyword arguments for callable. + + Returns + ------- + float + Integral of callable. + """ + + # if there is only one grid, perform the integration using the integrate method of the Grid + if len(grid_list) == 1: + vals = callable(grid_list[0].points, **kwcallargs) + return grid_list[0].integrate(vals) + else: + # The integration is performed by integrating the function over the last grid with all + # the other coordinates fixed for each possible combination of the other grids' points. + # + # Notes: + # ----- + # - The combination of the other grids' points is generated using a generator so that + # the memory usage is kept low. + # - The last grid is integrated using the integrate method of the Grid class. + + # generate all possible combinations for the first n-1 grids + points = itertools.product(*[grid.points for grid in grid_list[:-1]]) + # generate all corresponding weights to the generated points + weights = itertools.product(*[grid.weights for grid in grid_list[:-1]]) + # zip the points and weights together + data = zip(points, weights) + + integral = 0.0 + for i in data: + # define an auxiliar function that takes a single argument (a point of the last + # grid) but uses the other coordinates as fixed parameters i[0] and returns the + # value of the n particle function at that point (i.e. the value of the n particle + # function at the point defined by the last grid point and the other coordinates + # fixed by i[0]) + aux_func = lambda x: callable(*i[0], x, **kwcallargs) + + # calculate the value of the n particle function at each point of the last grid + vals = aux_func(grid_list[-1].points) + + # Integrate the function over the last grid with all the other coordinates fixed. + # The result is multiplied by the product of the weights corresponding to the other + # grids' points (stored in i[1]). + # This is equivalent to integrating the n particle function over the coordinates of + # the last particle with the other coordinates fixed. + integral += grid_list[-1].integrate(vals) * np.prod(i[1]) + return integral From 412879bd63e083db74607e09c24a0f04f1d6dc11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 16 Sep 2023 00:01:43 -0400 Subject: [PATCH 03/10] Improve documentation for Ngrid --- src/grid/ngrid.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py index dbc79eac0..2da79a6b6 100644 --- a/src/grid/ngrid.py +++ b/src/grid/ngrid.py @@ -27,23 +27,37 @@ class Ngrid(Grid): r""" Grid class for integration of N particle functions. + + The function to integrate must be a function of the form + f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> float where (xi, yi, zi) are the coordinates + of the i-th particle and n is the number of particles. Internally, one Grid is created for + each particle and the integration is performed by integrating the function over the space spanned + by the domain union of all the grids. """ def __init__(self, grid_list=None, n=None, **kwargs): r""" Initialize n particle grid. + At least one grid must be specified. If only one grid is specified, and a value for n bigger + than one is specified, the same grid will copied n times and one grid will used for each + particle. If more than one grid is specified, n will be ignored. In all cases, the function + to integrate must be a function of the form f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> + float where (xi, yi, zi) are the coordinates of the i-th particle and the function must + depend on a number of particles equal to the number of grids. + Parameters ---------- grid_list : list of Grid - List of grids for each particle. + List of grids, one Grid for each particle. n : int Number of particles. """ # check that grid_list is defined if grid_list is None: raise ValueError("The list must be specified") - # TODO check that grid_list is a list of grids + if not all(isinstance(grid, Grid) for grid in grid_list): + raise ValueError("The Grid list must contain only Grid objects") self.grid_list = grid_list self.n = n From 1a28e756bb4c0de1775c17c745a2d40c6c9beb9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Mon, 18 Sep 2023 12:55:28 -0400 Subject: [PATCH 04/10] Add checks and minos name fixes --- src/grid/ngrid.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py index 2da79a6b6..ee661d881 100644 --- a/src/grid/ngrid.py +++ b/src/grid/ngrid.py @@ -56,9 +56,26 @@ def __init__(self, grid_list=None, n=None, **kwargs): # check that grid_list is defined if grid_list is None: raise ValueError("The list must be specified") + + # check that grid_list is not empty + if len(grid_list) == 0: + raise ValueError("The list must contain at least one grid") + + # check that grid_list contains only Grid objects if not all(isinstance(grid, Grid) for grid in grid_list): raise ValueError("The Grid list must contain only Grid objects") + if n is not None: + # check that n is non negative + if n < 0: + raise ValueError("n must be non negative") + # check that for n > 1, the number of grids is equal to n or 1 + if len(grid_list) > 1 and len(grid_list) != n: + raise ValueError( + "Conflicting values for n and the number of grids. \n" + "If n is specified, the number of grids must be equal to n or 1." + ) + self.grid_list = grid_list self.n = n @@ -80,10 +97,11 @@ def integrate(self, callable, **call_kwargs): float Integral of callable. """ + # check that grid_list is not empty if len(self.grid_list) == 0: raise ValueError("The list must contain at least one grid") - if len(self.grid_list) == 1 and self.n > 1: + if len(self.grid_list) == 1 and self.n is not None and self.n > 1: return self._n_integrate(self.grid_list * self.n, callable, **call_kwargs) else: return self._n_integrate(self.grid_list, callable, **call_kwargs) @@ -100,7 +118,7 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): Callable to integrate. It must take a list of N three dimensional tuples as argument (one for each particle) and return a float (e.g. a function of the form f((x1,y1,z1), (x2,y2,z2)) -> float). - kwcallargs : dict + call_kwargs : dict Keyword arguments for callable. Returns @@ -111,7 +129,7 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): # if there is only one grid, perform the integration using the integrate method of the Grid if len(grid_list) == 1: - vals = callable(grid_list[0].points, **kwcallargs) + vals = callable(grid_list[0].points, **call_kwargs) return grid_list[0].integrate(vals) else: # The integration is performed by integrating the function over the last grid with all @@ -137,7 +155,7 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): # value of the n particle function at that point (i.e. the value of the n particle # function at the point defined by the last grid point and the other coordinates # fixed by i[0]) - aux_func = lambda x: callable(*i[0], x, **kwcallargs) + aux_func = lambda x: callable(*i[0], x, **call_kwargs) # calculate the value of the n particle function at each point of the last grid vals = aux_func(grid_list[-1].points) From 9e5eebba9a3e28691b3a074a255a05f8251c53da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Mon, 18 Sep 2023 12:57:17 -0400 Subject: [PATCH 05/10] Add test for Ngrid init --- src/grid/tests/test_ngrid.py | 64 ++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 src/grid/tests/test_ngrid.py diff --git a/src/grid/tests/test_ngrid.py b/src/grid/tests/test_ngrid.py new file mode 100644 index 000000000..d703b000c --- /dev/null +++ b/src/grid/tests/test_ngrid.py @@ -0,0 +1,64 @@ +# 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 +# -- +"""Ngrid tests file.""" + + +from unittest import TestCase +from grid.onedgrid import UniformInteger +from grid.rtransform import LinearInfiniteRTransform +from grid.ngrid import Ngrid + +import numpy as np +from numpy.testing import assert_allclose + + +class TestNgrid(TestCase): + """Ngrid tests class.""" + + def setUp(self): + """Set up the test.""" + # define the number of points + n = 500 + # create a linear grid with n points + self.linear_grid = UniformInteger(n) + # transform its boundaries to 0 and 1 + self.linear_grid = LinearInfiniteRTransform(rmin=1.0e-4, rmax=1.0).transform_1d_grid( + self.linear_grid + ) + # create a 3D grid with n equally spaced points between 0 and 1 along each axis + self.ngrid = Ngrid([self.linear_grid], 3) + + def test_init_raises(self): + """Assert that the init raises the correct error.""" + # case 1: the grid list is not given + with self.assertRaises(ValueError): + Ngrid(grid_list=None) + # case 2: the grid list is empty + with self.assertRaises(ValueError): + Ngrid(grid_list=[]) + # case 3: the grid list is not a list of Grid + with self.assertRaises(ValueError): + Ngrid(grid_list=[self.linear_grid, 1], n=None) + # case 4: n is negative + with self.assertRaises(ValueError): + Ngrid(grid_list=[self.linear_grid], n=-1) + # case 5: n and the grid list have different lengths + with self.assertRaises(ValueError): + Ngrid(grid_list=[self.linear_grid] * 3, n=2) \ No newline at end of file From e4154c623d9b4ddeedcafa33a9a0dab3985d7ef3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Mon, 18 Sep 2023 12:58:33 -0400 Subject: [PATCH 06/10] Add tests to Ngrid --- src/grid/tests/test_ngrid.py | 56 +++++++++++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/src/grid/tests/test_ngrid.py b/src/grid/tests/test_ngrid.py index d703b000c..da7edb9f9 100644 --- a/src/grid/tests/test_ngrid.py +++ b/src/grid/tests/test_ngrid.py @@ -61,4 +61,58 @@ def test_init_raises(self): Ngrid(grid_list=[self.linear_grid], n=-1) # case 5: n and the grid list have different lengths with self.assertRaises(ValueError): - Ngrid(grid_list=[self.linear_grid] * 3, n=2) \ No newline at end of file + Ngrid(grid_list=[self.linear_grid] * 3, n=2) + + def test_init(self): + """Assert that the init works as expected.""" + # case 1: the grid list is given and n is None + ngrid = Ngrid(grid_list=[self.linear_grid, self.linear_grid, self.linear_grid]) + self.assertEqual(len(ngrid.grid_list), 3) + self.assertEqual(ngrid.n, None) + + # case 2: the grid list is given (length 1) and n is not None + ngrid = Ngrid(grid_list=[self.linear_grid], n=3) + self.assertEqual(len(ngrid.grid_list), 1) + self.assertEqual(ngrid.n, 3) + + def test_single_grid_integration(self): + """Assert that the integration works as expected for a single grid.""" + + # define a function to integrate (x**2) + def f(x): + return x**2 + + # define a Ngrid with only one grid + ngrid = Ngrid(grid_list=[self.linear_grid]) + # integrate it + result = ngrid.integrate(f) + # check that the result is correct + self.assertAlmostEqual(result, 1.0 / 3.0, places=2) + + def test_2_grid_integration(self): + """Assert that the integration works as expected for two grids.""" + + # define a function to integrate (x**2+y**2) + def f(x, y): + return x**2 + y**2 + + # define a Ngrid with two grids + ngrid = Ngrid(grid_list=[self.linear_grid], n=2) + # integrate it + result = ngrid.integrate(f) + # check that the result is correct + self.assertAlmostEqual(result, 2.0 / 3.0, places=2) + + def test_3_grid_integration(self): + """Assert that the integration works as expected for three grids.""" + + # define a function to integrate (x**2+y**2+z**2) + def f(x, y, z): + return x * y * z + + # define a Ngrid with three grids + ngrid = Ngrid(grid_list=[self.linear_grid], n=3) + # integrate it + result = ngrid.integrate(f) + # check that the result is correct + self.assertAlmostEqual(result, 1.0 / 8.0, places=2) From 05a5b40d25fe06b23c92eb97c083dfa9247031c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Mon, 18 Sep 2023 12:58:55 -0400 Subject: [PATCH 07/10] Add ngrid to init file imports --- src/grid/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/grid/__init__.py b/src/grid/__init__.py index 6d020590e..e3e607fd9 100644 --- a/src/grid/__init__.py +++ b/src/grid/__init__.py @@ -34,3 +34,4 @@ from grid.onedgrid import * from grid.periodicgrid import * from grid.rtransform import * +from grid.ngrid import * From f565b5dcb1e494044f718d223258e6368b78417f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 23 May 2024 21:33:25 -0400 Subject: [PATCH 08/10] Modify input handling of Ngrids The function to integrate now has few requirements. It has to be a function of n variables where n is the number of grids, all input parameters must have the same dimensions as the corresponding grid points. Added two tests too. --- src/grid/ngrid.py | 47 +++++++------ src/grid/tests/test_ngrid.py | 123 ++++++++++++++++++++++++++++++++++- 2 files changed, 146 insertions(+), 24 deletions(-) diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py index ee661d881..46d85a437 100644 --- a/src/grid/ngrid.py +++ b/src/grid/ngrid.py @@ -22,16 +22,19 @@ import numpy as np from grid.basegrid import Grid import itertools +from numbers import Number class Ngrid(Grid): r""" Grid class for integration of N particle functions. - The function to integrate must be a function of the form + The function to integrate must be a function with all arguments with the same dimension as the + grid points. The function must depend on a number of particles equal to the number of grids. For + example, a function of the form f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> float where (xi, yi, zi) are the coordinates - of the i-th particle and n is the number of particles. Internally, one Grid is created for - each particle and the integration is performed by integrating the function over the space spanned + of the i-th particle and n is the number of particles. Internally, one Grid is created for each + particle and the integration is performed by integrating the function over the space spanned by the domain union of all the grids. """ @@ -41,10 +44,12 @@ def __init__(self, grid_list=None, n=None, **kwargs): At least one grid must be specified. If only one grid is specified, and a value for n bigger than one is specified, the same grid will copied n times and one grid will used for each - particle. If more than one grid is specified, n will be ignored. In all cases, the function - to integrate must be a function of the form f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> - float where (xi, yi, zi) are the coordinates of the i-th particle and the function must - depend on a number of particles equal to the number of grids. + particle. If more than one grid is specified, n will be ignored. In all cases, The function + to integrate must be a function with all arguments with the same dimension as the grid + points and must depend on a number of particles equal to the number of grids. For example, a + function of the form + f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> float where (xi, yi, zi) are the coordinates + of the i-th particle and n is the number of particles. Parameters ---------- @@ -86,9 +91,9 @@ def integrate(self, callable, **call_kwargs): Parameters ---------- callable : callable - Callable to integrate. It must take a list of N three dimensional tuples as argument - (one for each particle) and return a float (e.g. a function of the form - f((x1,y1,z1), (x2,y2,z2)) -> float). + Callable to integrate. It must take a list of arguments (one for each particle) with + the same dimension as the grid points and return a float (e.g. a function of the form + f([x1,y1,z1], [x2,y2,z2]) -> float). call_kwargs : dict Keyword arguments that will be passed to callable. @@ -115,9 +120,9 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): grid_list : list of Grid List of grids for each particle. callable : callable - Callable to integrate. It must take a list of N three dimensional tuples as argument - (one for each particle) and return a float (e.g. a function of the form - f((x1,y1,z1), (x2,y2,z2)) -> float). + Callable to integrate. It must take a list of arguments (one for each particle) with + the same dimension as the grid points and return a float (e.g. a function of the form + f([x1,y1,z1], [x2,y2,z2]) -> float). call_kwargs : dict Keyword arguments for callable. @@ -142,28 +147,28 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): # - The last grid is integrated using the integrate method of the Grid class. # generate all possible combinations for the first n-1 grids - points = itertools.product(*[grid.points for grid in grid_list[:-1]]) - # generate all corresponding weights to the generated points - weights = itertools.product(*[grid.weights for grid in grid_list[:-1]]) - # zip the points and weights together - data = zip(points, weights) + data = itertools.product(*[zip(grid.points, grid.weights) for grid in grid_list[:-1]]) integral = 0.0 for i in data: + to_point = lambda x: np.array([[x]]) if isinstance(x, Number) else x[None, :] + # extract points (convert to array, add dim) and corresponding weights combinations + points_comb = (to_point(j[0]) for j in i) + weights_comb = np.array([j[1] for j in i]) # define an auxiliar function that takes a single argument (a point of the last # grid) but uses the other coordinates as fixed parameters i[0] and returns the # value of the n particle function at that point (i.e. the value of the n particle # function at the point defined by the last grid point and the other coordinates # fixed by i[0]) - aux_func = lambda x: callable(*i[0], x, **call_kwargs) + aux_func = lambda x: callable(*points_comb, x, **call_kwargs) # calculate the value of the n particle function at each point of the last grid - vals = aux_func(grid_list[-1].points) + vals = aux_func(grid_list[-1].points).flatten() # Integrate the function over the last grid with all the other coordinates fixed. # The result is multiplied by the product of the weights corresponding to the other # grids' points (stored in i[1]). # This is equivalent to integrating the n particle function over the coordinates of # the last particle with the other coordinates fixed. - integral += grid_list[-1].integrate(vals) * np.prod(i[1]) + integral += grid_list[-1].integrate(vals) * np.prod(weights_comb) return integral diff --git a/src/grid/tests/test_ngrid.py b/src/grid/tests/test_ngrid.py index da7edb9f9..acba2b50c 100644 --- a/src/grid/tests/test_ngrid.py +++ b/src/grid/tests/test_ngrid.py @@ -21,15 +21,18 @@ from unittest import TestCase -from grid.onedgrid import UniformInteger -from grid.rtransform import LinearInfiniteRTransform +from grid.onedgrid import UniformInteger, GaussLegendre +from grid.rtransform import BeckeRTransform, LinearInfiniteRTransform +from grid.atomgrid import AtomGrid from grid.ngrid import Ngrid +import pytest +import itertools import numpy as np from numpy.testing import assert_allclose -class TestNgrid(TestCase): +class TestNgrid_linear(TestCase): """Ngrid tests class.""" def setUp(self): @@ -116,3 +119,117 @@ def f(x, y, z): result = ngrid.integrate(f) # check that the result is correct self.assertAlmostEqual(result, 1.0 / 8.0, places=2) + + +class TestNgrid_atom(TestCase): + """Ngrid tests class.""" + + def setUp(self): + """Set up atomgrid to use in the tests.""" + # construct a radial grid with 150 points + oned = GaussLegendre(npoints=150) + rgrid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned) + self.center1 = np.array([50.0, 0.0, -0.5], float) + self.center2 = np.array([0.0, 50.0, 0.5], float) + self.center3 = np.array([0.0, -0.5, 50.0], float) + self.atgrid1 = AtomGrid(rgrid, degrees=[22], center=self.center1) + self.atgrid2 = AtomGrid(rgrid, degrees=[22], center=self.center2) + self.atgrid3 = AtomGrid(rgrid, degrees=[22], center=self.center3) + + def make_gaussian(self, center, height, width): + """Make a gaussian function. + + Parameters + ---------- + center : np.ndarray + The center of the gaussian. + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + + def f(x): + r = np.linalg.norm(center - x, axis=1) + return height * np.exp(-(r**2) / (2 * width**2)) + + return f + + def gaussian_integral(self, height, width): + """Calculate the integral of a gaussian function. + + Parameters + ---------- + center : np.ndarray + The center of the gaussian. + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + return height * (width * np.sqrt(2 * np.pi)) ** 3 + + def test_single_grid_integration(self): + """Assert that the integration works as expected for a single grid.""" + + height, width = 1.0, 2.0 + + # define a function to integrate + g = self.make_gaussian(self.center1, height, width) + + # define a Ngrid with only one atom grid + ngrid = Ngrid(grid_list=[self.atgrid1]) + # integrate it + result = ngrid.integrate(g) + ref_val = self.gaussian_integral(height, width) + # check that the result is correct + self.assertAlmostEqual(result, ref_val, places=6) + + def test_2_grid_integration(self): + """Assert that the integration works as expected for two grids.""" + + height1, width1 = 1.0, 2.0 + height2, width2 = 0.7, 1.5 + + # define a function to integrate + g1 = self.make_gaussian(self.center1, height1, width1) + g2 = self.make_gaussian(self.center2, height2, width2) + + # function is product of two gaussians + func = lambda x, y: g1(x) * g2(y) + + # define a Ngrid with two grids + ngrid = Ngrid(grid_list=[self.atgrid1, self.atgrid2]) + # integrate it and compare with reference value + result = ngrid.integrate(func) + ref_val = self.gaussian_integral(height1, width1) * self.gaussian_integral(height2, width2) + + # check that the result is correct + self.assertAlmostEqual(result, ref_val, places=6) + + @pytest.mark.skip(reason="This is to slow.") + def test_3_grid_integration(self): + """Assert that the integration works as expected for three grids.""" + + height1, width1 = 1.0, 2.0 + height2, width2 = 0.7, 1.5 + height3, width3 = 0.5, 1.0 + + # define a function to integrate + g1 = self.make_gaussian(self.center1, height1, width1) + g2 = self.make_gaussian(self.center2, height2, width2) + g3 = self.make_gaussian(self.center3, height3, width3) + + # function is product of two gaussians + func = lambda x, y, z: g1(x) * g2(y) * g3(z) + + # define a Ngrid with two grids + ngrid = Ngrid(grid_list=[self.atgrid1, self.atgrid2, self.atgrid3]) + # integrate it and compare with reference value + result = ngrid.integrate(func) + ref_val = ( + self.gaussian_integral(height1, width1) + * self.gaussian_integral(height2, width2) + * self.gaussian_integral(height3, width3) + ) + + # check that the result is correct + self.assertAlmostEqual(result, ref_val, places=6) From 4773323412ffb9a77aeb4aab623cb7491fbf97cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Fri, 24 May 2024 00:27:50 -0400 Subject: [PATCH 09/10] Add mixed grids test --- src/grid/tests/test_ngrid.py | 104 ++++++++++++++++++++++++++++++++++- 1 file changed, 102 insertions(+), 2 deletions(-) diff --git a/src/grid/tests/test_ngrid.py b/src/grid/tests/test_ngrid.py index acba2b50c..edb75459e 100644 --- a/src/grid/tests/test_ngrid.py +++ b/src/grid/tests/test_ngrid.py @@ -159,8 +159,6 @@ def gaussian_integral(self, height, width): Parameters ---------- - center : np.ndarray - The center of the gaussian. height : float The height of the gaussian. width : float @@ -233,3 +231,105 @@ def test_3_grid_integration(self): # check that the result is correct self.assertAlmostEqual(result, ref_val, places=6) + + +class TestNgrid_mixed(TestCase): + """Ngrid tests class.""" + + def setUp(self): + """Set up atomgrid to use in the tests.""" + # construct a radial grid with 150 points + oned = GaussLegendre(npoints=150) + rgrid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned) + # the radial grid for points > 0, integrates half of the gaussian + self.center1d = 0.0 + self.center3d = np.array([0.0, 50.0, 0.5], float) + + self.onedgrid = rgrid + self.atgrid1 = AtomGrid(rgrid, degrees=[22], center=self.center3d) + + def make_gaussian_1d(self, center, height, width): + """Make a gaussian function. + + Parameters + ---------- + center : float + The center of the gaussian. + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + + def f(r): + return height * np.exp(-(r**2) / (2 * width**2)) + + return f + + def make_gaussian_3d(self, center, height, width): + """Make a gaussian function. + + Parameters + ---------- + center : np.ndarray + The center of the gaussian. + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + + def f(x): + r = np.linalg.norm(center - x, axis=1) + return height * np.exp(-(r**2) / (2 * width**2)) + + return f + + def gaussian_integral_1d(self, height, width): + """Calculate the integral of a gaussian function. + + Parameters + ---------- + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + return height * (width * np.sqrt(2 * np.pi)) + + def gaussian_integral_3d(self, height, width): + """Calculate the integral of a gaussian function. + + Parameters + ---------- + height : float + The height of the gaussian. + width : float + The width of the gaussian.""" + return height * (width * np.sqrt(2 * np.pi)) ** 3 + + def test_single_mixed_1d_3d(self): + """Assert that the integration works as expected for a single grid.""" + + height_1d, width_1d = 1.0, 2.0 + height_3d, width_3d = 3.1, 1.5 + + # define a function to integrate + g1d = self.make_gaussian_1d(self.center1d, height_1d, width_1d) + g3d = self.make_gaussian_3d(self.center3d, height_3d, width_3d) + + # define a function to integrate + g1 = self.make_gaussian_1d(self.center1d, height_1d, width_1d) + g2 = self.make_gaussian_3d(self.center3d, height_3d, width_3d) + + # function is product of two gaussians + func = lambda x, y: g1(x) * g2(y) + + # define a Ngrid with only one atom grid + ngrid = Ngrid(grid_list=[self.onedgrid, self.atgrid1]) + # integrate it + result = ngrid.integrate(func) + ref_val = ( + self.gaussian_integral_1d(height_1d, width_1d) + / 2 + * self.gaussian_integral_3d(height_3d, width_3d) + ) + # check that the result is correct + self.assertAlmostEqual(result, ref_val, places=6) From 06cfd108491841eb5ad3f77627338b0c1174f652 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Mon, 10 Feb 2025 15:46:02 -0500 Subject: [PATCH 10/10] Improve docstring --- src/grid/ngrid.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/grid/ngrid.py b/src/grid/ngrid.py index 46d85a437..78afb9aa3 100644 --- a/src/grid/ngrid.py +++ b/src/grid/ngrid.py @@ -27,15 +27,17 @@ class Ngrid(Grid): r""" - Grid class for integration of N particle functions. - - The function to integrate must be a function with all arguments with the same dimension as the - grid points. The function must depend on a number of particles equal to the number of grids. For - example, a function of the form - f((x1,y1,z1), (x2,y2,z2), ..., (xn, yn, zn)) -> float where (xi, yi, zi) are the coordinates - of the i-th particle and n is the number of particles. Internally, one Grid is created for each - particle and the integration is performed by integrating the function over the space spanned - by the domain union of all the grids. + Grid class for integration of N argument functions. + + This class is used for integrating functions of N arguments. + + ..math:: + \idotsint f(x_1, x_2, ..., x_N) dx_1 dx_2 ... dx_N + + The function to integrate must have all arguments :math:`\{x_i\}` with the same dimension as the + points of the corresponding grids (i.e. each of the arguments corresponds to a different grid). + For example for a function of the form f((x1,y1,z1), (x2,y2)) -> float the first argument must + be described by a 3D grid and the second argument by a 2D grid. """ def __init__(self, grid_list=None, n=None, **kwargs): @@ -151,6 +153,7 @@ def _n_integrate(self, grid_list, callable, **call_kwargs): integral = 0.0 for i in data: + # Add a dimension to the point (two if it is a number) to_point = lambda x: np.array([[x]]) if isinstance(x, Number) else x[None, :] # extract points (convert to array, add dim) and corresponding weights combinations points_comb = (to_point(j[0]) for j in i)