added optimal l2 projection in new ROM package

romtools/__init__.py

@@ -90,3 +90,4 @@
 from romtools.hyper_reduction import *
 from romtools.workflows import *
 from romtools.composite_vector_space import *
+from romtools.rom import *

romtools/rom/__init__.py

@@ -0,0 +1 @@
+from romtools.rom.projections import *

romtools/rom/projections.py

@@ -0,0 +1,97 @@
+#
+# ************************************************************************
+#
+#                         ROM Tools and Workflows
+# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
+#                              (NTESS)
+#
+# Under the terms of Contract DE-NA0003525 with NTESS, the
+# U.S. Government retains certain rights in this software.
+#
+# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# 3. Neither the name of the copyright holder nor the names of its
+# contributors may be used to endorse or promote products derived
+# from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+#
+# Questions? Contact Eric Parish ([email protected])
+#
+# ************************************************************************
+#
+import romtools
+import numpy as np
+from romtools.vector_space import VectorSpace
+def optimal_l2_projection(input_tensor : np.ndarray , vector_space : romtools.VectorSpace , weighting_matrix : np.ndarray = None, return_full_state=False):
+
+    '''
+    Compute L2 projection in the weighted inner product
+    arg min \| ( Phi \hat{x} + x_{ref}) - x \|_M^2 
+
+    Solution satisfies the linear system
+         Phi^T M Phi \hat{x} = \Phi^T M ( x - x_{ref} )
+    Args:
+        input_tensor (np.ndarray): 2d or 3d data array of size (n_vars, nx) or (n_vars, nx, n_snaps)
+        vector_space (romtools.VectorSpace): vector space class containing basis and affine offset
+        weighting_matrix (np.ndarray): 2d weighting matrix of size (nvars nx \times nvars nx) 
+    '''
+
+    basis = vector_space.get_basis()
+    shift_vector = vector_space.get_shift_vector()
+
+    nvars,nx,k = basis.shape
+    basis = np.reshape(basis,(nvars*nx,k))
+
+    # If input_tensor is just a single snapshot
+    if len(input_tensor.shape) == 2:
+        right_hand_side = (input_tensor - shift_vector).flatten()
+
+    # If input_tensor is many snapshots
+    elif len(input_tensor.shape) == 3:
+        right_hand_side = input_tensor - shift_vector[...,None]
+        right_hand_side = np.reshape(right_hand_side,(nvars*nx,input_tensor.shape[-1]))
+
+
+    if weighting_matrix is None:
+        left_hand_side = basis.transpose() @  basis 
+        right_hand_side = basis.transpose() @ right_hand_side
+
+    else:
+        left_hand_side = basis.transpose() @ ( weighting_matrix @ basis )
+        right_hand_side = basis.transpose() @ ( weighting_matrix  @ right_hand_side )
+
+    reduced_state = np.linalg.solve(left_hand_side,right_hand_side)
+    basis = np.reshape(basis,(nvars,nx,k))
+
+    if return_full_state == False:
+        return reduced_state 
+
+    else:
+        full_state = np.einsum('nik,k...->ni...',basis,reduced_state)
+        return reduced_state,full_state
+
+
+

tests/romtools/rom/test_projections.py

