Projection-based ROM example

.github/workflows/ci.yml

@@ -101,6 +101,20 @@ jobs:
             echo run pyMatrix unit test
             pytest test_pyMatrix.py --verbose
             mpirun -n 2 pytest test_pyMatrix.py --verbose
+            echo run pyOptions unit test
+            pytest test_pyOptions.py --verbose
+            echo run pySVD unit test
+            pytest test_pySVD.py --verbose
+            echo run pyStaticSVD unit test
+            pytest test_pyStaticSVD.py --verbose
+            echo run pyIncrementalSVD unit test
+            pytest test_pyIncrementalSVD.py --verbose
+            echo run pyBasisGenerator unit test
+            pytest test_pyBasisGenerator.py --verbose
+            echo run pyBasisReader unit test
+            pytest test_pyBasisReader.py --verbose
+            echo run pyBasisWriter unit test
+            pytest test_pyBasisWriter.py --verbose
   baseline:
     runs-on: ubuntu-latest
     needs: [docker-base-image]

CMakeLists.txt

@@ -1,19 +1,34 @@
 cmake_minimum_required(VERSION 3.12)
 
-project(pylibROM)
+project(_pylibROM)
 
 set(CMAKE_BUILD_TYPE Debug)
 set(PYBIND11_FINDPYTHON ON)
 
+#=================== ScaLAPACK (optional) ==================
+option(BUILD_SCALAPACK "Build static ScaLAPACK for libROM" OFF)
+
+if (BUILD_SCALAPACK)
+	set(WORK_DIR ${CMAKE_CURRENT_SOURCE_DIR}/extern)
+	add_custom_command(
+		OUTPUT SLPK_BUILD
+		WORKING_DIRECTORY ${WORK_DIR}
+		COMMAND ${WORK_DIR}/librom_scalapack.sh
+		COMMENT "Building Static ScaLAPACK..."
+	)
+	add_custom_target(RUN_SLPK_BUILD ALL DEPENDS SLPK_BUILD)
+endif(BUILD_SCALAPACK)
+
 #=================== libROM ================================
 
 #It is tedious to build libROM. option to not build
 set(LIBROM_DIR "" CACHE STRING "absolute path to the pre-installed libROM")
-set(BUILD_DEPS OFF)
+
+set(BUILD_LIBROM OFF)
 if (NOT LIBROM_DIR)
   set(LIBROM_DIR ${CMAKE_CURRENT_SOURCE_DIR}/extern/libROM)
   set(LIBROM_SCRIPTS_DIR ${LIBROM_DIR}/scripts)
-  set(BUILD_DEPS ON)
+  set(BUILD_LIBROM ON)
 endif()
 
 include(ExternalProject)
@@ -22,7 +37,7 @@ include(CMakePrintHelpers)
 cmake_print_variables(LIBROM_DIR)
 cmake_print_variables(LIBROM_SCRIPTS_DIR)
 
-if (BUILD_DEPS)
+if (BUILD_LIBROM)
 	# add_custom_command(
 	# 	OUTPUT LIBROM_BUILD
 	# 	WORKING_DIRECTORY ${LIBROM_DIR}
@@ -41,7 +56,7 @@ if (BUILD_DEPS)
         	INSTALL_COMMAND ""
 	)	
 	message("Building libROM dependency...")
-endif(BUILD_DEPS)
+endif(BUILD_LIBROM)
 
 #setup external dependency build and link paths for libROM
 set(LIBROM_INCLUDE_DIR ${LIBROM_DIR}/lib)
@@ -50,6 +65,36 @@ link_directories(${LIBROM_DIR}/build/lib) #this hack is the best way for non-cma
 #include mpi4py directory
 execute_process(COMMAND python3 -c "import mpi4py; print(mpi4py.get_include())" OUTPUT_VARIABLE MPI4PY)
 
