diff --git a/bindings/pylibROM/python_utils/cpp_utils.hpp b/bindings/pylibROM/python_utils/cpp_utils.hpp index a79a689..dffe431 100644 --- a/bindings/pylibROM/python_utils/cpp_utils.hpp +++ b/bindings/pylibROM/python_utils/cpp_utils.hpp @@ -50,4 +50,42 @@ T* getVectorPointer(py::array_t &u_in) return static_cast(buf_info.ptr); } + +template +py::buffer_info +get1DArrayBufferInfo(T *ptr, const int nelem) +{ + return py::buffer_info( + ptr, /* Pointer to buffer */ + sizeof(T), /* Size of one scalar */ + py::format_descriptor::format(), /* Python struct-style format descriptor */ + 1, /* Number of dimensions */ + { nelem }, /* Buffer dimensions */ + { sizeof(T) } /* Strides (in bytes) for each index */ + ); +} + +template +py::capsule +get1DArrayBufferHandle(T *ptr, const bool free_when_done=false) +{ + if (free_when_done) + return py::capsule(ptr, [](void *f){ + T *T_ptr = reinterpret_cast(f); + delete[] T_ptr; + }); + else + return py::capsule([](){}); +} + +template +py::array_t +get1DArrayFromPtr(T *ptr, const int nelem, bool free_when_done=false) +{ + // if empty array, no need to free when done. + free_when_done = free_when_done && (nelem > 0); + return py::array(get1DArrayBufferInfo(ptr, nelem), + get1DArrayBufferHandle(ptr, free_when_done)); +} + #endif diff --git a/bindings/pylibROM/utils/pyCSVDatabase.cpp b/bindings/pylibROM/utils/pyCSVDatabase.cpp index f475eb1..1f38b7a 100644 --- a/bindings/pylibROM/utils/pyCSVDatabase.cpp +++ b/bindings/pylibROM/utils/pyCSVDatabase.cpp @@ -99,5 +99,68 @@ void init_CSVDatabase(pybind11::module_ &m) { self.putIntegerArray(key, getVectorPointer(data), nelements); }); + csvdb.def("getIntegerArray", []( + CSVDatabase &self, const std::string& key, int nelements) + { + int *dataptr = new int[nelements]; + self.getIntegerArray(key, dataptr, nelements); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + csvdb.def("getIntegerVector", []( + CSVDatabase &self, const std::string& key, bool append) + { + std::vector *datavec = new std::vector; + self.getIntegerVector(key, *datavec, append); + return get1DArrayFromPtr(datavec->data(), datavec->size(), true); + }, + py::arg("key"), py::arg("append") = false); + + csvdb.def("getDoubleArray", []( + CSVDatabase &self, const std::string& key, int nelements) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + csvdb.def("getDoubleArray", []( + CSVDatabase &self, const std::string& key, int nelements, const std::vector& idx) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements, idx); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + csvdb.def("getDoubleArray", []( + CSVDatabase &self, const std::string& key, int nelements, + int offset, int block_size, int stride) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements, offset, block_size, stride); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + csvdb.def("getDoubleVector", []( + CSVDatabase &self, const std::string& key, bool append) + { + std::vector *datavec = new std::vector(); + self.getDoubleVector(key, *datavec, append); + return get1DArrayFromPtr(datavec->data(), datavec->size(), true); + }, + py::arg("key"), py::arg("append") = false); + + csvdb.def("getDoubleArraySize", &CSVDatabase::getDoubleArraySize); + + csvdb.def("getStringVector", []( + CSVDatabase &self, const std::string& file_name, bool append) + { + std::vector data; + self.getStringVector(file_name, data, append); + return py::cast(data); + }, py::arg("file_name"), py::arg("append") = false); + + csvdb.def("getLineCount", &CSVDatabase::getLineCount); + } diff --git a/bindings/pylibROM/utils/pyDatabase.cpp b/bindings/pylibROM/utils/pyDatabase.cpp index 6115382..b2d4d92 100644 --- a/bindings/pylibROM/utils/pyDatabase.cpp +++ b/bindings/pylibROM/utils/pyDatabase.cpp @@ -24,5 +24,11 @@ void init_Database(pybind11::module_ &m) { .export_values(); // TODO(kevin): finish binding of member functions. - + db.def("getInteger", []( + Database &self, const std::string& key) + { + int data; + self.getInteger(key, data); + return data; + }); } diff --git a/bindings/pylibROM/utils/pyHDFDatabase.cpp b/bindings/pylibROM/utils/pyHDFDatabase.cpp index 0d2f16c..769c2d4 100644 --- a/bindings/pylibROM/utils/pyHDFDatabase.cpp +++ b/bindings/pylibROM/utils/pyHDFDatabase.cpp @@ -38,6 +38,41 @@ void init_HDFDatabase(pybind11::module_ &m) { self.putIntegerArray(key, getVectorPointer(data), nelements); }); + hdfdb.def("getIntegerArray", []( + HDFDatabase &self, const std::string& key, int nelements) + { + int *dataptr = new int[nelements]; + self.getIntegerArray(key, dataptr, nelements); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + hdfdb.def("getDoubleArray", []( + HDFDatabase &self, const std::string& key, int nelements) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + hdfdb.def("getDoubleArray", []( + HDFDatabase &self, const std::string& key, int nelements, const std::vector& idx) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements, idx); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + hdfdb.def("getDoubleArray", []( + HDFDatabase &self, const std::string& key, int nelements, + int offset, int block_size, int stride) + { + double *dataptr = new double[nelements]; + self.getDoubleArray(key, dataptr, nelements, offset, block_size, stride); + return get1DArrayFromPtr(dataptr, nelements, true); + }); + + hdfdb.def("getDoubleArraySize", &HDFDatabase::getDoubleArraySize); + // hdfdb.def("__del__", [](HDFDatabase& self) { self.~HDFDatabase(); }); // Destructor } diff --git a/examples/data/inline-quad.mesh b/examples/data/inline-quad.mesh new file mode 100644 index 0000000..692520a --- /dev/null +++ b/examples/data/inline-quad.mesh @@ -0,0 +1,7 @@ +MFEM INLINE mesh v1.0 + +type = quad +nx = 4 +ny = 4 +sx = 1.0 +sy = 1.0 diff --git a/examples/dmd/heat_conduction_csv.sh b/examples/dmd/heat_conduction_csv.sh new file mode 100755 index 0000000..7d4d9f0 --- /dev/null +++ b/examples/dmd/heat_conduction_csv.sh @@ -0,0 +1,33 @@ +#!/bin/bash + +#SBATCH -N 1 +#SBATCH -t 0:05:00 +#SBATCH -p pdebug +#SBATCH -o sbatch.log +#SBATCH --open-mode truncate + +rm -rf parameters.txt dmd_list dmd_data run/dmd_data + +python3 ./parametric_heat_conduction.py -r 0.40 -save -csv -out dmd_data/dmd_par1 -no-vis > hc_par1.log +python3 ./parametric_heat_conduction.py -r 0.45 -save -csv -out dmd_data/dmd_par2 -no-vis > hc_par2.log +python3 ./parametric_heat_conduction.py -r 0.55 -save -csv -out dmd_data/dmd_par3 -no-vis > hc_par3.log +python3 ./parametric_heat_conduction.py -r 0.60 -save -csv -out dmd_data/dmd_par4 -no-vis > hc_par4.log +python3 ./parametric_heat_conduction.py -r 0.50 -save -csv -out dmd_data/dmd_par5 -no-vis > hc_par5.log + +mkdir dmd_list +mv -f run/dmd_data . + +awk 'END{print NR}' dmd_data/dmd_par1/step0/sol.csv > dmd_data/dim.csv + +mv dmd_data/dmd_par1/snap_list.csv dmd_list/dmd_par1.csv +mv dmd_data/dmd_par2/snap_list.csv dmd_list/dmd_par2.csv +mv dmd_data/dmd_par3/snap_list.csv dmd_list/dmd_par3.csv +mv dmd_data/dmd_par4/snap_list.csv dmd_list/dmd_par4.csv +mv dmd_data/dmd_par5/snap_list.csv dmd_list/dmd_par5.csv + +echo "dmd_par1,0.40,0.01,0,0" > dmd_list/dmd_train_parametric.csv +echo "dmd_par2,0.45,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "dmd_par3,0.55,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "dmd_par4,0.60,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "dmd_par5,0.50,0.01,0,0" > dmd_list/dmd_train_local.csv +echo "dmd_par5,0.50,0.01,0,0" > dmd_list/dmd_test.csv diff --git a/examples/dmd/heat_conduction_hdf.sh b/examples/dmd/heat_conduction_hdf.sh new file mode 100755 index 0000000..5583791 --- /dev/null +++ b/examples/dmd/heat_conduction_hdf.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +#SBATCH -N 1 +#SBATCH -t 0:05:00 +#SBATCH -p pdebug +#SBATCH -o sbatch.log +#SBATCH --open-mode truncate + +rm -rf parameters.txt dmd_list dmd_data run/dmd_data + +python3 ./parametric_heat_conduction.py -r 0.40 -save -hdf -out dmd_data/sim0 -no-vis > hc_par0.log +python3 ./parametric_heat_conduction.py -r 0.45 -save -hdf -out dmd_data/sim1 -no-vis > hc_par1.log +python3 ./parametric_heat_conduction.py -r 0.55 -save -hdf -out dmd_data/sim2 -no-vis > hc_par2.log +python3 ./parametric_heat_conduction.py -r 0.60 -save -hdf -out dmd_data/sim3 -no-vis > hc_par3.log +python3 ./parametric_heat_conduction.py -r 0.50 -save -hdf -out dmd_data/sim4 -no-vis > hc_par4.log + +mkdir dmd_list +mv -f run/dmd_data . + +echo "0,0.40,0.01,0,0" > dmd_list/dmd_train_parametric.csv +echo "1,0.45,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "2,0.55,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "3,0.60,0.01,0,0" >> dmd_list/dmd_train_parametric.csv +echo "4,0.50,0.01,0,0" > dmd_list/dmd_train_local.csv +echo "4,0.50,0.01,0,0" > dmd_list/dmd_test.csv diff --git a/examples/dmd/local_tw_csv.py b/examples/dmd/local_tw_csv.py new file mode 100644 index 0000000..21911e6 --- /dev/null +++ b/examples/dmd/local_tw_csv.py @@ -0,0 +1,622 @@ +# /****************************************************************************** +# * +# * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC +# * and other libROM project developers. See the top-level COPYRIGHT +# * file for details. +# * +# * SPDX-License-Identifier: (Apache-2.0 OR MIT) +# * +# *****************************************************************************/ + +# // Compile with: make local_tw_csv +# // +# // Generate CSV or HDF database on heat conduction with either +# // heat_conduction_csv.sh or heat_conduction_hdf.sh (HDF is more efficient). +# // +# // ================================================================================= +# // +# // Local serial DMD command for CSV or HDF: +# // mpirun -np 8 local_tw_csv -o hc_local_serial -rdim 16 -dtc 0.01 -csv +# // mpirun -np 8 local_tw_csv -o hc_local_serial -rdim 16 -dtc 0.01 -hdf +# // +# // Final-time prediction error (last line in run/hc_local_serial/dmd_par5_prediction_error.csv): +# // 0.0004063242226265 +# // +# // Local time windowing DMD command for CSV or HDF: +# // mpirun -np 8 local_tw_csv -o hc_local_tw -rdim 16 -nwinsamp 25 -dtc 0.01 -csv +# // mpirun -np 8 local_tw_csv -o hc_local_tw -nwinsamp 25 -dtc 0.01 -hdf +# // +# // Final-time prediction error (last line in run/hc_local_tw/dmd_par5_prediction_error.csv): +# // 0.0002458808673544 +# // +# // ================================================================================= +# // +# // Description: Local time windowing DMD on general CSV datasets. +# // +# // User specify file locations and names by -list LIST_DIR -train-set TRAIN_LIST -test-set TEST_LIST -data DATA_DIR -var VAR_NAME -o OUT_DIR +# // +# // File structure: +# // 1. LIST_DIR/TRAIN_LIST.csv -- each row specifies one training DATASET +# // 2. LIST_DIR/TEST_LIST.csv -- each row specifies one testing DATASET +# // 3. LIST_DIR/DATASET.csv -- each row specifies the suffix of one STATE in DATASET +# // 4. DATA_DIR/dim.csv -- specifies the dimension of VAR_NAME +# // 5. DATA_DIR/DATASET/tval.csv -- specifies the time instances +# // 6. DATA_DIR/DATASET/STATE/VAR_NAME.csv -- each row specifies one value of VAR_NAME of STATE +# // 7. DATA_DIR/DATASET/TEMPORAL_IDX.csv -- (optional) specifies the first and last temporal index in DATASET +# // 8. DATA_DIR/SPATIAL_IDX.csv -- (optional) each row specifies one spatial index of VAR_NAME +# // 9. run/OUT_DIR/indicator_val.csv -- (optional) each row specifies one indicator endpoint value + +import os +import io +import pathlib +import sys +try: + import mfem.par as mfem +except ModuleNotFoundError: + msg = "PyMFEM is not installed yet. Install PyMFEM:\n" + msg += "\tgit clone https://github.com/mfem/PyMFEM.git\n" + msg += "\tcd PyMFEM\n" + msg += "\tpython3 setup.py install --with-parallel\n" + raise ModuleNotFoundError(msg) + +from os.path import expanduser, join, dirname +import numpy as np +from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum + +sys.path.append("../../build") +import pylibROM.algo as algo +import pylibROM.linalg as linalg +import pylibROM.utils as utils +from pylibROM.python_utils import StopWatch + +if __name__ == "__main__": + infty = sys.maxsize + precision = 16 + + # 1. Initialize MPI. + from mpi4py import MPI + myid = MPI.COMM_WORLD.Get_rank() + num_procs = MPI.COMM_WORLD.Get_size() + + # 2. Parse command-line options. + from mfem.common.arg_parser import ArgParser + # bool train = true; + # double t_final = -1.0; + # double dtc = 0.0; + # double ddt = 0.0; + # int numWindows = 0; + # int windowNumSamples = infty; + # int windowOverlapSamples = 0; + # const char *rbf = "G"; + # const char *interp_method = "LS"; + # double admd_closest_rbf_val = 0.9; + # double ef = 0.9999; + # int rdim = -1; + # const char *list_dir = "dmd_list"; + # const char *data_dir = "dmd_data"; + # const char *sim_name = "sim"; + # const char *var_name = "sol"; + # const char *train_list = "dmd_train_local"; + # const char *test_list = "dmd_test"; + # const char *temporal_idx_list = "temporal_idx"; + # const char *spatial_idx_list = "spatial_idx"; + # const char *hdf_name = "dmd.hdf"; + # const char *snap_pfx = "step"; + # const char *basename = ""; + # bool save_csv = false; + # bool csvFormat = true; + + parser = ArgParser(description="local_tw_csv") + parser.add_argument('-vis', '--visualization', + action='store_true', default=True, + help='Enable GLVis visualization') + parser.add_argument("-train", "--train", + action='store_true', dest='train', default=True, + help="Enable DMD training.") + parser.add_argument("-no-train", "--no-train", + action='store_false', dest='train', + help="Disable DMD training.") + parser.add_argument("-tf", "--t-final", + action='store', default=-1.0, type=float, + help="Final time.") + parser.add_argument("-dtc", "--dtc", + action='store', default=0.0, type=float, + help="Fixed (constant) dt.") + parser.add_argument("-ddt", "--dtime-step", + action='store', default=0.0, type=float, + help="Desired Time step.") + parser.add_argument("-nwin", "--numwindows", + action='store', default=0, type=int, + help="Number of DMD windows.") + parser.add_argument("-nwinsamp", "--numwindowsamples", + action='store', default=infty, type=int, + help="Number of samples in DMD windows.") + parser.add_argument("-nwinover", "--numwindowoverlap", + action='store', default=0, type=int, + help="Number of samples for DMD window overlap.") + parser.add_argument("-rbf", "--radial-basis-function", + action='store', default="G", type=str, + help="Radial basis function used in interpolation. Options: \"G\", \"IQ\", \"IMQ\".") + parser.add_argument("-interp", "--interpolation-method", + action='store', default="LS", type=str, + help="Method of interpolation. Options: \"LS\", \"IDW\", \"LP\".") + parser.add_argument("-acrv", "--admd-crv", + action='store', default=0.9, type=float, + help="Adaptive DMD closest RBF value.") + parser.add_argument("-ef", "--energy-fraction", + action='store', default=0.9999, type=float, + help="Energy fraction for DMD.") + parser.add_argument("-rdim", "--rdim", + action='store', default=-1, type=int, + help="Reduced dimension for DMD.") + parser.add_argument("-list", "--list-directory", + action='store', default="dmd_list", type=str, + help="Location of training and testing data list.") + parser.add_argument("-data", "--data-directory", + action='store', default='dmd_data', type=str, + help="Location of training and testing data.") + parser.add_argument("-hdffile", "--hdf-file", + action='store', default='dmd.hdf', type=str, + help="Name of HDF file for training and testing data.") + parser.add_argument("-sim", "--sim-name", + action='store', default='sim', type=str, + help="Name of simulation.") + parser.add_argument("-var", "--variable-name", + action='store', default='sol', type=str, + help="Name of variable.") + parser.add_argument("-train-set", "--training-set-name", + action='store', default='dmd_train_local', type=str, + help="Name of the training datasets within the list directory.") + parser.add_argument("-test-set", "--testing-set-name", + action='store', default='dmd_test', type=str, + help="Name of the testing datasets within the list directory.") + parser.add_argument("-t-idx", "--temporal-index", + action='store', default='temporal_idx', type=str, + help="Name of the file indicating bound of temporal indices.") + parser.add_argument("-x-idx", "--spatial-index", + action='store', default='spatial_idx', type=str, + help="Name of the file indicating spatial indices.") + parser.add_argument("-snap-pfx", "--snapshot-prefix", + action='store', default='step', type=str, + help="Prefix of snapshots.") + parser.add_argument("-o", "--outputfile-name", + action='store', default='', type=str, + help="Name of the sub-folder to dump files within the run directory.") + parser.add_argument("-save", "--save", + action='store_true', dest='save', + help="Enable prediction result output (files in CSV format).") + parser.add_argument("-no-save", "--no-save", + action='store_false', dest='save', + help="Disable prediction result output (files in CSV format).") + parser.add_argument("-csv", "--csv", + action='store_true', dest='csvFormat', + help="Use CSV or HDF format for input files.") + parser.add_argument("-hdf", "--hdf", + action='store_false', dest='csvFormat', + help="Use CSV or HDF format for input files.") + + args = parser.parse_args() + if (myid == 0): + parser.print_options(args) + + dtc = args.dtc + ddt = args.dtime_step + t_final = args.t_final + basename = args.outputfile_name + csvFormat = args.csvFormat + var_name = args.variable_name + data_dir = args.data_directory + sim_name = args.sim_name + hdf_name = args.hdf_file + spatial_idx_list = args.spatial_index + temporal_idx_list = args.temporal_index + train = args.train + numWindows = args.numwindows + list_dir = args.list_directory + train_list = args.training_set_name + test_list = args.testing_set_name + windowNumSamples = args.numwindowsamples + windowOverlapSamples = args.numwindowoverlap + rbf = args.radial_basis_function + interp_method = args.interpolation_method + admd_closest_rbf_val = args.admd_crv + snap_pfx = args.snapshot_prefix + rdim = args.rdim + ef = args.energy_fraction + + + assert(not ((dtc > 0.0) and (ddt > 0.0))) + + if (t_final > 0.0): + save_csv = True + + outputPath = "run" + if (basename != ""): + outputPath += "/" + basename + + if (myid == 0): + pathlib.Path(outputPath).mkdir(parents=True, exist_ok=True) + + csv_db = utils.CSVDatabase() + prefix = "" + suffix = "" + if (csvFormat): + db = utils.CSVDatabase() + suffix = ".csv" + else: + db = utils.HDFDatabase() + + variable = var_name + nelements = 0 + + if (csvFormat): + nelements = db.getIntegerArray(data_dir + "/dim.csv", 1) + else: + db.open("%s/%s0/%s" % (data_dir, sim_name, hdf_name), "r") + nelements = db.getDoubleArraySize("step0sol") + db.close() + + assert(nelements > 0) + if (myid == 0): + print("Variable %s has dimension %d." % (var_name, nelements)) + + dim = nelements + idx_state_size = csv_db.getInteger("idx_state_size") + idx_state = csv_db.getIntegerArray("%s/%s.csv" % (data_dir, spatial_idx_list), idx_state_size) + if (idx_state.size > 0): + dim = idx_state.size + if (myid == 0): + print("Restricting on %d entries out of %d." % (dim, nelements)) + + if ((not train) or (numWindows > 0)): + indicator_val_size = db.getInteger("indicator_val_size") + indicator_val = db.getDoubleArray("%s/indicator_val.csv" % outputPath, + indicator_val_size) + if (indicator_val.size > 0): + if (numWindows > 0): + assert(numWindows == indicator_val.size) + else: + numWindows = indicator_val.size + if (myid == 0): + print("Read indicator range partition with %d windows." % numWindows) + + npar = 0 + # int num_train_snap_orig, num_train_snap; + # string training_par_dir; + # vector training_par_list; + # vector training_snap_bound; + if (train): + training_par_list = csv_db.getStringVector("%s/%s.csv" % (list_dir, train_list), False) + npar = training_par_list.size + assert(npar == 1) + + training_par_dir = training_par_list[0].split(',')[0] + + if (csvFormat): + num_train_snap_orig = csv_db.getLineCount("%s/%s.csv" % (list_dir, training_par_dir)) + else: + num_train_snap_orig = db.getInteger("numsnap") + assert(num_train_snap_orig > 0) + + snap_bound_size = 0 + if (csvFormat): + prefix = "%s/%s/" % (data_dir, training_par_dir) + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open("%s/%s/%s" % (data_dir, training_par_dir, hdf_name), "r") + snap_bound_size = db.getInteger("snap_bound_size") + + if (snap_bound_size > 0): + assert(snap_bound_size == 2) + training_snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + training_snap_bound[0] -= 1 + training_snap_bound[1] -= 1 + num_train_snap = training_snap_bound[1] - training_snap_bound[0] + 1 + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (training_snap_bound[0], training_snap_bound[1])) + else: + training_snap_bound = np.array([0, num_train_snap_orig - 1]) + num_train_snap = num_train_snap_orig + + db.close() + + assert(windowOverlapSamples < windowNumSamples) + numWindows = int(round(float(num_train_snap - 1) / float(windowNumSamples))) if (windowNumSamples < infty) else 1 + + assert(numWindows > 0) + if (myid == 0): + if (numWindows > 1): + print("Using time windowing DMD with %d windows." % numWindows) + else: + print("Using serial DMD.") + + dmd = [None] * numWindows + for window in range(numWindows): + if (train): + if (ddt > 0.0): + dmd[window] = algo.AdaptiveDMD(dim, ddt, rbf, + interp_method, admd_closest_rbf_val) + elif (dtc > 0.0): + dmd[window] = algo.DMD(dim, dtc) + else: + dmd[window] = algo.NonuniformDMD(dim) + else: + if (myid == 0): + print("Loading DMD model #%d." % window) + dmd[window] = algo.DMD("%s/window%d" % (outputPath, window)) + + dmd_training_timer, dmd_preprocess_timer, dmd_prediction_timer = StopWatch(), StopWatch(), StopWatch() + # double* sample = new double[dim]; + + if (train): + dmd_training_timer.Start() + + par_dir = training_par_list[0].split(',')[0] + + if (csvFormat): + snap_list = csv_db.getStringVector("%s/%s.csv" % (list_dir, par_dir), False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_train_snap_orig) + + if (myid == 0): + print("Loading samples for %s to train DMD." % par_dir) + + if (csvFormat): prefix = data_dir + "/" + par_dir + "/" + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_train_snap_orig) + + curr_window = 0 + overlap_count = 0 + for idx_snap in range(training_snap_bound[0], training_snap_bound[1] + 1): + snap = snap_pfx + tval = tvec[idx_snap] + + if (idx_snap == training_snap_bound[0]): + indicator_val = np.append(indicator_val, tval) + + if (csvFormat): + snap += snap_list[idx_snap]; # STATE + data_filename = "%s/%s/%s/%s.csv" % (data_dir, par_dir, snap, variable) # path to VAR_NAME.csv + sample = db.getDoubleArray(data_filename, nelements, idx_state) + else: + snap += "%d" % (snap_index_list[idx_snap]) # STATE + sample = db.getDoubleArray(snap + "sol", nelements, idx_state) + + dmd[curr_window].takeSample(sample, tval) + if (overlap_count > 0): + dmd[curr_window-1].takeSample(sample, tval) + overlap_count -= 1 + if ((curr_window+1 < numWindows) and (idx_snap+1 <= training_snap_bound[1])): + new_window = False + if (windowNumSamples < infty): + new_window = (idx_snap >= training_snap_bound[0] + (curr_window+1) * windowNumSamples) + else: + new_window = (tval >= indicator_val[curr_window+1]) + if (new_window): + overlap_count = windowOverlapSamples + curr_window += 1 + if (windowNumSamples < infty): + indicator_val = np.append(indicator_val, tval) + dmd[curr_window].takeSample(sample, tval) + + if (myid == 0): + print("Loaded %d samples for %s." % (num_train_snap, par_dir)) + if (windowNumSamples < infty): + print("Created new indicator range partition with %d windows." % numWindows) + csv_db.putDoubleVector("%s/indicator_val.csv" % outputPath, + indicator_val, numWindows) + + for window in range(numWindows): + if (rdim != -1): + if (myid == 0): + print("Creating DMD model #%d with rdim: %d" % (window, rdim)) + dmd[window].train(rdim) + elif (ef != -1): + if (myid == 0): + print("Creating DMD model #%d with energy fraction: %f" % (window, ef)) + dmd[window].train(ef) + dmd[window].save("%s/window%d" % (outputPath, window)) + if (myid == 0): + dmd[window].summary("%s/window%d" % (outputPath, window)) + + dmd_training_timer.Stop() + + if (ddt > 0.0): + admd = None + interp_snap = linalg.Vector(dim, True) + interp_error = [] + + for window in range(numWindows): + admd = dmd[window] + assert(admd is not None) + t_init = dmd[window].getTimeOffset() + + dtc = admd.getTrueDt() + f_snapshots = admd.getInterpolatedSnapshots() + if (myid == 0): + print("Verifying Adaptive DMD model #%d against interpolated snapshots." % window) + + for k in range(f_snapshots.numColumns()): + interp_snap = f_snapshots.getColumn(k) + result = admd.predict(t_init + k * dtc) + + dmd_solution = mfem.Vector(result.getData(), dim) + true_solution = mfem.Vector(interp_snap.getData(), dim) + diff = mfem.Vector(true_solution.Size()) + mfem.subtract_vector(dmd_solution, true_solution, diff) + + tot_diff_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff, diff)) + tot_true_solution_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, true_solution, + true_solution)) + rel_error = tot_diff_norm / tot_true_solution_norm + interp_error += [rel_error] + + if (myid == 0): + print("Norm of interpolated snapshot #%d is %f" % (k, tot_true_solution_norm)) + print("Absolute error of DMD prediction for interpolated snapshot #%d is %f" % (k, tot_diff_norm)) + print("Relative error of DMD prediction for interpolated snapshot #%d is %f" % (k, rel_error)) + del result + + if (myid == 0): + csv_db.putDoubleVector("%s/window%d_interp_error.csv" % (outputPath, window), + interp_error, f_snapshots.numColumns()) + interp_error.clear() + + db.close() + + testing_par_list = csv_db.getStringVector("%s/%s.csv" % (list_dir, test_list), False) + npar = testing_par_list.size + assert(npar > 0) + + num_tests = 0 + prediction_time, prediction_error = [], [] + + for idx_dataset in range(npar): + par_dir = testing_par_list[idx_dataset].split(',')[0] # testing DATASET + + if (myid == 0): + print("Predicting solution for %s using DMD." % par_dir) + + num_snap_orig = 0 + if (csvFormat): + num_snap_orig = csv_db.getLineCount("%s/%s.csv" % (list_dir, par_dir)) + else: + db.open(data_dir + "/" + par_dir + "/" + hdf_name, "r") + num_snap_orig = db.getInteger("numsnap") + + assert(num_snap_orig > 0) + if (csvFormat): + snap_list = csv_db.getStringVector("%s/%s.csv" % (list_dir, par_dir), False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_snap_orig) + + if (csvFormat): prefix = "%s/%s/" % (data_dir, par_dir) + + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_snap_orig) + + snap_bound_size = 0 + if (csvFormat): + prefix = "%s/%s/" % (data_dir, par_dir) + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open(data_dir + "/" + par_dir + "/" + hdf_name, "r") + snap_bound_size = db.getInteger("snap_bound_size") + + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + + if (snap_bound_size > 0): + assert(snap_bound_size == 2) + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + snap_bound[0] -= 1 + snap_bound[1] -= 1 + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (snap_bound[0], snap_bound[1])) + else: + snap_bound = np.array([0, num_snap_orig - 1]) + + num_snap = snap_bound[1] - snap_bound[0] + 1 + curr_window = 0 + for idx_snap in range(snap_bound[0], snap_bound[1] + 1): + snap = snap_pfx + tval = tvec[idx_snap] + if (csvFormat): + snap += snap_list[idx_snap]; # STATE + data_filename = "%s/%s/%s/%s.csv" % (data_dir, par_dir, snap, variable) # path to VAR_NAME.csv + sample = db.getDoubleArray(data_filename, nelements, idx_state) + else: + snap += "%d" % (snap_index_list[idx_snap]) # STATE + sample = db.getDoubleArray(snap + "sol", nelements, idx_state) + + if (myid == 0): + print("State %s read." % snap) + + if (idx_snap == 0): + dmd_preprocess_timer.Start() + init_cond = None + for window in range(numWindows): + if (myid == 0): + print("Projecting initial condition at t = %f for DMD model #%d" % (indicator_val[window], window)) + if (window == 0): + init_cond = linalg.Vector(dim, True) + for i in range(dim): + init_cond[i] = sample[i] + else: + init_cond = dmd[window-1].predict(indicator_val[window]) + dmd[window].projectInitialCondition(init_cond) + del init_cond + dmd_preprocess_timer.Stop() + + if (t_final > 0.0): # Actual prediction without true solution for comparison + num_tests += 1 + while ((curr_window+1 < numWindows) and (t_final > indicator_val[curr_window+1])): + curr_window += 1 + + if (myid == 0): + print("Predicting DMD solution at t = %f using DMD model #%d" % (t_final, curr_window)) + + dmd_prediction_timer.Start() + result = dmd[curr_window].predict(t_final) + dmd_prediction_timer.Stop() + if (myid == 0): + csv_db.putDoubleArray(outputPath + "/" + par_dir + + "_final_time_prediction.csv", + result.getData(), dim) + idx_snap = snap_bound[1]+1; # escape for-loop over idx_snap + del result + else: # Verify DMD prediction results against dataset + while ((curr_window+1 < numWindows) and (tval > indicator_val[curr_window+1])): + curr_window += 1 + + if (myid == 0): + print("Predicting DMD solution #%d at t = %f using DMD model #%d" % (idx_snap, tval, curr_window)) + + dmd_prediction_timer.Start() + result = dmd[curr_window].predict(tval) + dmd_prediction_timer.Stop() + + # Calculate the relative error between the DMD final solution and the true solution. + dmd_solution = mfem.Vector(result.getData(), result.dim()) + true_solution = mfem.Vector(sample, dim) + diff = mfem.Vector(true_solution.Size()) + mfem.subtract_vector(dmd_solution, true_solution, diff) + + tot_diff_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff, diff)) + tot_true_solution_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, true_solution, + true_solution)) + rel_error = tot_diff_norm / tot_true_solution_norm + + prediction_time += [tval] + prediction_error += [rel_error] + + if (myid == 0): + print("Norm of true solution at t = %f is %f" % (tval, tot_true_solution_norm)) + print("Absolute error of DMD solution at t = %f is %f" % (tval, tot_diff_norm)) + print("Relative error of DMD solution at t = %f is %f" % (tval, rel_error)) + if (save_csv): + csv_db.putDoubleArray("%s/%s_%s_prediction.csv" % (outputPath, par_dir, snap), + result.getData(), dim) + if (dim < nelements): + csv_db.putDoubleArray("%s/%s_%s_state.csv" % (outputPath, par_dir, snap), + sample, dim) + del result + + if ((myid == 0) and (t_final <= 0.0)): + csv_db.putDoubleVector(outputPath + "/" + par_dir + "_prediction_time.csv", + prediction_time, num_snap) + csv_db.putDoubleVector(outputPath + "/" + par_dir + "_prediction_error.csv", + prediction_error, num_snap) + + # prediction_time.clear() + # prediction_error.clear() + num_tests = num_tests + 1 if (t_final > 0.0) else num_tests + num_snap + + db.close() + + assert(num_tests > 0) + + if (myid == 0): + print("Elapsed time for training DMD: %e second\n" % dmd_training_timer.duration) + print("Elapsed time for preprocessing DMD: %e second\n" % dmd_preprocess_timer.duration) + print("Total elapsed time for predicting DMD: %e second\n" % dmd_prediction_timer.duration) + print("Average elapsed time for predicting DMD: %e second\n" % dmd_prediction_timer.duration / num_tests) + + del sample, dmd diff --git a/examples/dmd/parametric_heat_conduction.py b/examples/dmd/parametric_heat_conduction.py index 6ad9027..9e2a6e3 100644 --- a/examples/dmd/parametric_heat_conduction.py +++ b/examples/dmd/parametric_heat_conduction.py @@ -209,6 +209,9 @@ def EvalValue(self, x): parser.add_argument('-vis', '--visualization', action='store_true', default=True, help='Enable GLVis visualization') + parser.add_argument('-no-vis', '--no-visualization', + action='store_false', dest='visualization', + help='Enable GLVis visualization') parser.add_argument('-visit', '--visit-datafiles', action='store_true', default=False, help="Save data files for VisIt (visit.llnl.gov) visualization.") @@ -234,7 +237,10 @@ def EvalValue(self, x): action='store_true', default=False, help="Enable or disable MFEM DOF solution snapshot files).") parser.add_argument("-csv", "--csv", - action='store_true', default=False, + action='store_true', default=False, dest='csvFormat', + help="Use CSV or HDF format for files output by -save option.") + parser.add_argument("-hdf", "--hdf", + action='store_false', dest='csvFormat', help="Use CSV or HDF format for files output by -save option.") parser.add_argument("-out", "--outputfile-name", action='store', default="", type=str, @@ -278,7 +284,7 @@ def EvalValue(self, x): visit = args.visit_datafiles adios2 = args.adios2_streams visualization = args.visualization - csvFormat = args.csv + csvFormat = args.csvFormat dt = args.time_step t_final = args.t_final vis_steps = args.visualization_steps @@ -536,7 +542,7 @@ def EvalValue(self, x): pathlib.Path("%s/step%d" % (outputPath, ti)).mkdir(parents=True, exist_ok=True) db.putDoubleArray("%s/step%d/sol.csv" % (outputPath, ti), u.GetDataArray(), u.Size()) else: - db.putDoubleArray("step%dsol" % ti, u.GetData(), u.Size()) + db.putDoubleArray("step%dsol" % ti, u.GetDataArray(), u.Size()) ts += [t] snap_list += [ti] @@ -728,4 +734,4 @@ def EvalValue(self, x): del pwsnap_CAROM #endif - MPI.Finalize() \ No newline at end of file + MPI.Finalize() diff --git a/examples/dmd/parametric_tw_csv.py b/examples/dmd/parametric_tw_csv.py new file mode 100644 index 0000000..361ec2a --- /dev/null +++ b/examples/dmd/parametric_tw_csv.py @@ -0,0 +1,891 @@ +# ****************************************************************************** +# * +# * Copyright (c) 2013-2023, Lawrence Livermore National Security, LLC +# * and other libROM project developers. See the top-level COPYRIGHT +# * file for details. +# * +# * SPDX-License-Identifier: (Apache-2.0 OR MIT) +# * +# ***************************************************************************** +# // Compile with: make parametric_tw_csv +# // +# // Generate CSV or HDF database on heat conduction with either +# // heat_conduction_csv.sh or heat_conduction_hdf.sh (HDF is more efficient). +# // +# // ============================================================================= +# // +# // Parametric serial DMD command (for HDF version, append -hdf): +# // python3 parametric_tw_csv.py -o hc_parametric_serial -rdim 16 -dtc 0.01 -offline +# // python3 parametric_tw_csv.py -o hc_parametric_serial -rdim 16 -dtc 0.01 -online +# // +# // Final-time prediction error (Last line in run/hc_parametric_serial/dmd_par5_prediction_error.csv): +# // 0.0012598331433506 +# // +# // Parametric time windowing DMD command (for HDF version, append -hdf): +# // python3 parametric_tw_csv.py -o hc_parametric_tw -nwinsamp 25 -dtc 0.01 -offline +# // python3 parametric_tw_csv.py -o hc_parametric_tw -nwinsamp 25 -dtc 0.01 -online +# // +# // Final-time prediction error (Last line in run/hc_parametric_tw/dmd_par5_prediction_error.csv): +# // 0.0006507358659606 +# // +# // ============================================================================= +# // +# // Description: Parametric time windowing DMD on general CSV datasets. +# // +# // User specify file locations and names by -list LIST_DIR -train-set TRAIN_LIST -test-set TEST_LIST -data DATA_DIR -var VAR_NAME -o OUT_DIR +# // +# // File structure: +# // 1. LIST_DIR/TRAIN_LIST.csv -- each row specifies one training DATASET +# // 2. LIST_DIR/TEST_LIST.csv -- each row specifies one testing DATASET +# // 3. LIST_DIR/DATASET.csv -- each row specifies the suffix of one STATE in DATASET +# // 4. DATA_DIR/dim.csv -- specifies the dimension of VAR_NAME +# // 5. DATA_DIR/DATASET/tval.csv -- specifies the time instances +# // 6. DATA_DIR/DATASET/STATE/VAR_NAME.csv -- each row specifies one value of VAR_NAME of STATE +# // 7. DATA_DIR/DATASET/TEMPORAL_IDX.csv -- (optional) specifies the first and last temporal index in DATASET +# // 8. DATA_DIR/SPATIAL_IDX.csv -- (optional) each row specifies one spatial index of VAR_NAME +# // 9. run/OUT_DIR/indicator_val.csv -- (optional) each row specifies one indicator endpoint value + +import os +import io +import pathlib +import sys +try: + import mfem.par as mfem +except ModuleNotFoundError: + msg = "PyMFEM is not installed yet. Install PyMFEM:\n" + msg += "\tgit clone https://github.com/mfem/PyMFEM.git\n" + msg += "\tcd PyMFEM\n" + msg += "\tpython3 setup.py install --with-parallel\n" + raise ModuleNotFoundError(msg) + +from os.path import expanduser, join, dirname +import numpy as np +from numpy import sin, cos, exp, sqrt, pi, abs, array, floor, log, sum + +sys.path.append("../../build") +import pylibROM.algo as algo +import pylibROM.linalg as linalg +import pylibROM.utils as utils +from pylibROM.python_utils import StopWatch + +def GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, myid, mode): + if (mode == 1): + filename = "%s/%s%s_%d.hdf" % (data_dir, sim_name, par_dir, myid) + else: + filename = "%s/%s%s/%s_%d.hdf" % (data_dir, sim_name, par_dir, hdf_name, myid) + + print(filename) + return filename + +if __name__ == "__main__": + infty = sys.maxsize + precision = 16 + + # 1. Initialize MPI. + from mpi4py import MPI + myid = MPI.COMM_WORLD.Get_rank() + num_procs = MPI.COMM_WORLD.Get_size() + + # 2. Parse command-line options. + from mfem.common.arg_parser import ArgParser + + parser = ArgParser(description="parametric_tw_csv") + parser.add_argument("-offline", "--offline", + action='store_true', dest='offline', default=False, + help="Enable or disable the offline phase.") + parser.add_argument("-no-offline", "--no-offline", + action='store_false', dest='offline', + help="Enable or disable the offline phase.") + parser.add_argument("-online", "--online", + action='store_true', dest='online', default=False, + help="Enable or disable the online phase.") + parser.add_argument("-no-online", "--no-online", + action='store_false', dest='online', + help="Enable or disable the online phase.") + parser.add_argument("-predict", "--predict", + action='store_true', dest='predict', default=False, + help="Enable or disable DMD prediction.") + parser.add_argument("-no-predict", "--no-predict", + action='store_false', dest='predict', + help="Enable or disable DMD prediction.") + parser.add_argument("-tf", "--t-final", + action='store', dest='t_final', default=-1.0, type=float, + help="Final time.") + parser.add_argument("-dtc", "--dtc", + action='store', dest='dtc', default=0.0, type=float, + help="Fixed (constant) dt.") + parser.add_argument("-ddt", "--dtime-step", + action='store', dest='ddt', default=0.0, type=float, + help="Desired Time step.") + parser.add_argument("-nwin", "--numwindows", + action='store', dest='numWindows', default=0, type=int, + help="Number of DMD windows.") + parser.add_argument("-nwinsamp", "--numwindowsamples", + action='store', dest='windowNumSamples', default=infty, type=int, + help="Number of samples in DMD windows.") + parser.add_argument("-nwinover", "--numwindowoverlap", dest='windowOverlapSamples', + action='store', default=0, type=int, + help="Number of samples for DMD window overlap.") + parser.add_argument("-os", "--offset-indicator", + dest='offset_indicator', action='store_true', default=False, + help="Enable or distable the option of offset indicator.") + parser.add_argument("-no-os", "--no-offset-indicator", + dest='offset_indicator', action='store_false', + help="Enable or distable the option of offset indicator.") + parser.add_argument("-rbf", "--radial-basis-function", dest='rbf', + action='store', default='G', type=str, + help="Radial basis function used in interpolation. Options: \"G\", \"IQ\", \"IMQ\".") + parser.add_argument("-interp", "--interpolation-method", dest='interp_method', + action='store', default='LS', type=str, + help="Method of interpolation. Options: \"LS\", \"IDW\", \"LP\".") + parser.add_argument("-acrv", "--admd-crv", dest='admd_closest_rbf_val', + action='store', default=0.9, type=float, + help="Adaptive DMD closest RBF value.") + parser.add_argument("-pcrv", "--pdmd-crv", dest='pdmd_closest_rbf_val', + action='store', default=0.9, type=float, + help="Parametric DMD closest RBF value.") + parser.add_argument("-ef", "--energy-fraction", dest='ef', + action='store', default=0.9999, type=float, + help="Energy fraction for DMD.") + parser.add_argument("-rdim", "--rdim", dest='rdim', + action='store', default=-1, type=int, + help="Reduced dimension for DMD.") + parser.add_argument("-list", "--list-directory", dest='list_dir', + action='store', default='dmd_list', type=str, + help="Location of training and testing data list.") + parser.add_argument("-data", "--data-directory", dest='data_dir', + action='store', default='dmd_data', type=str, + help="Location of training and testing data.") + parser.add_argument("-hdffile", "--hdf-file", dest='hdf_name', + action='store', default='dmd', type=str, + help="Base of name of HDF file for training and testing data.") + parser.add_argument("-sim", "--sim-name", dest='sim_name', + action='store', default='sim', type=str, + help="Name of simulation.") + parser.add_argument("-var", "--variable-name", dest='var_name', + action='store', default='sol', type=str, + help="Name of variable.") + parser.add_argument("-train-set", "--training-set-name", dest='train_list', + action='store', default='dmd_train_parametric', type=str, + help="Name of the training datasets within the list directory.") + parser.add_argument("-test-set", "--testing-set-name", dest='test_list', + action='store', default='dmd_test', type=str, + help="Name of the testing datasets within the list directory.") + parser.add_argument("-t-idx", "--temporal-index", dest='temporal_idx_list', + action='store', default='temporal_idx', type=str, + help="Name of the file indicating bound of temporal indices.") + parser.add_argument("-x-idx", "--spatial-index", dest='spatial_idx_list', + action='store', default='spatial_idx', type=str, + help="Name of the file indicating spatial indices.") + parser.add_argument("-snap-pfx", "--snapshot-prefix", dest='snap_pfx', + action='store', default='step', type=str, + help="Prefix of snapshots.") + parser.add_argument("-o", "--outputfile-name", dest='basename', + action='store', default='', type=str, + help="Name of the sub-folder to dump files within the run directory.") + parser.add_argument("-save", "--save", dest='save_csv', + action='store_true', default=False, + help="Enable or disable prediction result output (files in CSV format).") + parser.add_argument("-no-save", "--no-save", dest='save_csv', + action='store_false', + help="Enable or disable prediction result output (files in CSV format).") + parser.add_argument("-save-hdf", "--save-hdf", dest='save_hdf', + action='store_true', default=False, + help="Enable or disable prediction result output (files in HDF format).") + parser.add_argument("-no-save-hdf", "--no-save-hdf", dest='save_hdf', + action='store_false', + help="Enable or disable prediction result output (files in HDF format).") + parser.add_argument("-csv", "--csv", dest='csvFormat', + action='store_true', default=True, + help="Use CSV or HDF format for input files.") + parser.add_argument("-hdf", "--hdf", dest='csvFormat', + action='store_false', + help="Use CSV or HDF format for input files.") + parser.add_argument("-wdim", "--wdim", dest='useWindowDims', + action='store_true', default=False, + help="Use DMD dimensions for each window, input from a CSV file.") + parser.add_argument("-no-wdim", "--no-wdim", dest='useWindowDims', + action='store_false', + help="Use DMD dimensions for each window, input from a CSV file.") + parser.add_argument("-subs", "--subsample", dest='subsample', + action='store', default=0, type=int, + help="Subsampling factor for training snapshots.") + parser.add_argument("-esubs", "--eval_subsample", dest='eval_subsample', + action='store', default=0, type=int, + help="Subsampling factor for evaluation.") + parser.add_argument("-hdfmode", "--hdfmodefilename", dest='fileNameMode', + action='store', default=0, type=int, + help="HDF filename mode.") + + args = parser.parse_args() + if (myid == 0): + parser.print_options(args) + + offline = args.offline + online = args.online + predict = args.predict + t_final = args.t_final + dtc = args.dtc + ddt = args.ddt + numWindows = args.numWindows + windowNumSamples = args.windowNumSamples + windowOverlapSamples = args.windowOverlapSamples + offset_indicator = args.offset_indicator + rbf = args.rbf + interp_method = args.interp_method + admd_closest_rbf_val = args.admd_closest_rbf_val + pdmd_closest_rbf_val = args.pdmd_closest_rbf_val + ef = args.ef + rdim = args.rdim + list_dir = args.list_dir + data_dir = args.data_dir + sim_name = args.sim_name + var_name = args.var_name + train_list = args.train_list + test_list = args.test_list + temporal_idx_list = args.temporal_idx_list + spatial_idx_list = args.spatial_idx_list + hdf_name = args.hdf_name + snap_pfx = args.snap_pfx + basename = args.basename + save_csv = args.save_csv + save_hdf = args.save_hdf + csvFormat = args.csvFormat + useWindowDims = args.useWindowDims + subsample = args.subsample + eval_subsample = args.eval_subsample + fileNameMode = args.fileNameMode + + assert((not (offline and online)) and (offline or online)) + assert(not ((dtc > 0.0) and (ddt > 0.0))) + dt_est = max(ddt, dtc) + + if (t_final > 0.0): + save_csv = True + + outputPath = "run" + if (basename != ""): + outputPath += "/" + basename + + if (myid == 0): + pathlib.Path(outputPath).mkdir(parents=True, exist_ok=True) + + csv_db = utils.CSVDatabase() + prefix = "" + suffix = "" + if (csvFormat): + db = utils.CSVDatabase() + suffix = ".csv" + else: + db = utils.HDFDatabase() + + if (save_hdf): + hdf_db = utils.HDFDatabase() + hdf_db.create("prediction" + (myid) + ".hdf") + + variable = var_name + nelements = 0 + + if (csvFormat): + nelements = db.getIntegerArray(data_dir + "/dim.csv", 1)[0] + else: + db.open(GetFilenameHDF(data_dir, sim_name, "0", hdf_name, + myid, fileNameMode), "r") + nelements = db.getDoubleArraySize("step0sol") + db.close() + + assert(nelements > 0) + if (myid == 0): + print("Variable %s has dimension %d." % (var_name, nelements)) + + dim = nelements + idx_state = csv_db.getIntegerVector("%s/%s.csv" % (data_dir, spatial_idx_list), False) + if (idx_state.size > 0): + dim = idx_state.size + if (myid == 0): + print("Restricting on %d entries out of %d." % (dim, nelements)) + + indicator_val = [] + if (online or (numWindows > 0)): + indicator_val = csv_db.getDoubleVector("%s/indicator_val.csv" % outputPath, False) + if (indicator_val.size > 0): + if (numWindows > 0): + assert(numWindows == indicator_val.size) + else: + numWindows = indicator_val.size + if (myid == 0): + print("Read indicator range partition with %d windows." % numWindows) + + par_dir_list = [] + par_vectors = [] + indicator_init = [] + indicator_last = [] + + training_par_list = csv_db.getStringVector((list_dir) + "/" + train_list + ".csv", False) + npar = len(training_par_list) + assert(npar > 0) + if (myid == 0): + print("Loading %d training datasets." % npar) + + dpar = -1 + for idx_dataset in range(npar): + par_info = training_par_list[idx_dataset].split(',') # training DATASET + + dpar = len(par_info) - 1 + curr_par = linalg.Vector(dpar, False) + + if (idx_dataset == 0): + assert(dpar > 0) + if (myid == 0): + print("Dimension of parametric space = %d." % dpar) + else: + assert(dpar == len(par_info) - 1) + + par_dir = par_info[0] + par_dir_list += [par_dir] + + if (not csvFormat): + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + + snap_bound_size = 0 + if (csvFormat): + prefix = (data_dir) + "/" + par_dir + "/" + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + snap_bound_size = db.getInteger("snap_bound_size") + + num_snap_orig = 0 + if (csvFormat): + num_snap_orig = csv_db.getLineCount((list_dir) + "/" + par_dir + ".csv") + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + num_snap_orig = db.getInteger("numsnap") + + assert(num_snap_orig > 0) + if (csvFormat): + snap_list = csv_db.getStringVector((list_dir) + "/" + par_dir + ".csv", False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_snap_orig) + + if (csvFormat): prefix = (data_dir) + "/" + par_dir + "/" + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + if (snap_bound.size > 0): + snap_bound[0] -= 1 + snap_bound[1] -= 1 + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (snap_bound[0], snap_bound[1])) + else: + snap_bound = np.array([0, num_snap_orig - 1]) + + for par_order in range(dpar): + curr_par[par_order] = float(par_info[par_order+1]) + par_vectors += [curr_par] + + if (offline): + if (csvFormat): prefix = (data_dir) + "/" + par_dir + "/" + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_snap_orig) + + if (offset_indicator): + indicator_init += [0.0] + indicator_last += [tvec[snap_bound[1]] - tvec[snap_bound[0]]] + else: + indicator_init += [tvec[snap_bound[0]]] + indicator_last += [tvec[snap_bound[1]]] + + db.close() + + assert(windowOverlapSamples < windowNumSamples) + if (offline and (len(indicator_val) == 0)): + indicator_min = min(indicator_init) + indicator_max = max(indicator_last) + numWindows = round((indicator_max - indicator_min) / (dt_est * windowNumSamples)) if (windowNumSamples < infty) else 1 + for window in range(numWindows): + indicator_val += [indicator_min + dt_est * windowNumSamples * window] + + if (myid == 0): + print("Created new indicator range partition with %d windows." % numWindows) + csv_db.putDoubleVector("%s/indicator_val.csv" % outputPath, + indicator_val, numWindows) + + assert(numWindows > 0) + if (myid == 0): + if (numWindows > 1): + print("Using time windowing DMD with %d windows." % numWindows) + else: + print("Using serial DMD.") + + # vector windowDim; + if (useWindowDims): + windowDim = csv_db.getIntegerVector("run/maxDim.csv") + assert(windowDim.size == numWindows) + + dmd_training_timer, dmd_preprocess_timer, dmd_prediction_timer = StopWatch(), StopWatch(), StopWatch() + dmd = [] + # vector dmd_curr_par; + # double* sample = new double[dim]; + + if (offline): + dmd_training_timer.Start() + + maxDim = [0] * numWindows + minSamp = [-1] * numWindows + + for idx_dataset in range(npar): + par_dir = par_dir_list[idx_dataset] + if (myid == 0): + print("Loading samples for %s to train DMD." % par_dir) + + dmd_curr_par = [None] * numWindows + for window in range(numWindows): + if (ddt > 0.0): + dmd_curr_par[window] = algo.AdaptiveDMD(dim, ddt, rbf, + interp_method, admd_closest_rbf_val) + elif (dtc > 0.0): + dmd_curr_par[window] = algo.DMD(dim, dtc, True) + else: + dmd_curr_par[window] = algo.NonuniformDMD(dim) + + dmd += [dmd_curr_par] + + num_snap_orig = 0 + if (csvFormat): + num_snap_orig = csv_db.getLineCount((list_dir) + "/" + par_dir + ".csv") + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + num_snap_orig = db.getInteger("numsnap") + + assert(num_snap_orig > 0) + # vector snap_list(num_snap_orig); + # vector snap_index_list(num_snap_orig); + if (csvFormat): + snap_list = csv_db.getStringVector((list_dir) + "/" + par_dir + ".csv", False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_snap_orig) + + # vector tvec(num_snap_orig); + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_snap_orig) + + snap_bound_size = 0 + if (csvFormat): + prefix = (data_dir) + "/" + par_dir + "/" + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + snap_bound_size = db.getInteger("snap_bound_size") + + # vector snap_bound(snap_bound_size); + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + + if (snap_bound.size > 0): + snap_bound[0] -= 1 + snap_bound[1] -= 1 + assert(snap_bound.size == 2) + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (snap_bound[0], snap_bound[1])) + else: + snap_bound = np.array([0, num_snap_orig - 1]) + + curr_window = 0 + overlap_count = 0 + for idx_snap in range(snap_bound[0], snap_bound[1] + 1): + if ((subsample > 1) and (idx_snap % subsample != 0) + and (idx_snap > snap_bound[0]) and (idx_snap < snap_bound[1])): + continue + + snap = snap_pfx + tval = tvec[idx_snap] + if (csvFormat): + snap += snap_list[idx_snap] # STATE + data_filename = (data_dir) + "/" + par_dir + "/" + snap + "/" + variable + ".csv" # path to VAR_NAME.csv + sample = db.getDoubleArray(data_filename, nelements, idx_state) + else: + snap += str(snap_index_list[idx_snap]) # STATE + sample = db.getDoubleArray(snap + "sol", nelements, idx_state) + + dmd[idx_dataset][curr_window].takeSample(sample, + tval - offset_indicator * tvec[snap_bound[0]]) + + if (overlap_count > 0): + dmd[idx_dataset][curr_window-1].takeSample(sample, + tval - offset_indicator * tvec[snap_bound[0]]) + overlap_count -= 1 + + # a rough estimate to correct the precision of the indicator range partition + indicator_snap = tval - offset_indicator * tvec[snap_bound[0]] + dt_est / 100.0 + + if ((curr_window+1 < numWindows) and (idx_snap+1 <= snap_bound[1]) + and (indicator_snap > indicator_val[curr_window+1])): + overlap_count = windowOverlapSamples + curr_window += 1 + dmd[idx_dataset][curr_window].takeSample(sample, + tval - offset_indicator * tvec[snap_bound[0]]) + + if (myid == 0): + print("Loaded %d samples for %s." % (snap_bound[1] - snap_bound[0] + 1, par_dir)) + + for window in range(numWindows): + if (useWindowDims): + rdim = windowDim[window] + + if (rdim != -1): + if (myid == 0): + print("Creating DMD model #%d with rdim: %d" % (window, rdim)) + dmd[idx_dataset][window].train(rdim) + elif (ef > 0.0): + if (myid == 0): + print("Creating DMD model #%d with energy fraction: %f" % (window, ef)) + dmd[idx_dataset][window].train(ef) + + if ((window > 0) and predict): + if (myid == 0): + print("Projecting initial condition at t = %f for DMD model #%d" % (indicator_val[window] + offset_indicator * tvec[snap_bound[0]], + window)) + init_cond = dmd[idx_dataset][window-1].predict(indicator_val[window]) + dmd[idx_dataset][window].projectInitialCondition(init_cond) + del init_cond + + # Make a directory for this window, only on the first parameter. + if (idx_dataset == 0): + outWindow = outputPath + "/window" + str(window) + pathlib.Path(outWindow).mkdir(parents=True, exist_ok=True) + + dmd[idx_dataset][window].save(outputPath + "/window" + + str(window) + "/par" + + str(idx_dataset)) + + if (myid == 0): + dmd[idx_dataset][window].summary(outputPath + "/window" + + str(window) + "/par" + + str(idx_dataset)) + + print("Window %d, DMD %d dim %d" % (window, idx_dataset, + dmd[idx_dataset][window].getDimension())) + + dim_w = min(dmd[idx_dataset][window].getDimension(), + dmd[idx_dataset][window].getNumSamples()-1) + maxDim[window] = max(maxDim[window], dim_w) + + if ((minSamp[window] < 0) or + (dmd[idx_dataset][window].getNumSamples() < minSamp[window])): + minSamp[window] = dmd[idx_dataset][window].getNumSamples() + # escape for-loop over window + + db.close() + + if ((not online) and (not predict)): + dmd[idx_dataset] = [None] * numWindows + # escape for-loop over idx_dataset + dmd_training_timer.Stop() + + # Limit maxDim by minSamp-1 + for window in range(numWindows): + maxDim[window] = min(maxDim[window], minSamp[window]-1) + + # Write maxDim to a CSV file + if (myid == 0): csv_db.putIntegerArray("run/maxDim.csv", maxDim, numWindows) + # escape if-statement of offline + + curr_par = linalg.Vector(dpar, False) + + if (online): + par_dir_list = [] + + dmd_preprocess_timer.Start() + testing_par_list = csv_db.getStringVector((list_dir) + "/" + test_list + ".csv", False) + npar = len(testing_par_list) + assert(npar > 0) + + if (myid == 0): + print("Loading %d testing datasets." % npar) + + dmd_curr_par = [None] * numWindows + dmd = [dmd_curr_par] * numWindows + + num_tests = 0 + # vector prediction_time, prediction_error; + + for idx_dataset in range(npar): + par_info = testing_par_list[idx_dataset].split(',') # testing DATASET + assert(dpar == len(par_info) - 1) + + par_dir = par_info[0] + par_dir_list += [par_dir] + if (myid == 0): + print("Interpolating DMD models for dataset %s" % par_dir) + + num_snap_orig = 0 + if (csvFormat): + num_snap_orig = csv_db.getLineCount((list_dir) + "/" + par_dir + ".csv") + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + num_snap_orig = db.getInteger("numsnap") + + assert(num_snap_orig > 0) + # vector snap_list(num_snap_orig); + # vector snap_index_list(num_snap_orig); + if (csvFormat): + snap_list = csv_db.getStringVector((list_dir) + "/" + par_dir + ".csv", False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_snap_orig) + + for par_order in range(dpar): + curr_par[par_order] = float(par_info[par_order+1]) + + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_snap_orig) + + snap_bound_size = 0 + if (csvFormat): + prefix = (data_dir) + "/" + par_dir + "/" + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + snap_bound_size = db.getInteger("snap_bound_size") + + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + + if (snap_bound.size > 0): + snap_bound[0] -= 1 + snap_bound[1] -= 1 + assert(snap_bound.size == 2) + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (snap_bound[0], snap_bound[1])) + else: + snap_bound = np.array([0, num_snap_orig - 1]) + + snap = snap_pfx + if (csvFormat): + snap += snap_list[snap_bound[0]]; # STATE + data_filename = (data_dir) + "/" + par_dir + "/" + snap + "/" + variable + ".csv" # path to VAR_NAME.csv + sample = db.getDoubleArray(data_filename, nelements, idx_state) + print("sample: ", sample) + print("sample size: ", sample.size) + else: + snap += str(snap_index_list[snap_bound[0]]) # STATE + sample = db.getDoubleArray(snap + "sol", nelements, idx_state) + + for window in range(numWindows): + dmd_paths = [] + for idx_trainset in range(len(par_vectors)): + dmd_paths += [outputPath + "/window" + str(window) + "/par" + str(idx_trainset)] + + if (len(par_vectors) > 1): + if (myid == 0): + print("Interpolating DMD model #%d" % window) + + dmd[idx_dataset][window] = algo.getParametricDMD(algo.DMD, par_vectors, + dmd_paths, curr_par, (rbf), + (interp_method), pdmd_closest_rbf_val) + elif ((len(par_vectors) == 1) and (dmd_paths.size() == 1)): + if (myid == 0): + print("Loading local DMD model #%d" % window) + dmd[idx_dataset][window] = algo.DMD(dmd_paths[0]) + + if (myid == 0): + print("Projecting initial condition at t = %f for DMD model #%d" + % (indicator_val[window] + offset_indicator * tvec[snap_bound[0]], window)) + + if (window == 0): + init_cond = linalg.Vector(dim, True) + for i in range(dim): + init_cond[i] = sample[i] + else: + init_cond = dmd[idx_dataset][window-1].predict(indicator_val[window]) + + dmd[idx_dataset][window].projectInitialCondition(init_cond) + + norm_init_cond = init_cond.norm() + if (myid == 0): + print("Initial condition norm %.6e for parameter %d, window %d" + % (norm_init_cond, idx_dataset, window - 1)) + + del init_cond + + if ((window > 0) and (indicator_val[window] < t_final)): + # To save memory, delete dmd[idx_dataset][window] for windows + # not containing t_final. + if (myid == 0): + print("Deleting DMD for parameter %d, window %d" + % (idx_dataset, window - 1)) + + del dmd[idx_dataset][window-1] + dmd[idx_dataset][window-1] = None + # escape for-loop over window + + db.close() + # escape for-loop over idx_dataset + dmd_preprocess_timer.Stop() + # escape if-statement of online + + if (online or predict): + num_tests = 0 + prediction_time, prediction_error = [], [] + + for idx_dataset in range(npar): + par_dir = par_dir_list[idx_dataset] + if (myid == 0): + print("Predicting solution for %s using DMD." % par_dir) + + num_snap_orig = 0 + if (csvFormat): + num_snap_orig = csv_db.getLineCount("%s/%s.csv" % (list_dir, par_dir)) + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + num_snap_orig = db.getInteger("numsnap") + + assert(num_snap_orig > 0) + if (csvFormat): + snap_list = csv_db.getStringVector("%s/%s.csv" % (list_dir, par_dir), False) + else: + snap_index_list = db.getIntegerArray("snap_list", num_snap_orig) + + tvec = db.getDoubleArray(prefix + "tval" + suffix, num_snap_orig) + + snap_bound_size = 0 + if (csvFormat): + prefix = "%s/%s/" % (data_dir, par_dir) + snap_bound_size = csv_db.getLineCount(prefix + temporal_idx_list + suffix) + else: + db.open(GetFilenameHDF(data_dir, sim_name, par_dir, hdf_name, + myid, fileNameMode), "r") + snap_bound_size = db.getInteger("snap_bound_size") + + snap_bound = db.getIntegerArray(prefix + temporal_idx_list + suffix, snap_bound_size) + + if (snap_bound_size > 0): + snap_bound[0] -= 1 + snap_bound[1] -= 1 + assert(snap_bound.size == 2) + if (myid == 0): + print("Restricting on snapshot #%d to %d." % (snap_bound[0], snap_bound[1])) + else: + snap_bound = np.array([0, num_snap_orig - 1]) + + num_snap = snap_bound[1] - snap_bound[0] + 1 + curr_window = 0 + for idx_snap in range(snap_bound[0], snap_bound[1] + 1): + if ((eval_subsample > 1) and (idx_snap % eval_subsample != 0) + and (idx_snap > snap_bound[0]) and (idx_snap < snap_bound[1])): + continue + + snap = snap_pfx + tval = tvec[idx_snap] + if (csvFormat): + snap += snap_list[idx_snap]; # STATE + data_filename = "%s/%s/%s/%s.csv" % (data_dir, par_dir, snap, variable) # path to VAR_NAME.csv + sample = db.getDoubleArray(data_filename, nelements, idx_state) + else: + snap += "%d" % (snap_index_list[idx_snap]) # STATE + sample = db.getDoubleArray(snap + "sol", nelements, idx_state) + + if (myid == 0): + print("State %s read." % snap) + + if (t_final > 0.0): # Actual prediction without true solution for comparison + num_tests += 1 + while ((curr_window+1 < numWindows) and + (t_final - offset_indicator * tvec[snap_bound[0]] > + indicator_val[curr_window+1])): + curr_window += 1 + + if (myid == 0): + print("Predicting DMD solution at t = %f using DMD model #%d" % (t_final, curr_window)) + + dmd_prediction_timer.Start() + result = dmd[idx_dataset][curr_window].predict(t_final - offset_indicator * tvec[snap_bound[0]]) + dmd_prediction_timer.Stop() + + if (myid == 0): + csv_db.putDoubleArray(outputPath + "/" + par_dir + + "_final_time_prediction.csv", + result.getData(), dim) + + if (save_hdf): + hdf_db.putDoubleArray(snap, result.getData(), dim) + + idx_snap = snap_bound[1]+1 # escape for-loop over idx_snap + del result + else: # Verify DMD prediction results against dataset + while ((curr_window+1 < numWindows) and + (tval - offset_indicator * tvec[snap_bound[0]] > + indicator_val[curr_window+1])): + curr_window += 1 + + if (myid == 0): + print("Predicting DMD solution #%d at t = %f using DMD model #%d" % (idx_snap, tval, curr_window)) + + dmd_prediction_timer.Start() + result = dmd[idx_dataset][curr_window].predict(tval - offset_indicator * tvec[snap_bound[0]]) + dmd_prediction_timer.Stop() + + # Calculate the relative error between the DMD final solution and the true solution. + dmd_solution = mfem.Vector(result.getData(), result.dim()) + true_solution = mfem.Vector(sample, dim) + diff = mfem.Vector(true_solution.Size()) + mfem.subtract_vector(dmd_solution, true_solution, diff) + + tot_diff_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, diff, diff)) + tot_true_solution_norm = sqrt(mfem.InnerProduct(MPI.COMM_WORLD, true_solution, + true_solution)) + rel_error = tot_diff_norm / tot_true_solution_norm + + prediction_time += [tval] + prediction_error += [rel_error] + + if (myid == 0): + print("Norm of true solution at t = %f is %f" % (tval, tot_true_solution_norm)) + print("Absolute error of DMD solution at t = %f is %f" % (tval, tot_diff_norm)) + print("Relative error of DMD solution at t = %f is %f" % (tval, rel_error)) + + if (idx_snap == snap_bound[1]): # Final step + with open("log.txt", "a") as f: + np.savetxt(f, [windowNumSamples, subsample, rel_error]) + + if (save_csv): + csv_db.putDoubleArray("%s/%s_%s_prediction.csv" % (outputPath, par_dir, snap), + result.getData(), dim) + if (dim < nelements): + csv_db.putDoubleArray("%s/%s_%s_state.csv" % (outputPath, par_dir, snap), + sample, dim) + + if (save_hdf): + hdf_db.putDoubleArray(snap, result.getData(), dim) + + del result + + if ((myid == 0) and (t_final <= 0.0)): + csv_db.putDoubleVector(outputPath + "/" + par_dir + "_prediction_time.csv", + prediction_time, num_snap) + csv_db.putDoubleVector(outputPath + "/" + par_dir + "_prediction_error.csv", + prediction_error, num_snap) + + # prediction_time.clear(); + # prediction_error.clear(); + num_tests = num_tests + 1 if (t_final > 0.0) else num_tests + num_snap + + db.close() + + assert(num_tests > 0) + + if (myid == 0): + print("Elapsed time for training DMD: %e second\n" % dmd_training_timer.duration) + print("Elapsed time for preprocessing DMD: %e second\n" % dmd_preprocess_timer.duration) + print("Total elapsed time for predicting DMD: %e second\n" % dmd_prediction_timer.duration) + print("Average elapsed time for predicting DMD: %e second\n" % (dmd_prediction_timer.duration / num_tests)) + # escape case (online || predict) + + if (save_hdf): + hdf_db.close() + del hdf_db + + del sample + del curr_par + del dmd