@@ -0,0 +1,199 @@
+import romtools
+import pytest
+import numpy as np
+import scipy.optimize
+
+@pytest.mark.mpi_skip
+def test_optimal_l2_projection_single_vector():
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = np.random.normal(size=(3,10))
+
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = None,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space)
+
+    phi = trial_space.get_basis()
+    def residual_for_ls_solver(xhat):
+      x = np.einsum('ijk,k...->ij...',phi,xhat)
+      error  = x.flatten() - data_to_project.flatten()
+      return error
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state)
+
+
+    ## Now test w/ an offset
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = np.random.normal(size=(3,10))
+
+    my_shifter = romtools.utils.create_average_shifter(data)
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = my_shifter,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space)
+
+    phi = trial_space.get_basis()
+    shift_vec = trial_space.get_shift_vector()
+    def residual_for_ls_solver(xhat):
+      x = np.einsum('ijk,k...->ij...',phi,xhat) + shift_vec
+      error  = x.flatten() - data_to_project.flatten()
+      return error
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state)
+
+    ## Now test w/ a weighted inner product 
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = np.random.normal(size=(3,10))
+    M = np.random.normal(size=((30,30)))
+    M = M @ M.transpose()
+    Mchol = np.linalg.cholesky(M)
+    my_shifter = romtools.utils.create_average_shifter(data)
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = my_shifter,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space, weighting_matrix=M)
+
+    phi = trial_space.get_basis()
+    shift_vec = trial_space.get_shift_vector()
+    def residual_for_ls_solver(xhat):
+      x = np.einsum('ijk,k...->ij...',phi,xhat) + shift_vec
+      error  = x.flatten() - data_to_project.flatten()
+      return Mchol.transpose() @ error
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state)
+
+
+@pytest.mark.mpi_skip
+def test_optimal_l2_projection_multiple_vectors():
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    n_data = 3
+    data_to_project = np.random.normal(size=(3,10,n_data))
+
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = None,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space)
+
+    phi = trial_space.get_basis()
+    def residual_for_ls_solver(xhat):
+      xhat = np.reshape(xhat,(rom_dim,n_data))
+      x = np.einsum('ijk,k...->ij...',phi,xhat)
+      error  = x.flatten() - data_to_project.flatten()
+      return error
+
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim*n_data),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state.flatten())
+
+
+    ## Now test w/ an offset
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = np.random.normal(size=(3,10,n_data))
+
+    my_shifter = romtools.utils.create_average_shifter(data)
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = my_shifter,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space)
+
+    phi = trial_space.get_basis()
+    shift_vec = trial_space.get_shift_vector()
+
+    def residual_for_ls_solver(xhat):
+      xhat = np.reshape(xhat,(rom_dim,n_data))
+      x = np.einsum('ijk,k...->ij...',phi,xhat) + shift_vec[...,None]
+      error  = x.flatten() - data_to_project.flatten()
+      return error
+
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim*n_data),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state.flatten())
+
+    ## Now test w/ a weighted inner product 
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = np.random.normal(size=(3,10,n_data))
+    M = np.random.normal(size=((30,30)))
+    M = M @ M.transpose()
+    Mchol = np.linalg.cholesky(M)
+    my_shifter = romtools.utils.create_average_shifter(data)
+    rom_dim = 4
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = my_shifter,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state = romtools.rom.optimal_l2_projection(data_to_project, trial_space, weighting_matrix=M)
+
+    phi = trial_space.get_basis()
+    shift_vec = trial_space.get_shift_vector()
+
+    def residual_for_ls_solver(xhat):
+      xhat = np.reshape(xhat,(rom_dim,n_data))
+      x = np.einsum('ijk,k...->ij...',phi,xhat) + shift_vec[...,None]
+      error  = x - data_to_project
+      error = np.reshape(error,((error.shape[0]*error.shape[1]),error.shape[2]))
+      return ( Mchol.transpose() @ error ).flatten()
+
+    scipy_reduced_state = scipy.optimize.least_squares(residual_for_ls_solver,np.zeros(rom_dim*n_data),jac='cs').x 
+    assert np.allclose(scipy_reduced_state,reduced_state.flatten())
+
+
+
+@pytest.mark.mpi_skip
+def test_optimal_l2_projection_full_return():
+    np.random.seed(1)
+    data = np.random.normal(size=(3,10,5))
+    data_to_project = data[...,0]
+
+    rom_dim = 5
+    my_truncater = romtools.vector_space.utils.BasisSizeTruncater(rom_dim)
+    trial_space = romtools.VectorSpaceFromPOD(snapshots=data,
+                                          truncater=my_truncater,
+                                          shifter = None,
+                                          orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
+                                          scaler = romtools.vector_space.utils.NoOpScaler())
+
+    reduced_state,projected_data = romtools.rom.optimal_l2_projection(data_to_project, trial_space,return_full_state=True)
+
+    assert np.allclose(projected_data,data_to_project)
+
+if __name__ == '__main__':
+    test_optimal_l2_projection_single_vector() 
+    test_optimal_l2_projection_multiple_vectors() 
+    test_optimal_l2_projection_full_return()