+# # MFEM is required.
+# # TODO(kevin): We do not bind mfem-related functions until we figure out how to type-cast SWIG Object.
+# #              Until then, mfem-related functions need to be re-implemented on python-end, using PyMFEM.
+
+# find_library(MFEM mfem
+# 	"$ENV{MFEM_DIR}/lib"
+# 	"$ENV{MFEM_DIR}"
+# 	"${LIBROM_DIR}/dependencies/mfem")
+# find_library(HYPRE HYPRE
+# 	"$ENV{HYPRE_DIR}/lib"
+# 	"${LIBROM_DIR}/dependencies/hypre/src/hypre/lib")
+# find_library(PARMETIS parmetis
+# 	"$ENV{PARMETIS_DIR}/lib"
+# 	"$ENV{PARMETIS_DIR}/build/lib/libparmetis"
+# 	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libparmetis")
+# find_library(METIS metis
+# 	"$ENV{METIS_DIR}/lib"
+# 	"$ENV{PARMETIS_DIR}/build/lib/libmetis"
+# 	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/build/lib/libmetis")
+# find_path(MFEM_INCLUDES mfem.hpp
+# 	"$ENV{MFEM_DIR}/include"
+# 	"$ENV{MFEM_DIR}"
+# 	"${LIBROM_DIR}/dependencies/mfem")
+# find_path(HYPRE_INCLUDES HYPRE.h
+# 	"$ENV{HYPRE_DIR}/include"
+# 	"${LIBROM_DIR}/dependencies/hypre/src/hypre/include")
+# find_path(PARMETIS_INCLUDES metis.h
+# 	"$ENV{PARMETIS_DIR}/metis/include"
+# 	"${LIBROM_DIR}/dependencies/parmetis-4.0.3/metis/include")
+
 #===================== pylibROM =============================
 
 
@@ -58,14 +103,26 @@ set(CMAKE_CXX_STANDARD 14)
 find_package(MPI REQUIRED)
 
 set(SOURCE_DIR "bindings/pylibROM") 
-include_directories(	${SOURCE_DIR} 
-			${LIBROM_INCLUDE_DIR} 
-			${MPI_INCLUDE_PATH}
-			${MPI4PY})
-add_subdirectory("extern/pybind11")
+include_directories(
+	${SOURCE_DIR} 
+	${LIBROM_INCLUDE_DIR} 
+	${MPI_INCLUDE_PATH}
+	${MPI4PY}
+	# ${MFEM_INCLUDES}
+	# ${HYPRE_INCLUDES}
+	# ${PARMETIS_INCLUDES}
+	# ${MFEM_C_INCLUDE_DIRS}
+)
+# link_libraries(
+# 	${MFEM}
+# 	${HYPRE}
+# 	${PARMETIS}
+# 	${METIS}
+# )
 
+add_subdirectory("extern/pybind11")
 
-pybind11_add_module(pylibROM
+pybind11_add_module(_pylibROM
     bindings/pylibROM/pylibROM.cpp
 
     bindings/pylibROM/linalg/pyMatrix.cpp
@@ -79,7 +136,14 @@ pybind11_add_module(pylibROM
     bindings/pylibROM/linalg/svd/pyIncrementalSVD.cpp
     bindings/pylibROM/algo/pyDMD.cpp 
 	bindings/pylibROM/utils/mpi_utils.cpp
+	bindings/pylibROM/utils/pyDatabase.cpp
+
+	# TODO(kevin): We do not bind mfem-related functions until we figure out how to type-cast SWIG Object.
+	#              Until then, mfem-related functions need to be re-implemented on python-end, using PyMFEM.
+	# bindings/pylibROM/mfem/Utilities.cpp
+
+	bindings/pylibROM/python_utils/cpp_utils.hpp
 )
 message("building pylibROM...")
 
-target_link_libraries(pylibROM PRIVATE ROM)
+target_link_libraries(_pylibROM PRIVATE ROM)

README.md

@@ -5,7 +5,7 @@ Python Interface for LLNL libROM
 
 1. Pull repository and all sub-module dependencies:
   ```
-  git clone --recurse-submodules https://github.com/sullan2/pylibROM.git
+  git clone --recurse-submodules https://github.com/llnl/pylibROM.git
   ```
 
 2. Compile and build pylibROM (from top-level pylibROM repo):
@@ -16,11 +16,15 @@ Python Interface for LLNL libROM
   ```
   pip install ./ --global-option="--librom_dir=/path/to/pre-installed-libROM"
   ```  
+  If you want to build static ScaLAPACK for libROM,
+  ```
+  pip install ./ --global-option="--install_scalapack"
+  ```
 
 3. Test python package (from top-level pylibROM repo):
   ```
   cd tests
-  python3.6 testVector.py
+  pytest test_pyVector.py
   ```
 
 ### Using PyMFEM

bindings/pylibROM/__init__.py

@@ -1,3 +1,5 @@
-from . import linalg
-from . import algo
-from . import utils
+# To add pure python routines to this module,
+# either define/import the python routine in this file.
+# This will combine both c++ bindings/pure python routines into this module.
+
+from _pylibROM import *

bindings/pylibROM/algo/__init__.py

@@ -1 +1,6 @@
-from .pyDMD import DMD
+# To add pure python routines to this module,
+# either define/import the python routine in this file.
+# This will combine both c++ bindings/pure python routines into this module.
+
+# For other c++ binding modules, change the module name accordingly.
+from _pylibROM.algo import *

bindings/pylibROM/algo/pyDMD.cpp

@@ -7,6 +7,7 @@
 #include <pybind11/stl.h>
 #include "algo/DMD.h"
 #include "linalg/Vector.h"
+#include "python_utils/cpp_utils.hpp"
 
 namespace py = pybind11;
 using namespace CAROM;
@@ -53,11 +54,9 @@ void init_DMD(pybind11::module_ &m) {
         }), py::arg("dim"), py::arg("dt"), py::arg("alt_output_basis") = false, py::arg("vec") = nullptr)
 
     // .def("setOffset", &PyDMD::setOffset, py::arg("offset_vector"), py::arg("order"))  //problem if we want to name the wrapper as DMD. Could get rid of the using namespace directive?
