Poly operator (#83)

* Parameterized test selection via `{posargs}` so that specific tests can be executed via `tox -- my/favorite/tests.py`

* fixed typo in documentation

* added operator_dimension for arbitrary polynomial operators and tests

* added _poly_operator to __init__.py

* renamed PolyOperator to PolynomialOperator

* added datablock to PolynomialOperator

* added apply method to PolynomialOperator

* start docstring for PolynomialOperator

* some additional tests, not passing yet

* Basis Backend Maintenance (#82)

* Basis.fit_compress() method

* update basis tests

* pragma: no cover to verify() methods

* time() -> process_time() in TimedBlock

* PODBasis handle sparse weights (w/ warning) if svdsolver='dense'

* models ok inferring multiple PolynomialOperators if they have different orders

* operator_dimension() now takes in r and m, not p, and not static

* fixed issue with test_lifting, test_05: changed dimension of PolynomialOperator.apply for p=0

* added .idea to gitignore, fixed failing test for polynomial operator

* Delete .idea directory

---------

Co-authored-by: Shane <smcquar@sandia.gov>

Author: nicolearetz

.gitignore

@@ -21,6 +21,7 @@ dist/
 .vscode/
 .pypirc
 .tox/
+.idea/
 
 # Extra Matlab files
 *.mat

src/opinf/models/mono/_nonparametric.py

@@ -57,7 +57,25 @@ def _check_operator_types_unique(ops):
         """Raise a ValueError if any two operators represent the same kind
         of operation (e.g., two constant operators).
         """
-        if len({type(op) for op in ops}) != len(ops):
+        optypes = []
+        for op in ops:
+            if isinstance(op, _operators.PolynomialOperator):
+                p = op.polynomial_order
+                if p == 0:
+                    optypes.append(_operators.ConstantOperator)
+                elif p == 1:
+                    optypes.append(_operators.LinearOperator)
+                elif p == 2:
+                    optypes.append(_operators.QuadraticOperator)
+                elif p == 3:
+                    optypes.append(_operators.CubicOperator)
+                elif p == 4:
+                    optypes.append(_operators.QuarticOperator)
+                else:
+                    optypes.append(f"PolynomialOperator_ord{p}")
+            else:
+                optypes.append(type(op))
+        if len(set(optypes)) != len(optypes):
             raise ValueError("duplicate type in list of operators to infer")
 
     def _get_operator_of_type(self, OpClass):

src/opinf/operators/__init__.py

@@ -5,3 +5,4 @@
 from ._nonparametric import *
 from ._affine import *
 from ._interpolate import *
+from ._polynomial_operator import *

src/opinf/operators/_base.py

@@ -753,7 +753,7 @@ def datablock(states: np.ndarray, inputs=None) -> np.ndarray:
            \end{array}\right]
            \in \RR^{d \times k}.
 
-        Here, ``states`` is the snapshot matrix
+        Here, ``states`` is the projected snapshot matrix
         :math:`[~\qhat_0~~\cdots~~\qhat_{k-1}~]`
         and ``inputs`` is the (optional) input matrix
         :math:`[~\u_0~~\cdots~~\u_{k-1}~]`.

src/opinf/operators/_polynomial_operator.py

@@ -0,0 +1,235 @@
+# from .. import utils
+from ._base import OpInfOperator
+
+import numpy as np
+
+import scipy.linalg as la
+
+# import scipy.sparse as sparse
+from scipy.special import comb
+
+__all__ = ["PolynomialOperator"]
+
+
+class PolynomialOperator(OpInfOperator):
+    r"""Polynomial state operator
+    :math:`\Ophat_{\ell}(\qhat,\u) = \Phat[\qhat^{\otimes p}]`
+    where :math:`\otimes p` indicates the Kronecker product of a vector with
+    itself :math:`p` times. The matrix :math:`\Phat` is
+    :math:`r \times \binom{r+p-1}{p}`.
+
+    This class is equivalent to the following.
+
+    * :math:`p = 1`: :class:`ConstantOperator`
+    * :math:`p = 2`: :class:`QuadraticOperator`
+    * :math:`p = 3`: :class:`CubicOperator`
+    * :math:`p = 4`: :class:`QuarticOperator`
+
+    Parameters
+    ----------
+    polynomial_order : int
+        Order of the Kronecker product.
+    entries : (r, (r + p - 1 choose p)) ndarray or None
+        Operator matrix :math:`\Phat`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> H = opinf.operators.QuadraticOperator()
+    >>> P = opinf.operators.PolynomialOperator(2)
+    >>> entries = np.random.random((10, 100))
+    >>> H.set_entries(entries)
+    >>> H.shape
+    (10, 55)
+    >>> P.set_entries(H.entries)
+    >>> q = np.random.random(10)
+    >>> outH = H.apply(q)
+    >>> out = P.apply(q)
+    >>> np.allclose(out, outH)
+    True
+    """
+
+    def __init__(self, polynomial_order: int, entries=None):
+        """Initialize an empty operator."""
+        if polynomial_order < 0 or (
+            not np.isclose(polynomial_order, p := int(polynomial_order))
+        ):
+            raise ValueError(
+                "expected non-negative integer polynomial order"
+                + f" polynomial_order. Got p={polynomial_order}"
+            )
+
+        self.polynomial_order = p
+
+        super().__init__(entries=entries)
+
+    def operator_dimension(self, r: int, m=None) -> int:
+        """
+        computes the number of non-redundant terms in a vector of length r
+        that is taken to the power p with the Kronecker product
+
+        Parameters
+        ----------
+        r : int
+            State dimension.
+        m : int or None
+            Input dimension -- currently not used
+        """
+        if r < 0 or (not np.isclose(r, int(r))):
+            raise ValueError(
+                f"expected non-negative integer reduced dimension r. Got r={r}"
+            )
+
+        # for constant operators the dimension does not matter
+        if (p := self.polynomial_order) == 0:
+            return 1
+
+        return comb(r, p, repetition=True, exact=True)
+
+    def datablock(self, states: np.ndarray, inputs=None) -> np.ndarray:
+        r"""Return the data matrix block corresponding to
+        this operator's polynomial order,
+        with ``states`` being the projected snapshots.
+
+        Parameters
+        ----------
+        states : (r, k) or (k,) ndarray
+            State vectors. Each column is a single state vector.
+            If one dimensional, it is assumed that :math:`r = 1`.
+        inputs : (m, k) or (k,) ndarray or None
+            Input vectors (not used).
+
+        Returns
+        -------
+        datablock : (self.operator_dimension(r), k) ndarray
+            where p is the polynomial order for this operator.
+        """
+        # if constant, we just return an array containing ones
+        # of shape 1 x <number of data points>
+        if self.polynomial_order == 0:
+            return np.ones((1, np.atleast_1d(states).shape[-1]))
+
+        # make sure data is in 2D
+        states = np.atleast_2d(states)
+
+        if states.shape[0] == 0:
+            return np.empty(shape=(0, states.shape[1]))
+
+        # compute data matrix
+        return PolynomialOperator.exp_p(states, self.polynomial_order)
+
+    @staticmethod
+    def keptIndices_p(r, p):
+        """
+        returns the non-redundant indices in a kronecker-product with
+        exponent p when the dimension of the vector is r
+        """
+        if p == 0:
+            return np.array([0])
+
+        dim_if_p_was_one_smaller = PolynomialOperator(
+            polynomial_order=(p - 1)
+        ).operator_dimension(r=r)
+        indexmatrix = np.reshape(
+            np.arange(r * dim_if_p_was_one_smaller),
+            (r, dim_if_p_was_one_smaller),
+        )
+        return np.hstack(
+            [
+                indexmatrix[
+                    i,
+                    : PolynomialOperator(
+                        polynomial_order=(p - 1)
+                    ).operator_dimension(i + 1),
+                ]
+                for i in range(r)
+            ]
+        )
+
+    @staticmethod
+    def exp_p(x, p, kept=None):
+        """
+        recursively computes x^p without the redundant terms
+        (it still computes them but then takes them out)
+        the result has shape
+
+        if x is 1-dimensional:
+        (PolynomialOperator.operator_dimension(x.shape[0]),)
+
+        otherwise:
+        (PolynomialOperator.operator_dimension(x.shape[0]), x.shape[1])
+        """
+        # for a constant operator, we just return 1 (x^0 = 1)
+        if p == 0:
+            return np.ones([1])
+
+        # for a linear operator, x^1 = 1
+        if p == 1:
+            return x
+
+        # identify kept entries in condensed Kronecker product for
+        # this reduced dimension
+        # for all polynomial orders up to self.polynomial order
+        if kept is None:
+            r = x.shape[0]
+            kept = [
+                PolynomialOperator.keptIndices_p(r=r, p=i)
+                for i in range(p + 1)
+            ]
+
+        # distinguish between the shapes of the input
+        if len(x.shape) == 1:
+            # this gets called when the ROM is run
+            return np.kron(x, PolynomialOperator.exp_p(x, p - 1, kept))[
+                kept[p]
+            ]
+        else:
+            # this gets called for constructing the data matrix
+            return la.khatri_rao(x, PolynomialOperator.exp_p(x, p - 1, kept))[
+                kept[p]
+            ]
+
+    def apply(self, state: np.ndarray, input_=None) -> np.ndarray:
+        r"""Apply the operator to the given state. Input is not used.
+        See OpInfOperator.apply for description.
+        """
+        if state.shape[0] != self.state_dimension:
+            raise ValueError(
+                f"Expected state of dimension r={self.state_dimension}."
+                + f"Got state.shape={state.shape}"
+            )
+
+        # constant
+        if self.polynomial_order == 0:
+            if np.ndim(state) == 2:  # r, k > 1.
+                return np.outer(self.entries, np.ones(state.shape[-1]))
+            if np.ndim(self.entries) == 2:
+                return self.entries[:, 0]
+            return self.entries
+        # note: no need to go through the trouble of identifying the
+        # non-redundant indices
+
+        # linear
+        if self.polynomial_order == 1:
+            return self.entries @ state
+        # note: no need to go through the trouble of identifying the
+        # non-redundant indices
+
+        # higher-order
+        restricted_kronecker_product = PolynomialOperator.exp_p(
+            x=state, p=self.polynomial_order, kept=self.nonredudant_entries
+        )
+        return self.entries @ restricted_kronecker_product
+
+    # Properties --------------------------------------------------------------
+    @property
+    def nonredudant_entries(self) -> list:
+        r"""list containing at index i a list of the indices that are kept
+        when restricting the i-times Kronecker product of a vector of
+        shape self.state_dimension() with itself.
+        """
+        # return self.__nonredudant_entries
+        return [
+            PolynomialOperator.keptIndices_p(r=self.state_dimension, p=i)
+            for i in range(self.polynomial_order + 1)
+        ]

tests/operators/test_nonparametric.py

@@ -33,6 +33,19 @@ def get_entries(self, r, m=None):
     def test_set_entries(self):
         raise NotImplementedError
 
+    def test_in_model(self, r=3, k=100, m=2):
+        """See if we can fit a model with this operator."""
+        model = opinf.models.ContinuousModel(operators=[self.get_operator()])
+
+        Q = np.random.random((r, k))
+        dQ = np.random.random((r, k))
+
+        if self.has_inputs:
+            U = np.random.random((m, k))
+            model.fit(states=Q, ddts=dQ, inputs=U)
+        else:
+            model.fit(states=Q, ddts=dQ, inputs=None)
+
 
 # No dependence on state or input =============================================
 class TestConstantOperator(_TestNonparametricOperator):