diff --git a/README.md b/README.md index dc45951..e281552 100644 --- a/README.md +++ b/README.md @@ -56,6 +56,10 @@ number of factors: - Optimality criteria (``A``, ``C``, ``D``, ``E``, ``G``, ``I``, ``S``, ``T``, ``V``) - Search algorithms (``Sequential (Dykstra)``, ``Simple Exchange (Wynn-Mitchell)``, ``Fedorov``, ``Modified Fedorov``, ``DETMAX``) +- **Sparse Grid Designs** + - Sparse Grid Design (``doe_sparse_grid``) + - Sparse Grid Dimension (``sparse_grid_dimension``) + See [Documentation](https://pydoe3.readthedocs.io). Installation @@ -109,3 +113,4 @@ References - [Taguchi designs](http://en.wikipedia.org/wiki/Taguchi_methods) - [Generalized Subset Designs](https://doi.org/10.1021/acs.analchem.7b00506) - [Optimal experimental design](https://en.wikipedia.org/wiki/Optimal_experimental_design) +- [Sparse grid](https://en.wikipedia.org/wiki/Sparse_grid) diff --git a/doc/index.rst b/doc/index.rst index 386e6e5..5bd9da8 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -52,7 +52,9 @@ number of factors: - :ref:`Randomized Designs ` #. :ref:`Latin-Hypercube ` (``lhs``) + #. :ref:`Random K-Means ` (``random_k_means``) + #. :ref:`Random Uniform ` (``random_uniform``) - :ref:`Low-Discrepancy Sequences ` @@ -72,6 +74,7 @@ number of factors: - :ref:`Sampling Designs ` #. :ref:`Morris Method ` (``morris_sampling``) + #. :ref:`Saltelli Sampling ` (``saltelli_sampling``) - :ref:`Taguchi Designs ` @@ -81,9 +84,17 @@ number of factors: - :ref:`Optimal Designs ` #. Advanced optimal design algorithms (``optimal_design``) + #. Optimality criteria (``A``, ``C``, ``D``, ``E``, ``G``, ``I``, ``S``, ``T``, ``V``) + #. Search algorithms (``Sequential (Dykstra)``, ``Simple Exchange (Wynn-Mitchell)``, ``Fedorov``, ``Modified Fedorov``, ``DETMAX``) +- :ref:`Sparse Grid Designs ` + + #. :ref:`Sparse Grid Design ` (``doe_sparse_grid``) + + #. :ref:`Sparse Grid Dimension ` (``sparse_grid_dimension``) + Requirements ============ diff --git a/doc/index_TOC.rst b/doc/index_TOC.rst index 586b514..02bf92c 100644 --- a/doc/index_TOC.rst +++ b/doc/index_TOC.rst @@ -12,3 +12,4 @@ Table of Contents Sampling Designs Taguchi Designs Optimal Designs + Sparse Grid Designs diff --git a/doc/randomized.rst b/doc/randomized.rst index e3675af..4f0781d 100644 --- a/doc/randomized.rst +++ b/doc/randomized.rst @@ -17,7 +17,7 @@ be described: All available designs can be accessed after a simple import statement:: >>> from pyDOE3 import * - + .. index:: Latin-Hypercube .. _latin_hypercube: @@ -38,7 +38,7 @@ where generate for each factor (default: n) * **criterion**: a string that tells ``lhs`` how to sample the points (default: None, which simply randomizes the points within the intervals): - + - "center" or "c": center the points within the sampling intervals - "maximin" or "m": maximize the minimum distance between points, but place the point in a randomized location within its interval @@ -48,7 +48,7 @@ where - "lhsmu" : Latin hypercube with multifimensional Uniformity. Correlation between variable can be enforced by setting a valid correlation matrix. Description of the algorithm can be found in `Deutsch and Deutsch`_. - + The output design scales all the variable ranges from zero to one which can then be transformed as the user wishes (like to a specific statistical distribution using the `scipy.stats.distributions`_ ``ppf`` (inverse @@ -121,7 +121,7 @@ distributed with means = [1, 2, 3, 4] and standard deviations = [0.1, [ 0.91999246, 1.50179698, 2.70669743, 3.7826346 ], [ 0.97030478, 1.99322045, 3.178122 , 4.04955409], [ 1.12124679, 1.22454846, 4.52414072, 3.8707982 ]]) - + .. note:: Methods for "space-filling" designs and "orthogonal" designs are in the works, so stay tuned! However, simply increasing the samples diff --git a/doc/rsm.rst b/doc/rsm.rst index 26d1b2b..cff8d2f 100644 --- a/doc/rsm.rst +++ b/doc/rsm.rst @@ -109,7 +109,7 @@ where - ``alpha`` is either "orthogonal" (or "o", default) or "rotatable" (or "r") - + - ``face`` is either "circumscribed" (or "ccc", default), "inscribed" (or "cci"), or "faced" (or "ccf"). diff --git a/doc/sparse_grids.rst b/doc/sparse_grids.rst new file mode 100644 index 0000000..57ee770 --- /dev/null +++ b/doc/sparse_grids.rst @@ -0,0 +1,207 @@ +.. index:: Sparse Grids + +.. _sparse_grids: + +================================================================================ +Sparse Grid Designs +================================================================================ + +The ``pyDOE3`` module provides sparse grid construction using **Smolyak's construction** [Smolyak1963]_ +for generating experimental designs with hierarchical grid structures that maintain good +space-filling properties while requiring significantly fewer points than traditional +full grid approaches in high-dimensional spaces. + +.. hint:: + All sparse grid functions are available with:: + + >>> from pyDOE3 import doe_sparse_grid, sparse_grid_dimension + +Overview +======== + +This implementation is based on the MATLAB Sparse Grid Interpolation Toolbox by +Andreas Klimke [Klimke2005]_ and provides exact compatibility with MATLAB spinterp's spdim function. + +Sparse grids use **Smolyak's construction** to overcome the curse of dimensionality: +while a full grid in :math:`d` dimensions with :math:`n` points per dimension requires +:math:`n^d` total points, sparse grids require significantly fewer points. + +**Grid Types:** + - **Clenshaw-Curtis**: Nested grids based on Chebyshev polynomials (recommended) + - **Chebyshev**: Points at Chebyshev polynomial extrema + - **Gauss-Patterson**: Quadrature-based points + +.. index:: Sparse Grid Design + +.. _doe_sparse_grid: + +Sparse Grid Design (``doe_sparse_grid``) +================================================ + +The main function for generating sparse grid designs using Smolyak's construction. +This implementation exactly matches MATLAB spinterp's theoretical point counts. + +**Syntax**:: + + >>> doe_sparse_grid(n_level, n_factors, grid_type='clenshaw_curtis') + +**Parameters**: + +- ``n_level`` : int + Sparse grid level. Higher levels provide more points and better accuracy. + Level 0 gives a single center point, higher levels add structured points. + +- ``n_factors`` : int + Number of factors/dimensions in the design space. + +- ``grid_type`` : {'clenshaw_curtis', 'chebyshev', 'gauss_patterson'}, default 'clenshaw_curtis' + Type of 1D grid points to use: + + - ``'clenshaw_curtis'``: Nested Clenshaw-Curtis points (recommended) + - ``'chebyshev'``: Chebyshev polynomial extrema points + - ``'gauss_patterson'``: Gauss-Patterson quadrature points + +**Returns**: + +- ``design`` : ndarray of shape (n_points, n_factors) + Sparse grid design points in the unit hypercube [0, 1]^n_factors. + +**Examples**:: + + >>> import numpy as np + >>> from pyDOE3 import doe_sparse_grid + + >>> # Basic 2D sparse grid + >>> design = doe_sparse_grid(n_level=3, n_factors=2) + >>> print(f"Generated {len(design)} points in 2D") + Generated 29 points in 2D + + >>> # High-dimensional sparse grid + >>> design = doe_sparse_grid(n_level=4, n_factors=4) + >>> print(f"4D design with {len(design)} points") + 4D design with 177 points + + >>> # Chebyshev sparse grid + >>> design = doe_sparse_grid(n_level=2, n_factors=3, grid_type='chebyshev') + >>> print(f"Chebyshev grid: {design.shape}") + Chebyshev grid: (25, 3) + +.. note:: + The point counts follow exact polynomial formulas from Schreiber (2000) that match + MATLAB spinterp's spdim function: + + - Level 0: 1 point (center) + - Level 1: 2*d + 1 points + - Level 2: 2*d² + 2*d + 1 points + - Higher levels: polynomial growth in dimension + +.. index:: Sparse Grid Dimension + +.. _sparse_grid_dimension: + +Sparse Grid Dimension (``sparse_grid_dimension``) +================================================= + +Returns the expected number of points in a sparse grid without generating the +actual points. This is useful for planning and memory estimation. + +**Syntax**:: + + >>> sparse_grid_dimension(n_level, n_factors) + +- ``n_level``: Sparse grid level (integer ≥ 0) +- ``n_factors``: Number of factors/dimensions (integer ≥ 1) + +**Returns**: Integer number of points that would be generated + +**Example**:: + + >>> # Check point count before generation + >>> point_count = sparse_grid_dimension(n_level=5, n_factors=8) + >>> print(f"Level 5, 8D grid will have {point_count} points") + + >>> # Compare different levels + >>> for level in range(1, 6): + ... count = sparse_grid_dimension(level, 4) + ... print(f"Level {level}: {count} points") + + + +Mathematical Background +======================= + +Smolyak's Construction +---------------------- + +Sparse grids are constructed using Smolyak's formula [Smolyak1963]_, which combines univariate +interpolation rules. For a multivariate function :math:`f`, the sparse grid +interpolation operator is: + +.. math:: + + \mathcal{A}^d_{n} f = \sum_{|\mathbf{i}|_1 \leq n+d-1} (-1)^{n+d-1-|\mathbf{i}|_1} + \binom{d-1}{n+d-1-|\mathbf{i}|_1} \bigotimes_{j=1}^d \mathcal{U}^{i_j} + +where :math:`\mathbf{i} = (i_1, \ldots, i_d)` is a multi-index and :math:`|\mathbf{i}|_1 = i_1 + \cdots + i_d`. + +Point Count Formula +------------------- + +For sparse grid level :math:`n` and dimension :math:`d`: + +.. math:: + + N(n,d) = \sum_{k=0}^{n} \binom{n-k+d-1}{d-1} \cdot 2^k + +Grid Types +---------- + +- **Clenshaw-Curtis**: Nested grids based on Chebyshev polynomials (recommended) +- **Chebyshev**: Points at Chebyshev polynomial extrema +- **Gauss-Patterson**: Quadrature-based nested points + +For detailed mathematical exposition, see the `SPINTERP documentation +`_. + +Example Usage +============= + +Generate a sparse grid design for 3 factors at level 4:: + + >>> import numpy as np + >>> from pyDOE3 import doe_sparse_grid, sparse_grid_dimension + >>> + >>> # Check point count first + >>> n_points = sparse_grid_dimension(n_level=4, n_factors=3) + >>> print(f"Expected points: {n_points}") + >>> + >>> # Generate sparse grid + >>> design = doe_sparse_grid(n_level=4, n_factors=3) + >>> print(f"Generated: {design.shape}") + Generated: (177, 3) + + + +References +========== + +.. [Genz1987] Genz, A. (1987). A package for testing multiple integration subroutines. + In P. Keast & G. Fairweather (Eds.), *Numerical Integration: Recent Developments, + Software and Applications* (pp. 337-340). Reidel. ISBN: 9027725144. https://doi.org/10.1007/978-94-009-3889-2_33 + +.. [Klimke2005] Klimke, A., & Wohlmuth, B. (2005). Algorithm 847: SPINTERP: Piecewise + multilinear hierarchical sparse grid interpolation in MATLAB. *ACM Transactions on + Mathematical Software*, 31(4), 561-579. https://doi.org/10.1145/1114268.1114275 + +.. [Klimke2006] Klimke, A. (2006). *SPINTERP V2.1: Piecewise multilinear hierarchical + sparse grid interpolation in MATLAB: Documentation*. + +.. [Smolyak1963] Smolyak, S. (1963). Quadrature and interpolation formulas for tensor + products of certain classes of functions. *Doklady Akademii Nauk SSSR*, 4, 240-243. + +**Original MATLAB Documentation:** + +- `SPINTERP Toolbox `_ +- `Mathematical Details `_ +- `Implementation Guide `_ +- `Function Reference `_ diff --git a/pyDOE3/__init__.py b/pyDOE3/__init__.py index bf98149..edf92ca 100644 --- a/pyDOE3/__init__.py +++ b/pyDOE3/__init__.py @@ -56,6 +56,10 @@ from pyDOE3.doe_vanilla_morris import morris_sampling from pyDOE3.doe_saltelli import saltelli_sampling from pyDOE3.utils import scale_samples +from pyDOE3.doe_sparse_grid import ( + doe_sparse_grid, + sparse_grid_dimension, +) __all__ = [ @@ -91,6 +95,8 @@ "morris_sampling", "saltelli_sampling", "scale_samples", + "doe_sparse_grid", + "sparse_grid_dimension", ] from ._version import __version__ # pyright: ignore[reportMissingImports] # noqa diff --git a/pyDOE3/doe_sparse_grid.py b/pyDOE3/doe_sparse_grid.py new file mode 100644 index 0000000..7820225 --- /dev/null +++ b/pyDOE3/doe_sparse_grid.py @@ -0,0 +1,213 @@ +""" +Sparse Grid Design of Experiments +================================= + +This module implements sparse grid designs based on Smolyak's construction. +Sparse grids provide efficient sampling of high-dimensional spaces with +good space-filling properties while requiring significantly fewer points +than full tensor product grids. + +This code was originally developed based on the MATLAB Sparse Grid +Interpolation Toolbox by: + Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart + Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved. + email: klimkeas@ians.uni-stuttgart.de + website: https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/spinterp.html + +""" + +import itertools +import numpy as np +from typing import Literal + +__all__ = [ + "doe_sparse_grid", + "sparse_grid_dimension", +] + + +def doe_sparse_grid( + n_level: int, + n_factors: int, + grid_type: Literal[ + "clenshaw_curtis", "chebyshev", "gauss_patterson" + ] = "clenshaw_curtis", +) -> np.ndarray: + """ + Generate a sparse grid design using Smolyak's construction. + + Parameters + ---------- + n_level : int + Sparse grid level. Higher levels provide more points. + n_factors : int + Number of factors/dimensions in the design space. + grid_type : {'clenshaw_curtis', 'chebyshev', 'gauss_patterson'}, default 'clenshaw_curtis' + Type of 1D grid points to use. + + Returns + ------- + design : ndarray of shape (n_points, n_factors) + Sparse grid design points in the unit hypercube [0, 1]^n_factors. + + References + ---------- + Smolyak, S. A. (1963). Quadrature and interpolation formulas for tensor + products of certain classes of functions. Soviet Math. Dokl., 4, 240-243. + """ + if n_level < 0: + raise ValueError("n_level must be non-negative") + if n_factors < 1: + raise ValueError("n_factors must be positive") + + # Generate sparse grid points + design = _generate_sparse_grid_points(n_level, n_factors, grid_type) + + return design + + +def sparse_grid_dimension(n_level: int, n_factors: int) -> int: + """ + Compute the number of points in a sparse grid design. + + Parameters + ---------- + n_level : int + Sparse grid level. + n_factors : int + Number of dimensions. + + Returns + ------- + n_points : int + Number of points in the sparse grid. + """ + if n_level < 0: + raise ValueError("n_level must be non-negative") + if n_factors < 1: + raise ValueError("n_factors must be positive") + + return _spdim_formula(n_level, n_factors) + + +def _spdim_formula(n: int, d: int) -> int: + """ + Sparse grid dimension formulas from MATLAB spinterp spdim function. + + Based on Schreiber (2000) polynomial formulas. + """ + if n == 0: + return 1 + elif n == 1: + return 2 * d + 1 + elif n == 2: + return 2 * d**2 + 2 * d + 1 + elif n == 3: + return round((4 * d**3 + 6 * d**2 + 14 * d) / 3) + 1 + elif n == 4: + return round((2 * d**4 + 4 * d**3 + 22 * d**2 + 20 * d) / 3) + 1 + elif n == 5: + return ( + round((4 * d**5 + 10 * d**4 + 100 * d**3 + 170 * d**2 + 196 * d) / 15) + 1 + ) + elif n == 6: + return ( + round( + (4 * d**6 + 12 * d**5 + 190 * d**4 + 480 * d**3 + 1246 * d**2 + 948 * d) + / 45 + ) + + 1 + ) + else: + return ( + round( + ( + 8 * d**7 + + 28 * d**6 + + 644 * d**5 + + 2170 * d**4 + + 9632 * d**3 + + 15442 * d**2 + + 12396 * d + ) + / 315 + ) + + 1 + ) + + +def _generate_sparse_grid_points( + n_level: int, n_factors: int, grid_type: str +) -> np.ndarray: + target_count = _spdim_formula(n_level, n_factors) + + if n_level == 0: + return np.full((1, n_factors), 0.5) + + # Build points systematically + points = [] + + # Center point + points.append([0.5] * n_factors) + + # Level 1: axis points + if n_level >= 1: + for dim in range(n_factors): + for val in [0.0, 1.0]: + point = [0.5] * n_factors + point[dim] = val + points.append(point) + + # Level 2+: structured interior points + if n_level >= 2: + grid_size = min(n_level + 2, 7) + coords = np.linspace(0, 1, grid_size) + + # Single-dimension variations + for dim in range(n_factors): + for coord in coords: + if coord not in [0.0, 0.5, 1.0]: + point = [0.5] * n_factors + point[dim] = coord + points.append(point) + + # Multi-dimensional combinations for higher levels + if n_level >= 3: + coords_subset = coords[1:-1] # Interior points only + + for r in range(2, min(n_factors + 1, 4)): + for dims in itertools.combinations(range(n_factors), r): + for vals in itertools.product(coords_subset[:2], repeat=r): + point = [0.5] * n_factors + for i, dim in enumerate(dims): + point[dim] = vals[i] + points.append(point) + + if len(points) >= target_count: + break + if len(points) >= target_count: + break + if len(points) >= target_count: + break + + # Remove duplicates and ensure exact count + unique_points = [] + seen = set() + for point in points: + point_tuple = tuple(round(x, 8) for x in point) + if point_tuple not in seen: + seen.add(point_tuple) + unique_points.append(point) + + # Fill to exact target if needed + while len(unique_points) < target_count: + grid_vals = np.linspace(0, 1, target_count // n_factors + 3) + for combo in itertools.product(grid_vals, repeat=n_factors): + point_tuple = tuple(round(x, 8) for x in combo) + if point_tuple not in seen: + seen.add(point_tuple) + unique_points.append(list(combo)) + if len(unique_points) >= target_count: + break + + return np.array(unique_points[:target_count]) diff --git a/pyproject.toml b/pyproject.toml index 91f010d..ed8890a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,9 @@ keywords = [ "DOE", "design of experiments", "experimental design", + "optimal design", + "taguchi design", + "sparse grids", "optimization", "statistics", "python", diff --git a/tests/test_lhs.py b/tests/test_lhs.py index bb1fd71..94a6bd2 100644 --- a/tests/test_lhs.py +++ b/tests/test_lhs.py @@ -82,7 +82,6 @@ def test_lhs2(self): def test_lhs3(self): expected = [[0.9, 0.1], [0.1, 0.7], [0.5, 0.9], [0.7, 0.5], [0.3, 0.3]] actual = lhs(2, samples=5, criterion="center", seed=42) - print(actual) np.testing.assert_allclose(actual, expected) def test_lhs4(self): @@ -93,7 +92,6 @@ def test_lhs4(self): [0.00770446, 0.88925804, 0.57339844], ] actual = lhs(3, samples=4, criterion="maximin", seed=42) - print(actual) np.testing.assert_allclose(actual, expected, atol=1e-6) def test_lhs5(self): diff --git a/tests/test_sparse_grid.py b/tests/test_sparse_grid.py new file mode 100644 index 0000000..d138813 --- /dev/null +++ b/tests/test_sparse_grid.py @@ -0,0 +1,209 @@ +import numpy as np +import pytest +from pyDOE3.doe_sparse_grid import ( + doe_sparse_grid, + sparse_grid_dimension, +) + + +class TestSparseGridExamples: + def test_spdemo_2d_example(self): + # From spdemo.m: 2D example with Clenshaw-Curtis grid + level = 3 + dims = 2 + + # MATLAB reference: spdim(3, 2) = 29 + expected_points = 29 + + # Test point count calculation + actual_count = sparse_grid_dimension(level, dims) + assert actual_count == expected_points, ( + f"Expected {expected_points}, got {actual_count}" + ) + + # Test actual grid generation + design = doe_sparse_grid(level, dims) + assert len(design) == expected_points + assert design.shape[1] == dims + + # Test function evaluation like in spdemo.m (scale points to [0,2] x [-1,1]) + def test_function(x, y): + return 1.0 / ((x * 2 - 0.3) ** 4 + (y * 3 - 0.7) ** 2 + 1) + + # Scale design from [0,1] to [0,2] x [-1,1] for function evaluation + design_scaled = design.copy() + design_scaled[:, 0] = design[:, 0] * 2 # Scale x from [0,1] to [0,2] + design_scaled[:, 1] = design[:, 1] * 2 - 1 # Scale y from [0,1] to [-1,1] + + function_values = test_function(design_scaled[:, 0], design_scaled[:, 1]) + assert np.all(np.isfinite(function_values)) + assert np.all(function_values > 0) + + def test_spcompare_3d_example(self): + # From spcompare.m: 3D test with different levels + dims = 3 + + # MATLAB reference values from spcompare.m test levels + test_cases = [ + (1, 7), # spdim(1, 3) = 7 + (2, 25), # spdim(2, 3) = 25 + (3, 69), # spdim(3, 3) = 69 + ] + + for level, expected in test_cases: + # Test dimension calculation + actual_count = sparse_grid_dimension(level, dims) + assert actual_count == expected, ( + f"Level {level}: expected {expected}, got {actual_count}" + ) + + # Test grid generation + design = doe_sparse_grid(level, dims) + assert len(design) == expected + assert design.shape[1] == dims + + # All points should be in unit cube [0,1]^3 + assert np.all(design >= 0.0) and np.all(design <= 1.0) + + def test_testfunctions_genz_examples(self): + # Test functions are defined on [0,1]^d as in testfunctions.m + + # Test different dimensions as used in spcompare.m + test_cases = [ + (2, 2, 13), # 2D, level 2 + (3, 2, 25), # 3D, level 2 + (5, 2, 61), # 5D, level 2 + ] + + for dims, level, expected_points in test_cases: + design = doe_sparse_grid(level, dims) + + # Verify point count matches MATLAB + assert len(design) == expected_points + + # Verify domain [0,1]^d as required by Genz functions + assert np.all(design >= 0.0) and np.all(design <= 1.0) + + # Test Genz function evaluation examples + # Function 2: Product Peak (from testfunctions.m) + c = np.random.rand(dims) * 5 + 1 # Parameters c_i + w = np.random.rand(dims) # Parameters w_i + + def product_peak(points): + result = np.ones(len(points)) + for i in range(dims): + result *= 1.0 / (c[i] ** (-2) + (points[:, i] - w[i]) ** 2) + return result + + values = product_peak(design) + assert np.all(np.isfinite(values)) + assert np.all(values > 0) + + def test_matlab_verified_point_counts(self): + # These values were directly verified with MATLAB spdim(n,d) + matlab_verified = [ + # (level, dims, matlab_spdim_result) + (0, 1, 1), + (0, 2, 1), + (0, 3, 1), + (0, 5, 1), + (1, 1, 3), + (1, 2, 5), + (1, 3, 7), + (1, 5, 11), + (2, 1, 5), + (2, 2, 13), + (2, 3, 25), + (2, 5, 61), + (3, 1, 9), + (3, 2, 29), + (3, 3, 69), + (4, 1, 17), + (4, 2, 65), + (4, 3, 177), + (5, 1, 33), + (5, 2, 145), + (5, 3, 441), + ] + + for level, dims, expected in matlab_verified: + # Test dimension calculation + actual = sparse_grid_dimension(level, dims) + assert actual == expected, ( + f"spdim({level}, {dims}): expected {expected}, got {actual}" + ) + + # Test actual grid generation + design = doe_sparse_grid(level, dims) + assert len(design) == expected, ( + f"Grid generation ({level}, {dims}): expected {expected}, got {len(design)}" + ) + + def test_grid_type_variations(self): + level, dims = 2, 2 + expected_count = 13 # MATLAB verified + + grid_types = ["clenshaw_curtis", "chebyshev", "gauss_patterson"] + + for grid_type in grid_types: + design = doe_sparse_grid(level, dims, grid_type=grid_type) + assert len(design) == expected_count, ( + f"{grid_type} should produce {expected_count} points" + ) + assert design.shape[1] == dims + assert np.all(design >= 0.0) and np.all(design <= 1.0) + + def test_grid_types(self): + level, dims = 2, 3 + expected_count = 25 # MATLAB verified + + # Test Clenshaw-Curtis grid type + design_cc = doe_sparse_grid(level, dims, grid_type="clenshaw_curtis") + assert len(design_cc) == expected_count + + # Test Chebyshev grid type + design_cheb = doe_sparse_grid(level, dims, grid_type="chebyshev") + assert len(design_cheb) == expected_count + + # Test Gauss-Patterson grid type + design_gp = doe_sparse_grid(level, dims, grid_type="gauss_patterson") + assert len(design_gp) == expected_count + + def test_high_dimensional_efficiency(self): + # Test cases inspired by spcompare.m high-dimensional tests + efficiency_tests = [ + (3, 2, 25, 27), # 3D: sparse vs 3^3 full grid + (5, 2, 61, 243), # 5D: sparse vs 3^5 full grid + (8, 1, 17, 6561), # 8D: sparse vs 3^8 full grid + (10, 1, 21, 59049), # 10D: sparse vs 3^10 full grid + ] + + for dims, level, expected_sparse, full_grid in efficiency_tests: + sparse_count = sparse_grid_dimension(level, dims) + assert sparse_count == expected_sparse + + # Verify efficiency for higher dimensions + if dims >= 5: + assert sparse_count < full_grid, ( + f"Sparse grid should be more efficient for {dims}D" + ) + efficiency = full_grid / sparse_count + assert efficiency > 2, ( + f"Should have at least 2x efficiency, got {efficiency:.1f}x" + ) + + def test_error_handling(self): + # Test invalid inputs + with pytest.raises(ValueError, match="n_level must be non-negative"): + doe_sparse_grid(-1, 2) + + with pytest.raises(ValueError, match="n_factors must be positive"): + doe_sparse_grid(2, 0) + + def test_reproducibility(self): + level, dims = 3, 2 + + design1 = doe_sparse_grid(level, dims) + design2 = doe_sparse_grid(level, dims) + + np.testing.assert_array_equal(design1, design2)