Add sparse grid designs

README.md

@@ -6,6 +6,10 @@ pyDOE3: An experimental design package for python
 [![DOI](https://zenodo.org/badge/709347557.svg)](https://zenodo.org/doi/10.5281/zenodo.10958492)
 [![Ruff](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/astral-sh/ruff/main/assets/badge/v2.json)](https://github.com/astral-sh/ruff)
 
+[![PyPI Downloads](https://img.shields.io/pypi/dm/pyDOE3.svg?label=PyPI%20downloads)](https://pypi.org/project/pyDOE3/)
+[![Conda Downloads](https://img.shields.io/conda/dn/conda-forge/pydoe3.svg?label=Conda%20downloads)](https://anaconda.org/conda-forge/pydoe3)
+[![Python versions](https://img.shields.io/pypi/pyversions/pyDOE3.svg)](https://pypi.org/project/pyDOE3/)
+
 This package is designed to help the scientist, engineer, statistician, etc., to
 construct appropriate experimental designs.
 
@@ -52,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
@@ -105,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)

doc/index.rst

@@ -52,7 +52,9 @@ number of factors:
 - :ref:`Randomized Designs <randomized>`
 
   #. :ref:`Latin-Hypercube <latin_hypercube>` (``lhs``)
+
   #. :ref:`Random K-Means <random_k_means>` (``random_k_means``)
+
   #. :ref:`Random Uniform <random_uniform>` (``random_uniform``)
 
 - :ref:`Low-Discrepancy Sequences <low_discrepancy>`
@@ -72,6 +74,7 @@ number of factors:
 - :ref:`Sampling Designs <sampling_designs>`
 
   #. :ref:`Morris Method <morris_method>` (``morris_sampling``)
+
   #. :ref:`Saltelli Sampling <saltelli_sampling>` (``saltelli_sampling``)
 
 - :ref:`Taguchi Designs <taguchi_designs>`
@@ -81,9 +84,17 @@ number of factors:
 - :ref:`Optimal Designs <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 <sparse_grids>`
+
+  #. :ref:`Sparse Grid Design <doe_sparse_grid>` (``doe_sparse_grid``)
+
+  #. :ref:`Sparse Grid Dimension <sparse_grid_dimension>` (``sparse_grid_dimension``)
+
 Requirements
 ============
 

doc/index_TOC.rst

@@ -12,3 +12,4 @@ Table of Contents
    Sampling Designs <sampling_designs>
    Taguchi Designs <taguchi>
    Optimal Designs <doe_optimal>
+   Sparse Grid Designs <sparse_grids>

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 

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").
 

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 
+<https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/help/whatis.html>`_.
+
+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 <https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/spinterp.html>`_
+- `Mathematical Details <https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/help/whatis.html>`_
+- `Implementation Guide <https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/help/getting_started.html>`_
+- `Function Reference <https://people.sc.fsu.edu/~jburkardt/m_src/spinterp/help/functions_list.html>`_

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