-    .def("takeSample", [](DMD &self, py::array_t<double> u_in, double t) {
-            py::buffer_info buf_info = u_in.request();
-            double* data = static_cast<double*>(buf_info.ptr);
-            self.takeSample(data, t);
-        })
+    .def("takeSample", [](DMD &self, py::array_t<double> &u_in, double t) {
+        self.takeSample(getVectorPointer(u_in), t);
+    })
     .def("train", py::overload_cast<double, const Matrix*, double>(&DMD::train),
             py::arg("energy_fraction"), py::arg("W0") = nullptr, py::arg("linearity_tol") = 0.0)
     .def("train", py::overload_cast<int, const Matrix*, double>(&DMD::train),

bindings/pylibROM/linalg/__init__.py

@@ -1,7 +1,6 @@
-from .pyVector import Vector
-from .pyMatrix import Matrix 
-from .pyBasisWriter import BasisWriter
-from .pyBasisReader import BasisReader
-from .pyBasisGenerator import BasisGenerator 
-from .pyOptions import  Options 
-from . import svd 
+# To add pure python routines to this module,
+# either define/import the python routine in this file.
+# This will combine both c++ bindings/pure python routines into this module.
+
+# For other c++ binding modules, change the module name accordingly.
+from _pylibROM.linalg import *

bindings/pylibROM/linalg/pyBasisGenerator.cpp

@@ -3,42 +3,36 @@
 #include <pybind11/operators.h>
 #include <pybind11/stl.h>
 #include "linalg/BasisGenerator.h"
-
+#include "python_utils/cpp_utils.hpp"
 
 namespace py = pybind11;
 using namespace CAROM;
 
 void init_BasisGenerator(pybind11::module_ &m) {
-    py::enum_<Database::formats>(m, "Formats")
-        .value("HDF5", Database::formats::HDF5);
 
     py::class_<BasisGenerator>(m, "BasisGenerator")
         // .def(py::init<Options, bool, const std::string&, Database::formats>())
-        .def(py::init<const Options&, bool, const std::string&, Database::formats>(),py::arg("options"),py::arg("incremental"),py::arg("basis_file_name") = "",py::arg("file_format") = static_cast<int>(Database::formats::HDF5))
+        .def(py::init<const Options&, bool, const std::string&, Database::formats>(),
+            py::arg("options"),
+            py::arg("incremental"),
+            py::arg("basis_file_name") = "",
+            py::arg("file_format") = Database::formats::HDF5
+        )
         .def("isNextSample", (bool (BasisGenerator::*)(double)) &BasisGenerator::isNextSample)
         .def("updateRightSV", (bool (BasisGenerator::*)()) &BasisGenerator::updateRightSV)
-        .def("takeSample", [](BasisGenerator& self, py::array_t<double> u_in, double time, double dt, bool add_without_increase = false) {
-            py::buffer_info buf_info = u_in.request();
-            if (buf_info.ndim != 1)
-                throw std::runtime_error("Input array must be 1-dimensional");
-
-            double* u_in_data = static_cast<double*>(buf_info.ptr);
-            return self.takeSample(u_in_data, time, dt, add_without_increase);
+        .def("takeSample", [](BasisGenerator& self, py::array_t<double> &u_in, double time, double dt, bool add_without_increase = false) {
+            return self.takeSample(getVectorPointer(u_in), time, dt, add_without_increase);
         }, py::arg("u_in"), py::arg("time"), py::arg("dt"), py::arg("add_without_increase") = false)
         .def("endSamples", &BasisGenerator::endSamples, py::arg("kind") = "basis")
         .def("writeSnapshot", (void (BasisGenerator::*)()) &BasisGenerator::writeSnapshot)
-        .def("loadSamples", &BasisGenerator::loadSamples, py::arg("base_file_name"), py::arg("kind") = "basis", py::arg("cut_off") = 1e9, py::arg("db_format") = static_cast<int>(Database::formats::HDF5))
-        .def("computeNextSampleTime", [](BasisGenerator& self, py::array_t<double> u_in, py::array_t<double> rhs_in, double time) {
-             py::buffer_info buf_info_u = u_in.request();
-             py::buffer_info buf_info_rhs = rhs_in.request();
-
-             if (buf_info_u.ndim != 1 || buf_info_rhs.ndim != 1)
-               throw std::runtime_error("Input arrays must be 1-dimensional");
-
-             double* u_in_data = static_cast<double*>(buf_info_u.ptr);
-             double* rhs_in_data = static_cast<double*>(buf_info_rhs.ptr);
-
-             return self.computeNextSampleTime(u_in_data, rhs_in_data, time);
+        .def("loadSamples", (void (BasisGenerator::*)(const std::string&, const std::string&, int, Database::formats)) &BasisGenerator::loadSamples,
+            py::arg("base_file_name"),
+            py::arg("kind") = "basis",
+            py::arg("cut_off") = static_cast<int>(1e9),
+            py::arg("db_format") = Database::formats::HDF5
+        )
+        .def("computeNextSampleTime", [](BasisGenerator& self, py::array_t<double> &u_in, py::array_t<double> &rhs_in, double time) {
+             return self.computeNextSampleTime(getVectorPointer(u_in), getVectorPointer(rhs_in), time);
         }, py::arg("u_in"), py::arg("rhs_in"), py::arg("time"))
 
         .def("getSpatialBasis", (const Matrix* (BasisGenerator::*)()) &BasisGenerator::getSpatialBasis,py::return_value_policy::reference)

bindings/pylibROM/linalg/pyBasisReader.cpp

@@ -12,7 +12,10 @@ using namespace CAROM;
 
 void init_BasisReader(pybind11::module_ &m) {
      py::class_<BasisReader>(m, "BasisReader")
-        .def(py::init<const std::string&, Database::formats>(),py::arg("base_file_name"),py::arg("file_format") = static_cast<int>(Database::formats::HDF5))
+        .def(py::init<const std::string&, Database::formats>(),
+            py::arg("base_file_name"),
+            py::arg("file_format") = Database::formats::HDF5
+        )
         .def("isNewBasis",(bool (BasisReader::*)(double)) &BasisReader::isNewBasis)
         .def("getSpatialBasis",(Matrix* (BasisReader::*)(double)) &BasisReader::getSpatialBasis)
         .def("getSpatialBasis",(Matrix* (BasisReader::*)(double,int)) &BasisReader::getSpatialBasis)
@@ -30,4 +33,4 @@ void init_BasisReader(pybind11::module_ &m) {
         .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(double,int)) &BasisReader::getSnapshotMatrix)
         .def("getSnapshotMatrix",(Matrix* (BasisReader::*)(double,int,int)) &BasisReader::getSnapshotMatrix)
         .def("__del__", [](BasisReader& self) { self.~BasisReader(); }); // Destructor
-}
+}

bindings/pylibROM/linalg/pyMatrix.cpp

@@ -32,7 +32,7 @@ void init_matrix(pybind11::module_ &m) {
         .def(py::init<>())
         .def(py::init<int, int, bool, bool>(),
             py::arg("num_rows"), py::arg("num_cols"), py::arg("distributed"), py::arg("randomized") = false)
-        .def(py::init([](py::array_t<double> mat, bool distributed, bool copy_data = true) {
+        .def(py::init([](py::array_t<double> &mat, bool distributed, bool copy_data = true) {
             py::buffer_info buf_info = mat.request();
             int num_rows = buf_info.shape[0];     
             int num_cols = buf_info.shape[1];     

bindings/pylibROM/linalg/pyVector.cpp

@@ -35,7 +35,7 @@ void init_vector(pybind11::module_ &m) {
         .def(py::init<int, bool>())
 
         // Constructor
-        .def(py::init([](py::array_t<double> vec, bool distributed, bool copy_data = true) {
+        .def(py::init([](py::array_t<double> &vec, bool distributed, bool copy_data = true) {
             py::buffer_info buf_info = vec.request();
             int dim = buf_info.shape[0];
             double* data = static_cast<double*>(buf_info.ptr);

bindings/pylibROM/linalg/svd/__init__.py

@@ -1,3 +1,6 @@
-from .pySVD import SVD
-from .pyStaticSVD import StaticSVD 
-from .pyIncrementalSVD import IncrementalSVD
+# To add pure python routines to this module,
+# either define/import the python routine in this file.
+# This will combine both c++ bindings/pure python routines into this module.
+
+# For other c++ binding modules, change the module name accordingly.
+from _pylibROM.linalg.svd import *