From 83b21a22d47a30cb1bccd61ef504a3913075d640 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 16 Apr 2024 22:19:20 -0700 Subject: [PATCH 01/17] Implemented snapshot interpolation, and initial pass at integrating with DMD. Pivoting the DMD changes to derived class SnapshotDMD in future commits. Turned off "window overlap" in parametric_dw_csv. --- CMakeLists.txt | 2 + examples/dmd/parametric_dw_csv.cpp | 11 +- lib/CMakeLists.txt | 1 + lib/algo/DMD.cpp | 36 +++ lib/algo/DMD.h | 2 + lib/algo/SnapshotDMD.cpp | 77 +++++++ lib/algo/SnapshotDMD.h | 20 ++ .../manifold_interp/SnapshotInterpolator.cpp | 209 ++++++++++++++++++ .../manifold_interp/SnapshotInterpolator.h | 51 +++++ 9 files changed, 407 insertions(+), 2 deletions(-) create mode 100644 lib/algo/SnapshotDMD.cpp create mode 100644 lib/algo/SnapshotDMD.h create mode 100644 lib/algo/manifold_interp/SnapshotInterpolator.cpp create mode 100644 lib/algo/manifold_interp/SnapshotInterpolator.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a9d9ffa3..f11498fe8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,6 +155,7 @@ if (ENABLE_EXAMPLES) local_dw_csv parametric_tw_csv parametric_dw_csv + SnapshotInterpolatorTest parametric_dmdc_heat_conduction) set(example_directories prom @@ -182,6 +183,7 @@ if (ENABLE_EXAMPLES) dmd dmd dmd + dmd dmd) list(LENGTH examples len1) diff --git a/examples/dmd/parametric_dw_csv.cpp b/examples/dmd/parametric_dw_csv.cpp index 7222ea85f..6adccbf2b 100644 --- a/examples/dmd/parametric_dw_csv.cpp +++ b/examples/dmd/parametric_dw_csv.cpp @@ -504,8 +504,9 @@ int main(int argc, char *argv[]) } int ns = dmd[idx_dataset][curr_window]->getNumSamples(); - overlap_count.push_back(windowOverlapSamples + - max(0, rdim+1 - ns)); + //overlap_count.push_back(windowOverlapSamples + + // max(0, rdim+1 - ns)); + overlap_count.push_back(0); curr_window += 1; dmd[idx_dataset][curr_window]->takeSample(sample, tval); @@ -550,6 +551,12 @@ int main(int argc, char *argv[]) { if (rdim != -1) { + if(!myid) + { + cout << "DMD model #" << window << "currently has " << dmd[idx_dataset][window]->getNumSamples() << " samples" << std::endl; + cout << "Interpolating DMD model #" << window << " to contain " << rdim+1 << " snapshots" << endl; + } + dmd[idx_dataset][window]->interpolateToNSnapshots(rdim+1); if (myid == 0) { cout << "Creating DMD model #" << window << " with rdim: " << rdim << endl; diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 0e546061d..ac3e66a8e 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -45,6 +45,7 @@ set(module_list algo/manifold_interp/Interpolator algo/manifold_interp/MatrixInterpolator algo/manifold_interp/VectorInterpolator + algo/manifold_interp/SnapshotInterpolator hyperreduction/DEIM hyperreduction/GNAT hyperreduction/QDEIM diff --git a/lib/algo/DMD.cpp b/lib/algo/DMD.cpp index a66f3b65b..a91d2a594 100755 --- a/lib/algo/DMD.cpp +++ b/lib/algo/DMD.cpp @@ -15,6 +15,7 @@ #include "linalg/Matrix.h" #include "linalg/Vector.h" #include "linalg/scalapack_wrapper.h" +#include "manifold_interp/SnapshotInterpolator.h" #include "utils/Utilities.h" #include "utils/CSVDatabase.h" #include "utils/HDFDatabase.h" @@ -211,6 +212,41 @@ void DMD::train(double energy_fraction, const Matrix* W0, double linearity_tol) delete f_snapshots; } +void DMD::interpolateToNSnapshots(int n) +{ + SnapshotInterpolator* interp = new SnapshotInterpolator(); + std::vector new_snapshots; + std::vector new_times; + CSVDatabase* csv_db(new CSVDatabase); + int dim = d_snapshots[0]->dim(); + +/* + std::string debugPath = "/usr/workspace/ptranq/dmdStuff"; + if(!d_rank) + { + std::cout << "Interpolating DMD snapshots with dimension " << dim << std::endl; + for(int i = 0; i < d_snapshots.size(); ++i) + { + csv_db->putDoubleArray(debugPath+"/snapshot_" + std::to_string(i) + "pre_interp.csv",d_snapshots[i]->getData(),dim); + } + } +*/ + new_snapshots = interp->interpolate(d_sampled_times,d_snapshots,n,&new_times); + d_snapshots = new_snapshots; + d_sampled_times = new_times; + d_dt = d_sampled_times[2]->getData()[0]-d_sampled_times[1]->getData()[0]; +/* + if(!d_rank) + { + for(int i = 0; i < d_snapshots.size(); ++i) + { + csv_db->putDoubleArray(debugPath+"/snapshot_" + std::to_string(i) + "post_interp.csv",d_snapshots[i]->getData(),dim); + } + } +*/ + +} + void DMD::train(int k, const Matrix* W0, double linearity_tol) { const Matrix* f_snapshots = getSnapshotMatrix(); diff --git a/lib/algo/DMD.h b/lib/algo/DMD.h index ba3750251..6511fc1b8 100644 --- a/lib/algo/DMD.h +++ b/lib/algo/DMD.h @@ -178,6 +178,8 @@ class DMD */ const Matrix* getSnapshotMatrix(); + void interpolateToNSnapshots(int n); + /** * @brief Load the object state from a file. * diff --git a/lib/algo/SnapshotDMD.cpp b/lib/algo/SnapshotDMD.cpp new file mode 100644 index 000000000..77788a6ee --- /dev/null +++ b/lib/algo/SnapshotDMD.cpp @@ -0,0 +1,77 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, 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) + * + *****************************************************************************/ + +// Description: Implementation of the AdaptiveDMD algorithm. + +#include "AdaptiveDMD.h" +#include "manifold_interp/VectorInterpolator.h" +#include "manifold_interp/SnapshotInterpolator.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include + +namespace CAROM { + +AdaptiveDMD::~AdaptiveDMD() +{ + for (auto interp_snapshot : d_interp_snapshots) + { + delete interp_snapshot; + } +} + +void AdaptiveDMD::train(double energy_fraction, const Matrix* W0, + double linearity_tol) +{ + const Matrix* f_snapshots = getInterpolatedSnapshots(); + CAROM_VERIFY(f_snapshots->numColumns() > 1); + CAROM_VERIFY(energy_fraction > 0 && energy_fraction <= 1); + d_energy_fraction = energy_fraction; + constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol); + + delete f_snapshots; +} + +void AdaptiveDMD::train(int k, const Matrix* W0, double linearity_tol) +{ + const Matrix* f_snapshots = getInterpolatedSnapshots(); + CAROM_VERIFY(f_snapshots->numColumns() > 1); + CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1); + d_energy_fraction = -1.0; + d_k = k; + constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol); + + delete f_snapshots; +} + +void AdaptiveDMD::interpolateSnapshots() +{ + CAROM_VERIFY(d_interp_snapshots.size() == 0); + CAROM_VERIFY(d_snapshots.size() == d_sampled_times.size()); + CAROM_VERIFY(d_sampled_times.size() > 1); + + SnapshotInterpolator* interp = new SnapshotInterpolator(); + std::vector new_times; + + d_interp_snapshots = interp->interpolate(d_sampled_times,d_snapshots,n,&new_times); +} + +double AdaptiveDMD::getTrueDt() const +{ + return d_dt; +} + +const Matrix* AdaptiveDMD::getInterpolatedSnapshots() +{ + if (d_interp_snapshots.size() == 0) interpolateSnapshots(); + return createSnapshotMatrix(d_interp_snapshots); +} + +} diff --git a/lib/algo/SnapshotDMD.h b/lib/algo/SnapshotDMD.h new file mode 100644 index 000000000..ce4527db3 --- /dev/null +++ b/lib/algo/SnapshotDMD.h @@ -0,0 +1,20 @@ +#ifndef included_SnapshotDMD_h +#define included_SnapshotDMD_h + +#include "DMD.h" +#include + +namespace CAROM { + +class Vector; +class SnapshotDMD : public DMD +{ +public: + +private: + std::vector d_interp_snapshots; + + void interpolateSnapshots(); +}; +} +#endif \ No newline at end of file diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp new file mode 100644 index 000000000..b511e7d18 --- /dev/null +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -0,0 +1,209 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, 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) + * + *****************************************************************************/ + +// Description: Implementation of the VectorInterpolator algorithm. + +#include "SnapshotInterpolator.h" + +#include +#include +#include "linalg/Vector.h" +#include "mpi.h" + +/* Use C++11 built-in shared pointers if available; else fallback to Boost. */ +#if __cplusplus >= 201103L +#include +#else +#include +#endif +using namespace std; + +namespace CAROM { + +SnapshotInterpolator::SnapshotInterpolator() +{ + +} + + +std::vector SnapshotInterpolator::interpolate(std::vector snapshot_ts, + std::vector snapshots, + std::vector output_ts) +{ + CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); + CAROM_VERIFY(snapshot_ts.size() > 0); + CAROM_VERIFY(output_ts.size() > 0); + + int n_out = output_ts.size(); + int n_snap = snapshots.size(); + int n_dim = snapshots[0]->dim(); + + std::vector output_snapshots; + for(int i = 0; i < n_out; ++i) + { + Vector* temp_snapshot =new Vector(snapshots[0]->dim(), + snapshots[0]->distributed()); + output_snapshots.push_back(temp_snapshot); + } + + for(int i = 0; i < n_dim; ++i) + { + double h_temp,delta_temp, t; + std::vector h,d,delta,t_in,y_in; + for(int j = 0; j < n_snap-1; ++j) + { + h_temp = snapshot_ts[j+1]->getData()[0] - snapshot_ts[j]->getData()[0]; + h.push_back(h_temp); + delta_temp = (snapshots[j+1]->getData()[i] - snapshots[j]->getData()[i])/h_temp; + delta.push_back(delta_temp); + t_in.push_back(snapshot_ts[j]->getData()[0]); + y_in.push_back(snapshots[j]->getData()[i]); + } + t_in.push_back(snapshot_ts[n_snap-1]->getData()[0]); + y_in.push_back(snapshots[n_snap-1]->getData()[i]); + + double d_temp = ((2*h[0] + h[1])*delta[0] - h[0]*delta[1])/(h[0]+h[1]); + if(sign(d_temp)!=sign(delta[0])) + { + d_temp = 0; + } + else if (sign(delta[0]) != sign(delta[1]) && abs(d_temp) > abs(3*delta[0])) + { + d_temp = 3*delta[0]; + } + d.push_back(d_temp); + int counter = 0; + for(int j = 0; j < n_snap-2; ++j) + { + d_temp = computeDerivative(delta[j],delta[j+1],h[j],h[j+1]); + d.push_back(d_temp); + while(output_ts[counter]->getData()[0] <= t_in[j+1]) + { + t = output_ts[counter]->getData()[0]; + output_snapshots[counter]->getData()[i] = y_in[j]*computeH1(t,t_in[j],t_in[j+1]) + + y_in[j+1]*computeH2(t,t_in[j],t_in[j+1]) + + d[j]*computeH3(t,t_in[j],t_in[j+1]) + + d[j+1]*computeH4(t,t_in[j],t_in[j+1]); + counter += 1; + } + } + + d_temp = ((2*h[n_snap-2]+h[n_snap-3])*delta[n_snap-2] - h[n_snap-2]*delta[n_snap-3])/(h[n_snap-2]+h[n_snap-3]); + if(sign(d_temp) != sign(delta[n_snap-2])) + { + d_temp = 0; + } + else if (sign(delta[n_snap-2]) != sign(delta[n_snap-3]) && abs(d_temp) > abs(3*delta[n_snap-2])) + { + d_temp = 3*delta[n_snap-2]; + } + d.push_back(d_temp); + while(counter < n_out && output_ts[counter]->getData()[0] <= t_in[n_snap-1]) + { + t = output_ts[counter]->getData()[0]; + output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t,t_in[n_snap-2],t_in[n_snap-1]) + + y_in[n_snap-1]*computeH2(t,t_in[n_snap-2],t_in[n_snap-1]) + + d[n_snap-2]*computeH3(t,t_in[n_snap-2],t_in[n_snap-1]) + + d[n_snap-1]*computeH4(t,t_in[n_snap-2],t_in[n_snap-1]); + counter += 1; + } + } + return output_snapshots; +} + +std::vector SnapshotInterpolator::interpolate(std::vector snapshot_ts, + std::vector snapshots, + int n_out, + std::vector* output_ts) +{ + CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); + CAROM_VERIFY(snapshot_ts.size() > 0); + CAROM_VERIFY(n_out > 0); + + int n_snap = snapshots.size(); + int n_dim = snapshots[0]->dim(); + + + + double t_min = snapshot_ts[0]->getData()[0]; + double t_max = snapshot_ts[n_snap-1]->getData()[0]; + double dt = (t_max-t_min)/(n_out-1); + output_ts->clear(); + + + for(int i = 0; i < n_out; ++i) + { + Vector* temp_t = new Vector(1, false); + temp_t->getData()[0] = t_min + i*dt; + output_ts->push_back(temp_t); + } + std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap << std::endl; + return interpolate(snapshot_ts,snapshots,*output_ts); +} + +double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1, double h2) +{ + double d = 0.0; + double alpha = (h1 + 2*h2)/(3*(h1+h2)); + if(S1*S2 > 0) + { + d = S1*S2/(alpha*S2 + (1-alpha)*S1); + } + else + { + d = 0.0; + } + return d; +} + +double SnapshotInterpolator::computeH1(double x, double xl, double xr) +{ + double h = xr - xl; + return computePhi((xr-x)/h); +} + +double SnapshotInterpolator::computeH2(double x, double xl, double xr) +{ + double h = xr - xl; + return computePhi((x-xl)/h); +} + +double SnapshotInterpolator::computeH3(double x, double xl, double xr) +{ + double h = xr-xl; + return -h*computePsi((xr-x)/h); +} + +double SnapshotInterpolator::computeH4(double x, double xl, double xr) +{ + double h = xr-xl; + return h*computePsi((x-xl)/h); +} + +double SnapshotInterpolator::computePhi(double t) +{ + return 3.*pow(t,2.) - 2*pow(t,3.); +} + +double SnapshotInterpolator::computePsi(double t) +{ + return pow(t,3.) - pow(t,2.); +} + +int SnapshotInterpolator::sign(double a) +{ + double TOL = 1e-15; + if(abs(a) < TOL)return 0; + else if(a > 0) return 1; + else if (a < 0) return -1; + return 0; +} + +} \ No newline at end of file diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.h b/lib/algo/manifold_interp/SnapshotInterpolator.h new file mode 100644 index 000000000..943afb973 --- /dev/null +++ b/lib/algo/manifold_interp/SnapshotInterpolator.h @@ -0,0 +1,51 @@ +#ifndef included_SnapshotInterpolator_h +#define included_SnapshotInterpolator_h + +#include "Interpolator.h" +#include +#include + +namespace CAROM { + +class Matrix; +class Vector; + +/** + * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) + * + */ +class SnapshotInterpolator +{ +public: + SnapshotInterpolator(); + ~SnapshotInterpolator(); + + std::vector interpolate(std::vector snapshot_ts, + std::vector snapshots, + std::vector output_ts); + + std::vector interpolate(std::vector snapshot_ts, + std::vector snapshots, + int n_out, + std::vector* output_ts); + + std::vector interpolate(); + + + +private: + + double computeDerivative(double S1, double S2, double h1, double h2); + double computeH1(double x, double xl, double xr); + double computeH2(double x, double xl, double xr); + double computeH3(double x, double xl, double xr); + double computeH4(double x, double xl, double xr); + double computePhi(double t); + double computePsi(double t); + int sign(double a); +}; + +} +#endif \ No newline at end of file From b91ee833b4d5cb6a24cc2b787204d488d35e33d0 Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 12:58:48 -0700 Subject: [PATCH 02/17] Reduced scope to only the interpolator and a unit test. Will construct SnapshotDMD as a separat PR. --- CMakeLists.txt | 5 +- examples/dmd/parametric_dw_csv.cpp | 11 +- lib/algo/DMD.cpp | 36 ----- lib/algo/DMD.h | 2 - lib/algo/SnapshotDMD.cpp | 77 ----------- lib/algo/SnapshotDMD.h | 20 --- unit_tests/test_SnapshotInterpolator.cpp | 165 +++++++++++++++++++++++ 7 files changed, 169 insertions(+), 147 deletions(-) delete mode 100644 lib/algo/SnapshotDMD.cpp delete mode 100644 lib/algo/SnapshotDMD.h create mode 100644 unit_tests/test_SnapshotInterpolator.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f11498fe8..b5424af3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,7 +155,6 @@ if (ENABLE_EXAMPLES) local_dw_csv parametric_tw_csv parametric_dw_csv - SnapshotInterpolatorTest parametric_dmdc_heat_conduction) set(example_directories prom @@ -183,7 +182,6 @@ if (ENABLE_EXAMPLES) dmd dmd dmd - dmd dmd) list(LENGTH examples len1) @@ -283,7 +281,8 @@ if(GTEST_FOUND) IncrementalSVDBrand GreedyCustomSampler NNLS - basis_conversion) + basis_conversion + SnapshotInterpolator) foreach(stem IN LISTS unit_test_stems) add_executable(test_${stem} unit_tests/test_${stem}.cpp) target_link_libraries(test_${stem} PRIVATE ROM diff --git a/examples/dmd/parametric_dw_csv.cpp b/examples/dmd/parametric_dw_csv.cpp index 6adccbf2b..7222ea85f 100644 --- a/examples/dmd/parametric_dw_csv.cpp +++ b/examples/dmd/parametric_dw_csv.cpp @@ -504,9 +504,8 @@ int main(int argc, char *argv[]) } int ns = dmd[idx_dataset][curr_window]->getNumSamples(); - //overlap_count.push_back(windowOverlapSamples + - // max(0, rdim+1 - ns)); - overlap_count.push_back(0); + overlap_count.push_back(windowOverlapSamples + + max(0, rdim+1 - ns)); curr_window += 1; dmd[idx_dataset][curr_window]->takeSample(sample, tval); @@ -551,12 +550,6 @@ int main(int argc, char *argv[]) { if (rdim != -1) { - if(!myid) - { - cout << "DMD model #" << window << "currently has " << dmd[idx_dataset][window]->getNumSamples() << " samples" << std::endl; - cout << "Interpolating DMD model #" << window << " to contain " << rdim+1 << " snapshots" << endl; - } - dmd[idx_dataset][window]->interpolateToNSnapshots(rdim+1); if (myid == 0) { cout << "Creating DMD model #" << window << " with rdim: " << rdim << endl; diff --git a/lib/algo/DMD.cpp b/lib/algo/DMD.cpp index a91d2a594..a66f3b65b 100755 --- a/lib/algo/DMD.cpp +++ b/lib/algo/DMD.cpp @@ -15,7 +15,6 @@ #include "linalg/Matrix.h" #include "linalg/Vector.h" #include "linalg/scalapack_wrapper.h" -#include "manifold_interp/SnapshotInterpolator.h" #include "utils/Utilities.h" #include "utils/CSVDatabase.h" #include "utils/HDFDatabase.h" @@ -212,41 +211,6 @@ void DMD::train(double energy_fraction, const Matrix* W0, double linearity_tol) delete f_snapshots; } -void DMD::interpolateToNSnapshots(int n) -{ - SnapshotInterpolator* interp = new SnapshotInterpolator(); - std::vector new_snapshots; - std::vector new_times; - CSVDatabase* csv_db(new CSVDatabase); - int dim = d_snapshots[0]->dim(); - -/* - std::string debugPath = "/usr/workspace/ptranq/dmdStuff"; - if(!d_rank) - { - std::cout << "Interpolating DMD snapshots with dimension " << dim << std::endl; - for(int i = 0; i < d_snapshots.size(); ++i) - { - csv_db->putDoubleArray(debugPath+"/snapshot_" + std::to_string(i) + "pre_interp.csv",d_snapshots[i]->getData(),dim); - } - } -*/ - new_snapshots = interp->interpolate(d_sampled_times,d_snapshots,n,&new_times); - d_snapshots = new_snapshots; - d_sampled_times = new_times; - d_dt = d_sampled_times[2]->getData()[0]-d_sampled_times[1]->getData()[0]; -/* - if(!d_rank) - { - for(int i = 0; i < d_snapshots.size(); ++i) - { - csv_db->putDoubleArray(debugPath+"/snapshot_" + std::to_string(i) + "post_interp.csv",d_snapshots[i]->getData(),dim); - } - } -*/ - -} - void DMD::train(int k, const Matrix* W0, double linearity_tol) { const Matrix* f_snapshots = getSnapshotMatrix(); diff --git a/lib/algo/DMD.h b/lib/algo/DMD.h index 6511fc1b8..ba3750251 100644 --- a/lib/algo/DMD.h +++ b/lib/algo/DMD.h @@ -178,8 +178,6 @@ class DMD */ const Matrix* getSnapshotMatrix(); - void interpolateToNSnapshots(int n); - /** * @brief Load the object state from a file. * diff --git a/lib/algo/SnapshotDMD.cpp b/lib/algo/SnapshotDMD.cpp deleted file mode 100644 index 77788a6ee..000000000 --- a/lib/algo/SnapshotDMD.cpp +++ /dev/null @@ -1,77 +0,0 @@ -/****************************************************************************** - * - * Copyright (c) 2013-2024, 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) - * - *****************************************************************************/ - -// Description: Implementation of the AdaptiveDMD algorithm. - -#include "AdaptiveDMD.h" -#include "manifold_interp/VectorInterpolator.h" -#include "manifold_interp/SnapshotInterpolator.h" -#include "linalg/Matrix.h" -#include "linalg/Vector.h" -#include - -namespace CAROM { - -AdaptiveDMD::~AdaptiveDMD() -{ - for (auto interp_snapshot : d_interp_snapshots) - { - delete interp_snapshot; - } -} - -void AdaptiveDMD::train(double energy_fraction, const Matrix* W0, - double linearity_tol) -{ - const Matrix* f_snapshots = getInterpolatedSnapshots(); - CAROM_VERIFY(f_snapshots->numColumns() > 1); - CAROM_VERIFY(energy_fraction > 0 && energy_fraction <= 1); - d_energy_fraction = energy_fraction; - constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol); - - delete f_snapshots; -} - -void AdaptiveDMD::train(int k, const Matrix* W0, double linearity_tol) -{ - const Matrix* f_snapshots = getInterpolatedSnapshots(); - CAROM_VERIFY(f_snapshots->numColumns() > 1); - CAROM_VERIFY(k > 0 && k <= f_snapshots->numColumns() - 1); - d_energy_fraction = -1.0; - d_k = k; - constructDMD(f_snapshots, d_rank, d_num_procs, W0, linearity_tol); - - delete f_snapshots; -} - -void AdaptiveDMD::interpolateSnapshots() -{ - CAROM_VERIFY(d_interp_snapshots.size() == 0); - CAROM_VERIFY(d_snapshots.size() == d_sampled_times.size()); - CAROM_VERIFY(d_sampled_times.size() > 1); - - SnapshotInterpolator* interp = new SnapshotInterpolator(); - std::vector new_times; - - d_interp_snapshots = interp->interpolate(d_sampled_times,d_snapshots,n,&new_times); -} - -double AdaptiveDMD::getTrueDt() const -{ - return d_dt; -} - -const Matrix* AdaptiveDMD::getInterpolatedSnapshots() -{ - if (d_interp_snapshots.size() == 0) interpolateSnapshots(); - return createSnapshotMatrix(d_interp_snapshots); -} - -} diff --git a/lib/algo/SnapshotDMD.h b/lib/algo/SnapshotDMD.h deleted file mode 100644 index ce4527db3..000000000 --- a/lib/algo/SnapshotDMD.h +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef included_SnapshotDMD_h -#define included_SnapshotDMD_h - -#include "DMD.h" -#include - -namespace CAROM { - -class Vector; -class SnapshotDMD : public DMD -{ -public: - -private: - std::vector d_interp_snapshots; - - void interpolateSnapshots(); -}; -} -#endif \ No newline at end of file diff --git a/unit_tests/test_SnapshotInterpolator.cpp b/unit_tests/test_SnapshotInterpolator.cpp new file mode 100644 index 000000000..b20b75379 --- /dev/null +++ b/unit_tests/test_SnapshotInterpolator.cpp @@ -0,0 +1,165 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, 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): +// parametric_tw_csv -o hc_parametric_serial -rdim 16 -dtc 0.01 -offline +// parametric_tw_csv -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): +// parametric_tw_csv -o hc_parametric_tw -nwinsamp 25 -dtc 0.01 -offline +// parametric_tw_csv -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 + +#include "mfem.hpp" +#include "algo/DMD.h" +#include "algo/AdaptiveDMD.h" +#include "algo/NonuniformDMD.h" +#include "algo/manifold_interp/SnapshotInterpolator.h" +#include "linalg/Vector.h" +#include "linalg/Matrix.h" +#include "utils/HDFDatabase.h" +#include "utils/CSVDatabase.h" +#include +#include +#include + +#ifndef _WIN32 +#include // mkdir +#else +#include // _mkdir +#define mkdir(dir, mode) _mkdir(dir) +#endif + +using namespace std; +using namespace CAROM; + +int main(int argc, char *argv[]) +{ + int n_out = 25; + int n_snap = 11; + //input times + vector t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15}; + //Test function from original PCHIP paper + vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; + //sin(x) + vector y2{0, 0.909297426825682, 0.141120008059867, -0.958924274663138, -0.279415498198926, + 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, + 0.990607355694870, 0.650287840157117}; + + //query points + vector tq{0, 0.6000, 1.2000, 1.8000, 2.4000, 3.0000, 3.6000, 4.2000, 4.8000, 5.4000, 6.0000, 6.6000, + 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, + 13.2000, 13.8000, 14.4000, 15.000}; + //pchip reference solution for original test function. + vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, + 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, + 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, + 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, + 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, + 50.000000000000000, 53.832363636363638, 55.720727272727267, 58.433818181818161, + 67.055999999999969, 85.000000000000000}; + //pchip reference solution for sin(x) + vector yq2{0, 0.569748887844739, 0.833039030477175, 0.906694689302138, 0.701592409322499, 0.141120008059867, + -0.288488198334352, -0.697095554949408, -0.939878053603591, -0.782971083698006, -0.279415498198926, + 0.188293444380395, 0.669217685146469, 0.965688937709033, 0.846474732609843, 0.412118485241757, + -0.077580693288345, -0.623537711025344, -0.971758328554163, -0.890773577229575, -0.536572918000435, + -0.041614069121016, 0.560852411851254, 0.957953731078007, 0.938810668593539, 0.650287840157117}; + + std::vector snapshots; + std::vector out_snapshots; + std::vector reference_snapshots; + std::vector times; + std::vector out_times; + for(int i = 0; i < t.size(); ++i) + { + Vector* temp_t = new Vector(1, false); + temp_t->item(0) = t[i]; + times.push_back(temp_t); + Vector* temp_y = new Vector(2,false); + temp_y->item(0) = y[i]; + temp_y->item(1) = y2[i]; + snapshots.push_back(temp_y); + } + + for(int i = 0; i < tq.size(); ++i) + { + Vector* temp_t = new Vector(1, false); + temp_t->item(0) = tq[i]; + out_times.push_back(temp_t); + Vector* temp_y = new Vector(2,false); + temp_y->item(0) = yq1[i]; + temp_y->item(1) = yq2[i]; + reference_snapshots.push_back(temp_y); + } + + SnapshotInterpolator* interp = new SnapshotInterpolator(); + + std::cout << "Beginning base interpolation function" << std::endl; + out_snapshots = interp->interpolate(times,snapshots,out_times); + std::cout << "Finished interpolation " << std::endl; +/* + for(int i = 0; i < out_snapshots.size(); ++i) + { + std::cout << "Time " << i << " is " << out_times[i]->item(0) << std::endl; + std::cout << "Reference at " << i << " is (" << reference_snapshots[i]->item(0) << + "," << reference_snapshots[i]->item(1) << ")" << std::endl; + std::cout << "Snapshot interpolation at " << i << " is (" << out_snapshots[i]->item(0) << + "," << out_snapshots[i]->item(1) << ")" << std::endl; + } + */ + for(int i = 0; i < out_snapshots.size(); ++i) + { + std::cout << "Error at time[" << i << "] = " << out_times[i]->item(0) << " is (" << + reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," + << reference_snapshots[i]->item(1) - out_snapshots[i]->item(1) << ") " << std::endl; + } + + + std::cout << "Beginning variant interpolation function" << std::endl; + out_snapshots = interp->interpolate(times,snapshots,26,&out_times); + std::cout << "Finished interpolation " << std::endl; + + for(int i = 0; i < out_snapshots.size(); ++i) + { + std::cout << "Error at time[" << i << "] = " << out_times[i]->item(0) << " is (" << + reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," + << reference_snapshots[i]->item(1) - out_snapshots[i]->item(1) << ") " << std::endl; + } + +} \ No newline at end of file From 36a487cad74133dba79918205881d24a3a8fd6cb Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 13:40:41 -0700 Subject: [PATCH 03/17] Fixed unit test compilation. --- CMakeLists.txt | 6 +-- ...tor.cpp => snapshot_interpolator_test.cpp} | 47 ------------------- 2 files changed, 3 insertions(+), 50 deletions(-) rename unit_tests/{test_SnapshotInterpolator.cpp => snapshot_interpolator_test.cpp} (69%) diff --git a/CMakeLists.txt b/CMakeLists.txt index b5424af3f..f34a8e490 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,7 +228,8 @@ if (ENABLE_EXAMPLES) weak_scaling random_test smoke_static - load_samples) + load_samples + snapshot_interpolator_test) if (USE_MFEM) set(regression_test_names @@ -281,8 +282,7 @@ if(GTEST_FOUND) IncrementalSVDBrand GreedyCustomSampler NNLS - basis_conversion - SnapshotInterpolator) + basis_conversion) foreach(stem IN LISTS unit_test_stems) add_executable(test_${stem} unit_tests/test_${stem}.cpp) target_link_libraries(test_${stem} PRIVATE ROM diff --git a/unit_tests/test_SnapshotInterpolator.cpp b/unit_tests/snapshot_interpolator_test.cpp similarity index 69% rename from unit_tests/test_SnapshotInterpolator.cpp rename to unit_tests/snapshot_interpolator_test.cpp index b20b75379..fcaa424fc 100644 --- a/unit_tests/test_SnapshotInterpolator.cpp +++ b/unit_tests/snapshot_interpolator_test.cpp @@ -8,56 +8,10 @@ * *****************************************************************************/ -// 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): -// parametric_tw_csv -o hc_parametric_serial -rdim 16 -dtc 0.01 -offline -// parametric_tw_csv -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): -// parametric_tw_csv -o hc_parametric_tw -nwinsamp 25 -dtc 0.01 -offline -// parametric_tw_csv -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 - -#include "mfem.hpp" -#include "algo/DMD.h" -#include "algo/AdaptiveDMD.h" -#include "algo/NonuniformDMD.h" #include "algo/manifold_interp/SnapshotInterpolator.h" #include "linalg/Vector.h" -#include "linalg/Matrix.h" -#include "utils/HDFDatabase.h" -#include "utils/CSVDatabase.h" #include #include -#include #ifndef _WIN32 #include // mkdir @@ -161,5 +115,4 @@ int main(int argc, char *argv[]) reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," << reference_snapshots[i]->item(1) - out_snapshots[i]->item(1) << ") " << std::endl; } - } \ No newline at end of file From 6f4c57a23ed06c2b118f0664f41059e45546ef85 Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 13:42:55 -0700 Subject: [PATCH 04/17] Stylize --- .../manifold_interp/SnapshotInterpolator.cpp | 52 +++++++++++-------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index b511e7d18..c27e3c175 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -33,9 +33,10 @@ SnapshotInterpolator::SnapshotInterpolator() } -std::vector SnapshotInterpolator::interpolate(std::vector snapshot_ts, - std::vector snapshots, - std::vector output_ts) +std::vector SnapshotInterpolator::interpolate(std::vector + snapshot_ts, + std::vector snapshots, + std::vector output_ts) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); @@ -49,7 +50,7 @@ std::vector SnapshotInterpolator::interpolate(std::vector snap for(int i = 0; i < n_out; ++i) { Vector* temp_snapshot =new Vector(snapshots[0]->dim(), - snapshots[0]->distributed()); + snapshots[0]->distributed()); output_snapshots.push_back(temp_snapshot); } @@ -73,7 +74,7 @@ std::vector SnapshotInterpolator::interpolate(std::vector snap if(sign(d_temp)!=sign(delta[0])) { d_temp = 0; - } + } else if (sign(delta[0]) != sign(delta[1]) && abs(d_temp) > abs(3*delta[0])) { d_temp = 3*delta[0]; @@ -87,20 +88,23 @@ std::vector SnapshotInterpolator::interpolate(std::vector snap while(output_ts[counter]->getData()[0] <= t_in[j+1]) { t = output_ts[counter]->getData()[0]; - output_snapshots[counter]->getData()[i] = y_in[j]*computeH1(t,t_in[j],t_in[j+1]) + - y_in[j+1]*computeH2(t,t_in[j],t_in[j+1]) + - d[j]*computeH3(t,t_in[j],t_in[j+1]) + - d[j+1]*computeH4(t,t_in[j],t_in[j+1]); + output_snapshots[counter]->getData()[i] = y_in[j]*computeH1(t,t_in[j], + t_in[j+1]) + + y_in[j+1]*computeH2(t,t_in[j],t_in[j+1]) + + d[j]*computeH3(t,t_in[j],t_in[j+1]) + + d[j+1]*computeH4(t,t_in[j],t_in[j+1]); counter += 1; } } - d_temp = ((2*h[n_snap-2]+h[n_snap-3])*delta[n_snap-2] - h[n_snap-2]*delta[n_snap-3])/(h[n_snap-2]+h[n_snap-3]); + d_temp = ((2*h[n_snap-2]+h[n_snap-3])*delta[n_snap-2] - h[n_snap-2]*delta[n_snap + -3])/(h[n_snap-2]+h[n_snap-3]); if(sign(d_temp) != sign(delta[n_snap-2])) { d_temp = 0; } - else if (sign(delta[n_snap-2]) != sign(delta[n_snap-3]) && abs(d_temp) > abs(3*delta[n_snap-2])) + else if (sign(delta[n_snap-2]) != sign(delta[n_snap-3]) + && abs(d_temp) > abs(3*delta[n_snap-2])) { d_temp = 3*delta[n_snap-2]; } @@ -108,20 +112,22 @@ std::vector SnapshotInterpolator::interpolate(std::vector snap while(counter < n_out && output_ts[counter]->getData()[0] <= t_in[n_snap-1]) { t = output_ts[counter]->getData()[0]; - output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t,t_in[n_snap-2],t_in[n_snap-1]) + - y_in[n_snap-1]*computeH2(t,t_in[n_snap-2],t_in[n_snap-1]) + - d[n_snap-2]*computeH3(t,t_in[n_snap-2],t_in[n_snap-1]) + - d[n_snap-1]*computeH4(t,t_in[n_snap-2],t_in[n_snap-1]); + output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t, + t_in[n_snap-2],t_in[n_snap-1]) + + y_in[n_snap-1]*computeH2(t,t_in[n_snap-2],t_in[n_snap-1]) + + d[n_snap-2]*computeH3(t,t_in[n_snap-2],t_in[n_snap-1]) + + d[n_snap-1]*computeH4(t,t_in[n_snap-2],t_in[n_snap-1]); counter += 1; } } return output_snapshots; } -std::vector SnapshotInterpolator::interpolate(std::vector snapshot_ts, - std::vector snapshots, - int n_out, - std::vector* output_ts) +std::vector SnapshotInterpolator::interpolate(std::vector + snapshot_ts, + std::vector snapshots, + int n_out, + std::vector* output_ts) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); @@ -137,18 +143,20 @@ std::vector SnapshotInterpolator::interpolate(std::vector snap double dt = (t_max-t_min)/(n_out-1); output_ts->clear(); - + for(int i = 0; i < n_out; ++i) { Vector* temp_t = new Vector(1, false); temp_t->getData()[0] = t_min + i*dt; output_ts->push_back(temp_t); } - std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap << std::endl; + std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap << + std::endl; return interpolate(snapshot_ts,snapshots,*output_ts); } -double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1, double h2) +double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1, + double h2) { double d = 0.0; double alpha = (h1 + 2*h2)/(3*(h1+h2)); From aa160e99d5e537cee86683d4952322d24de8b274 Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 13:45:06 -0700 Subject: [PATCH 05/17] Improved description of snapshot interpolator to include references --- lib/algo/manifold_interp/SnapshotInterpolator.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index c27e3c175..432d7ee4b 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -8,8 +8,12 @@ * *****************************************************************************/ -// Description: Implementation of the VectorInterpolator algorithm. - +/** + * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) + * + */ #include "SnapshotInterpolator.h" #include From a12a6456ebab5a1424b3726c0def24d8db7d3004 Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 14:06:26 -0700 Subject: [PATCH 06/17] Stylize --- lib/algo/manifold_interp/SnapshotInterpolator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index 432d7ee4b..d9aa3c873 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -10,7 +10,7 @@ /** * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE - * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ From 706876b2121a933132046460b577f3d012dddb25 Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 25 Apr 2024 14:21:04 -0700 Subject: [PATCH 07/17] Stylize --- .../manifold_interp/SnapshotInterpolator.h | 8 ++-- unit_tests/snapshot_interpolator_test.cpp | 46 ++++++++++--------- 2 files changed, 29 insertions(+), 25 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.h b/lib/algo/manifold_interp/SnapshotInterpolator.h index 943afb973..df2517fd6 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.h +++ b/lib/algo/manifold_interp/SnapshotInterpolator.h @@ -12,7 +12,7 @@ class Vector; /** * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE - * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ @@ -22,18 +22,18 @@ class SnapshotInterpolator SnapshotInterpolator(); ~SnapshotInterpolator(); - std::vector interpolate(std::vector snapshot_ts, + std::vector interpolate(std::vector snapshot_ts, std::vector snapshots, std::vector output_ts); - std::vector interpolate(std::vector snapshot_ts, + std::vector interpolate(std::vector snapshot_ts, std::vector snapshots, int n_out, std::vector* output_ts); std::vector interpolate(); - + private: diff --git a/unit_tests/snapshot_interpolator_test.cpp b/unit_tests/snapshot_interpolator_test.cpp index fcaa424fc..dd3386911 100644 --- a/unit_tests/snapshot_interpolator_test.cpp +++ b/unit_tests/snapshot_interpolator_test.cpp @@ -30,7 +30,7 @@ int main(int argc, char *argv[]) //input times vector t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15}; //Test function from original PCHIP paper - vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; + vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; //sin(x) vector y2{0, 0.909297426825682, 0.141120008059867, -0.958924274663138, -0.279415498198926, 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, @@ -38,11 +38,11 @@ int main(int argc, char *argv[]) //query points vector tq{0, 0.6000, 1.2000, 1.8000, 2.4000, 3.0000, 3.6000, 4.2000, 4.8000, 5.4000, 6.0000, 6.6000, - 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, + 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, 13.2000, 13.8000, 14.4000, 15.000}; //pchip reference solution for original test function. - vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, - 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, + vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, + 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, @@ -83,25 +83,27 @@ int main(int argc, char *argv[]) } SnapshotInterpolator* interp = new SnapshotInterpolator(); - + std::cout << "Beginning base interpolation function" << std::endl; out_snapshots = interp->interpolate(times,snapshots,out_times); std::cout << "Finished interpolation " << std::endl; -/* - for(int i = 0; i < out_snapshots.size(); ++i) - { - std::cout << "Time " << i << " is " << out_times[i]->item(0) << std::endl; - std::cout << "Reference at " << i << " is (" << reference_snapshots[i]->item(0) << - "," << reference_snapshots[i]->item(1) << ")" << std::endl; - std::cout << "Snapshot interpolation at " << i << " is (" << out_snapshots[i]->item(0) << - "," << out_snapshots[i]->item(1) << ")" << std::endl; - } - */ + /* + for(int i = 0; i < out_snapshots.size(); ++i) + { + std::cout << "Time " << i << " is " << out_times[i]->item(0) << std::endl; + std::cout << "Reference at " << i << " is (" << reference_snapshots[i]->item(0) << + "," << reference_snapshots[i]->item(1) << ")" << std::endl; + std::cout << "Snapshot interpolation at " << i << " is (" << out_snapshots[i]->item(0) << + "," << out_snapshots[i]->item(1) << ")" << std::endl; + } + */ for(int i = 0; i < out_snapshots.size(); ++i) { - std::cout << "Error at time[" << i << "] = " << out_times[i]->item(0) << " is (" << - reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," - << reference_snapshots[i]->item(1) - out_snapshots[i]->item(1) << ") " << std::endl; + std::cout << "Error at time[" << i << "] = " << out_times[i]->item( + 0) << " is (" << + reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," + << reference_snapshots[i]->item(1) - out_snapshots[i]->item( + 1) << ") " << std::endl; } @@ -111,8 +113,10 @@ int main(int argc, char *argv[]) for(int i = 0; i < out_snapshots.size(); ++i) { - std::cout << "Error at time[" << i << "] = " << out_times[i]->item(0) << " is (" << - reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," - << reference_snapshots[i]->item(1) - out_snapshots[i]->item(1) << ") " << std::endl; + std::cout << "Error at time[" << i << "] = " << out_times[i]->item( + 0) << " is (" << + reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," + << reference_snapshots[i]->item(1) - out_snapshots[i]->item( + 1) << ") " << std::endl; } } \ No newline at end of file From 841b47f140835167698488ad4400e4b128a0a727 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 7 May 2024 16:43:37 -0700 Subject: [PATCH 08/17] Addressed Dylan's Comments from PR#283. --- .../manifold_interp/SnapshotInterpolator.cpp | 65 +++++++++--------- .../manifold_interp/SnapshotInterpolator.h | 67 +++++++++++++------ unit_tests/snapshot_interpolator_test.cpp | 66 +++++++++++------- 3 files changed, 119 insertions(+), 79 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index d9aa3c873..bee830b71 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -10,12 +10,13 @@ /** * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE - * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ #include "SnapshotInterpolator.h" +#include #include #include #include "linalg/Vector.h" @@ -31,30 +32,24 @@ using namespace std; namespace CAROM { -SnapshotInterpolator::SnapshotInterpolator() -{ - -} - - -std::vector SnapshotInterpolator::interpolate(std::vector - snapshot_ts, - std::vector snapshots, - std::vector output_ts) +void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, + std::vector& snapshots, + std::vector& output_ts, + std::vector& output_snapshots) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); CAROM_VERIFY(output_ts.size() > 0); + CAROM_VERIFY(output_snapshots.size() == 0); int n_out = output_ts.size(); int n_snap = snapshots.size(); int n_dim = snapshots[0]->dim(); - std::vector output_snapshots; for(int i = 0; i < n_out; ++i) { - Vector* temp_snapshot =new Vector(snapshots[0]->dim(), - snapshots[0]->distributed()); + Vector* temp_snapshot = new Vector(snapshots[0]->dim(), + snapshots[0]->distributed()); output_snapshots.push_back(temp_snapshot); } @@ -113,7 +108,9 @@ std::vector SnapshotInterpolator::interpolate(std::vector d_temp = 3*delta[n_snap-2]; } d.push_back(d_temp); - while(counter < n_out && output_ts[counter]->getData()[0] <= t_in[n_snap-1]) + + while(counter < n_out + && output_ts[counter]->getData()[0] - t_in[n_snap-1] <= FLT_EPSILON ) { t = output_ts[counter]->getData()[0]; output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t, @@ -124,18 +121,20 @@ std::vector SnapshotInterpolator::interpolate(std::vector counter += 1; } } - return output_snapshots; } -std::vector SnapshotInterpolator::interpolate(std::vector - snapshot_ts, - std::vector snapshots, - int n_out, - std::vector* output_ts) +void SnapshotInterpolator::interpolate(std::vector& + snapshot_ts, + std::vector& snapshots, + int n_out, + std::vector& output_ts, + std::vector& output_snapshots) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); CAROM_VERIFY(n_out > 0); + CAROM_VERIFY(output_ts.size() == 0 && output_snapshots.size() == 0); + int n_snap = snapshots.size(); int n_dim = snapshots[0]->dim(); @@ -145,21 +144,21 @@ std::vector SnapshotInterpolator::interpolate(std::vector double t_min = snapshot_ts[0]->getData()[0]; double t_max = snapshot_ts[n_snap-1]->getData()[0]; double dt = (t_max-t_min)/(n_out-1); - output_ts->clear(); - + output_ts.clear(); for(int i = 0; i < n_out; ++i) { Vector* temp_t = new Vector(1, false); temp_t->getData()[0] = t_min + i*dt; - output_ts->push_back(temp_t); + output_ts.push_back(temp_t); } std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap << std::endl; - return interpolate(snapshot_ts,snapshots,*output_ts); + interpolate(snapshot_ts,snapshots,output_ts, output_snapshots); } -double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1, +const double SnapshotInterpolator::computeDerivative(double S1, double S2, + double h1, double h2) { double d = 0.0; @@ -175,41 +174,41 @@ double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1, return d; } -double SnapshotInterpolator::computeH1(double x, double xl, double xr) +const double SnapshotInterpolator::computeH1(double x, double xl, double xr) { double h = xr - xl; return computePhi((xr-x)/h); } -double SnapshotInterpolator::computeH2(double x, double xl, double xr) +const double SnapshotInterpolator::computeH2(double x, double xl, double xr) { double h = xr - xl; return computePhi((x-xl)/h); } -double SnapshotInterpolator::computeH3(double x, double xl, double xr) +const double SnapshotInterpolator::computeH3(double x, double xl, double xr) { double h = xr-xl; return -h*computePsi((xr-x)/h); } -double SnapshotInterpolator::computeH4(double x, double xl, double xr) +const double SnapshotInterpolator::computeH4(double x, double xl, double xr) { double h = xr-xl; return h*computePsi((x-xl)/h); } -double SnapshotInterpolator::computePhi(double t) +const double SnapshotInterpolator::computePhi(double t) { return 3.*pow(t,2.) - 2*pow(t,3.); } -double SnapshotInterpolator::computePsi(double t) +const double SnapshotInterpolator::computePsi(double t) { return pow(t,3.) - pow(t,2.); } -int SnapshotInterpolator::sign(double a) +const int SnapshotInterpolator::sign(double a) { double TOL = 1e-15; if(abs(a) < TOL)return 0; diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.h b/lib/algo/manifold_interp/SnapshotInterpolator.h index df2517fd6..ab272ab7a 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.h +++ b/lib/algo/manifold_interp/SnapshotInterpolator.h @@ -12,39 +12,66 @@ class Vector; /** * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE - * PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ class SnapshotInterpolator { public: - SnapshotInterpolator(); - ~SnapshotInterpolator(); - std::vector interpolate(std::vector snapshot_ts, - std::vector snapshots, - std::vector output_ts); - std::vector interpolate(std::vector snapshot_ts, - std::vector snapshots, - int n_out, - std::vector* output_ts); - - std::vector interpolate(); + SnapshotInterpolator() + {} + ~SnapshotInterpolator() + {} + /** + * @brief Compute new snapshots interpolated from snapshot_ts to + * output_ts. + * + * @param[in] snapshot_ts The parameter points. + * @param[in] snapshots The rotation matrices associated with + * each parameter point. + * @param[in] output_ts Requested times for interpolated + * snapshots + * @param[out] output_snapshots snapshots at output_ts interpolated + * from snapshot_ts + */ + void interpolate(std::vector& snapshot_ts, + std::vector& snapshots, + std::vector& output_ts, + std::vector&output_snapshots); + /** + * @brief Compute new snapshots interpolated from snapshot_ts to + * output_ts. + * + * @param[in] snapshot_ts The parameter points. + * @param[in] snapshots The rotation matrices associated with + * each parameter point. + * @param[in] n_out Number of output snapshots requested + * @param[out] output_ts std::vector of CAROM::Vectors that are + the times of the interpolated snapshots. + * @param[out] output_snapshots snapshots at output_ts interpolated + * from snapshot_ts + */ + void interpolate(std::vector& snapshot_ts, + std::vector& snapshots, + int n_out, + std::vector& output_ts, + std::vector& output_snapshots); private: - double computeDerivative(double S1, double S2, double h1, double h2); - double computeH1(double x, double xl, double xr); - double computeH2(double x, double xl, double xr); - double computeH3(double x, double xl, double xr); - double computeH4(double x, double xl, double xr); - double computePhi(double t); - double computePsi(double t); - int sign(double a); + const double computeDerivative(double S1, double S2, double h1, double h2); + const double computeH1(double x, double xl, double xr); + const double computeH2(double x, double xl, double xr); + const double computeH3(double x, double xl, double xr); + const double computeH4(double x, double xl, double xr); + const double computePhi(double t); + const double computePsi(double t); + const int sign(double a); }; } diff --git a/unit_tests/snapshot_interpolator_test.cpp b/unit_tests/snapshot_interpolator_test.cpp index dd3386911..2de21e589 100644 --- a/unit_tests/snapshot_interpolator_test.cpp +++ b/unit_tests/snapshot_interpolator_test.cpp @@ -10,6 +10,7 @@ #include "algo/manifold_interp/SnapshotInterpolator.h" #include "linalg/Vector.h" +#include #include #include @@ -27,6 +28,7 @@ int main(int argc, char *argv[]) { int n_out = 25; int n_snap = 11; + bool SUCCESS = true; //input times vector t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15}; //Test function from original PCHIP paper @@ -84,39 +86,51 @@ int main(int argc, char *argv[]) SnapshotInterpolator* interp = new SnapshotInterpolator(); - std::cout << "Beginning base interpolation function" << std::endl; - out_snapshots = interp->interpolate(times,snapshots,out_times); - std::cout << "Finished interpolation " << std::endl; - /* - for(int i = 0; i < out_snapshots.size(); ++i) - { - std::cout << "Time " << i << " is " << out_times[i]->item(0) << std::endl; - std::cout << "Reference at " << i << " is (" << reference_snapshots[i]->item(0) << - "," << reference_snapshots[i]->item(1) << ")" << std::endl; - std::cout << "Snapshot interpolation at " << i << " is (" << out_snapshots[i]->item(0) << - "," << out_snapshots[i]->item(1) << ")" << std::endl; - } - */ + + interp->interpolate(times,snapshots,out_times,out_snapshots); + for(int i = 0; i < out_snapshots.size(); ++i) { - std::cout << "Error at time[" << i << "] = " << out_times[i]->item( - 0) << " is (" << - reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," - << reference_snapshots[i]->item(1) - out_snapshots[i]->item( - 1) << ") " << std::endl; + if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( + 0)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + std::cout << "Test failed." << std::endl; + return SUCCESS; + } + else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( + 1)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + std::cout << "Test failed." << std::endl; + return SUCCESS; + } } - std::cout << "Beginning variant interpolation function" << std::endl; - out_snapshots = interp->interpolate(times,snapshots,26,&out_times); - std::cout << "Finished interpolation " << std::endl; + + out_snapshots.clear(); + out_times.clear(); + + interp->interpolate(times,snapshots,26,out_times,out_snapshots); for(int i = 0; i < out_snapshots.size(); ++i) { - std::cout << "Error at time[" << i << "] = " << out_times[i]->item( - 0) << " is (" << - reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << "," - << reference_snapshots[i]->item(1) - out_snapshots[i]->item( - 1) << ") " << std::endl; + if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( + 0)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + std::cout << "Test failed." << std::endl; + return SUCCESS; + } + else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( + 1)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + std::cout << "Test failed." << std::endl; + return SUCCESS; + } } + std::cout << "Test Succeeded" << std::endl; + return SUCCESS; } \ No newline at end of file From 52003a0a9c0576b92d6fe34a4bc3d7fcd82777b8 Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:26:20 -0700 Subject: [PATCH 09/17] Added snapshot interpolator unit test to CMakeLists --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f34a8e490..a967d30f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -282,7 +282,8 @@ if(GTEST_FOUND) IncrementalSVDBrand GreedyCustomSampler NNLS - basis_conversion) + basis_conversion + SnapshotInterpolator) foreach(stem IN LISTS unit_test_stems) add_executable(test_${stem} unit_tests/test_${stem}.cpp) target_link_libraries(test_${stem} PRIVATE ROM From 61ff4b4c5ca7490e86d77a2d52284e958ab0604c Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:26:54 -0700 Subject: [PATCH 10/17] Added snapshot_interpolation unit test --- unit_tests/test_SnapshotInterpolator.cpp | 148 +++++++++++++++++++++++ 1 file changed, 148 insertions(+) create mode 100644 unit_tests/test_SnapshotInterpolator.cpp diff --git a/unit_tests/test_SnapshotInterpolator.cpp b/unit_tests/test_SnapshotInterpolator.cpp new file mode 100644 index 000000000..cb0de429e --- /dev/null +++ b/unit_tests/test_SnapshotInterpolator.cpp @@ -0,0 +1,148 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2024, 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) + * + *****************************************************************************/ + +// Description: This source file is a test runner that uses the Google Test +// Framework to run unit tests on the CAROM::Matrix class. + +#include +#include + +#ifdef CAROM_HAS_GTEST +#include +#include "algo/manifold_interp/SnapshotInterpolator.h" +#include "linalg/Vector.h" +#include + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +TEST(InterpolationTest,test_accuracy) +{ + int n_out = 25; + int n_snap = 11; + bool SUCCESS = true; + //input times + std::vector t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15}; + //Test function from original PCHIP paper + std::vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; + //sin(x) + std::vector y2{0, 0.909297426825682, 0.141120008059867, -0.958924274663138, -0.279415498198926, + 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, + 0.990607355694870, 0.650287840157117}; + + //query points + std::vector tq{0, 0.6000, 1.2000, 1.8000, 2.4000, 3.0000, 3.6000, 4.2000, 4.8000, 5.4000, 6.0000, 6.6000, + 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, + 13.2000, 13.8000, 14.4000, 15.000}; + //pchip reference solution for original test function. + std::vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, + 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, + 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, + 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, + 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, + 50.000000000000000, 53.832363636363638, 55.720727272727267, 58.433818181818161, + 67.055999999999969, 85.000000000000000}; + //pchip reference solution for sin(x) + std::vector yq2{0, 0.569748887844739, 0.833039030477175, 0.906694689302138, 0.701592409322499, 0.141120008059867, + -0.288488198334352, -0.697095554949408, -0.939878053603591, -0.782971083698006, -0.279415498198926, + 0.188293444380395, 0.669217685146469, 0.965688937709033, 0.846474732609843, 0.412118485241757, + -0.077580693288345, -0.623537711025344, -0.971758328554163, -0.890773577229575, -0.536572918000435, + -0.041614069121016, 0.560852411851254, 0.957953731078007, 0.938810668593539, 0.650287840157117}; + + std::vector snapshots; + std::vector out_snapshots; + std::vector reference_snapshots; + std::vector times; + std::vector out_times; + for(int i = 0; i < t.size(); ++i) + { + CAROM::Vector* temp_t = new CAROM::Vector(1, false); + temp_t->item(0) = t[i]; + times.push_back(temp_t); + CAROM::Vector* temp_y = new CAROM::Vector(2,false); + temp_y->item(0) = y[i]; + temp_y->item(1) = y2[i]; + snapshots.push_back(temp_y); + } + + for(int i = 0; i < tq.size(); ++i) + { + CAROM::Vector* temp_t = new CAROM::Vector(1, false); + temp_t->item(0) = tq[i]; + out_times.push_back(temp_t); + CAROM::Vector* temp_y = new CAROM::Vector(2,false); + temp_y->item(0) = yq1[i]; + temp_y->item(1) = yq2[i]; + reference_snapshots.push_back(temp_y); + } + + CAROM::SnapshotInterpolator* interp = new CAROM::SnapshotInterpolator(); + + + interp->interpolate(times,snapshots,out_times,out_snapshots); + + for(int i = 0; i < out_snapshots.size(); ++i) + { + if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( + 0)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + break; + } + else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( + 1)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + break; + } + } + + out_snapshots.clear(); + out_times.clear(); + + interp->interpolate(times,snapshots,26,out_times,out_snapshots); + + for(int i = 0; i < out_snapshots.size(); ++i) + { + if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( + 0)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + break; + } + else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( + 1)) > 10.*FLT_EPSILON) + { + SUCCESS = false; + break; + } + } + EXPECT_TRUE(SUCCESS); +} + + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + return result; +} +#else // #ifndef CAROM_HAS_GTEST +int main() +{ + std::cout << "libROM was compiled without Google Test support, so unit " + << "tests have been disabled. To enable unit tests, compile " + << "libROM with Google Test support." << std::endl; +} +#endif // #endif CAROM_HAS_GTEST From 21afbeaf8d10529bbd96f00ffd3e4e43cdb17f6c Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:28:40 -0700 Subject: [PATCH 11/17] Removed old snapshot interpolator test --- unit_tests/snapshot_interpolator_test.cpp | 136 ---------------------- 1 file changed, 136 deletions(-) delete mode 100644 unit_tests/snapshot_interpolator_test.cpp diff --git a/unit_tests/snapshot_interpolator_test.cpp b/unit_tests/snapshot_interpolator_test.cpp deleted file mode 100644 index 2de21e589..000000000 --- a/unit_tests/snapshot_interpolator_test.cpp +++ /dev/null @@ -1,136 +0,0 @@ -/****************************************************************************** - * - * Copyright (c) 2013-2024, 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) - * - *****************************************************************************/ - -#include "algo/manifold_interp/SnapshotInterpolator.h" -#include "linalg/Vector.h" -#include -#include -#include - -#ifndef _WIN32 -#include // mkdir -#else -#include // _mkdir -#define mkdir(dir, mode) _mkdir(dir) -#endif - -using namespace std; -using namespace CAROM; - -int main(int argc, char *argv[]) -{ - int n_out = 25; - int n_snap = 11; - bool SUCCESS = true; - //input times - vector t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15}; - //Test function from original PCHIP paper - vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; - //sin(x) - vector y2{0, 0.909297426825682, 0.141120008059867, -0.958924274663138, -0.279415498198926, - 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, - 0.990607355694870, 0.650287840157117}; - - //query points - vector tq{0, 0.6000, 1.2000, 1.8000, 2.4000, 3.0000, 3.6000, 4.2000, 4.8000, 5.4000, 6.0000, 6.6000, - 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, - 13.2000, 13.8000, 14.4000, 15.000}; - //pchip reference solution for original test function. - vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, - 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, - 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, - 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, - 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, - 50.000000000000000, 53.832363636363638, 55.720727272727267, 58.433818181818161, - 67.055999999999969, 85.000000000000000}; - //pchip reference solution for sin(x) - vector yq2{0, 0.569748887844739, 0.833039030477175, 0.906694689302138, 0.701592409322499, 0.141120008059867, - -0.288488198334352, -0.697095554949408, -0.939878053603591, -0.782971083698006, -0.279415498198926, - 0.188293444380395, 0.669217685146469, 0.965688937709033, 0.846474732609843, 0.412118485241757, - -0.077580693288345, -0.623537711025344, -0.971758328554163, -0.890773577229575, -0.536572918000435, - -0.041614069121016, 0.560852411851254, 0.957953731078007, 0.938810668593539, 0.650287840157117}; - - std::vector snapshots; - std::vector out_snapshots; - std::vector reference_snapshots; - std::vector times; - std::vector out_times; - for(int i = 0; i < t.size(); ++i) - { - Vector* temp_t = new Vector(1, false); - temp_t->item(0) = t[i]; - times.push_back(temp_t); - Vector* temp_y = new Vector(2,false); - temp_y->item(0) = y[i]; - temp_y->item(1) = y2[i]; - snapshots.push_back(temp_y); - } - - for(int i = 0; i < tq.size(); ++i) - { - Vector* temp_t = new Vector(1, false); - temp_t->item(0) = tq[i]; - out_times.push_back(temp_t); - Vector* temp_y = new Vector(2,false); - temp_y->item(0) = yq1[i]; - temp_y->item(1) = yq2[i]; - reference_snapshots.push_back(temp_y); - } - - SnapshotInterpolator* interp = new SnapshotInterpolator(); - - - interp->interpolate(times,snapshots,out_times,out_snapshots); - - for(int i = 0; i < out_snapshots.size(); ++i) - { - if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( - 0)) > 10.*FLT_EPSILON) - { - SUCCESS = false; - std::cout << "Test failed." << std::endl; - return SUCCESS; - } - else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( - 1)) > 10.*FLT_EPSILON) - { - SUCCESS = false; - std::cout << "Test failed." << std::endl; - return SUCCESS; - } - } - - - - out_snapshots.clear(); - out_times.clear(); - - interp->interpolate(times,snapshots,26,out_times,out_snapshots); - - for(int i = 0; i < out_snapshots.size(); ++i) - { - if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item( - 0)) > 10.*FLT_EPSILON) - { - SUCCESS = false; - std::cout << "Test failed." << std::endl; - return SUCCESS; - } - else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item( - 1)) > 10.*FLT_EPSILON) - { - SUCCESS = false; - std::cout << "Test failed." << std::endl; - return SUCCESS; - } - } - std::cout << "Test Succeeded" << std::endl; - return SUCCESS; -} \ No newline at end of file From 712f4e1319628a2b23d3c77d44666bee68fe4282 Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:48:29 -0700 Subject: [PATCH 12/17] Added improvements for floating point comparison --- lib/algo/manifold_interp/SnapshotInterpolator.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index bee830b71..ef0e61d64 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -41,6 +41,8 @@ void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, CAROM_VERIFY(snapshot_ts.size() > 0); CAROM_VERIFY(output_ts.size() > 0); CAROM_VERIFY(output_snapshots.size() == 0); + CAROM_VERIFY(snapshot_ts[0]->getData()[0] - FLT_EPSILON <= output_ts[0]->getData()[0] && + output_ts[output_ts.size()-1]->getData()[0] <= snapshot_ts[snapshot_ts.size()-1]->getData()[0] + FLT_EPSILON); int n_out = output_ts.size(); int n_snap = snapshots.size(); @@ -110,7 +112,7 @@ void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, d.push_back(d_temp); while(counter < n_out - && output_ts[counter]->getData()[0] - t_in[n_snap-1] <= FLT_EPSILON ) + && output_ts[counter]->getData()[0] <= t_in[n_snap-1] + FLT_EPSILON ) { t = output_ts[counter]->getData()[0]; output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t, @@ -152,8 +154,6 @@ void SnapshotInterpolator::interpolate(std::vector& temp_t->getData()[0] = t_min + i*dt; output_ts.push_back(temp_t); } - std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap << - std::endl; interpolate(snapshot_ts,snapshots,output_ts, output_snapshots); } From 98197e6d26a1a280fb8dd3c86c6afbaa9d0e1cca Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:51:11 -0700 Subject: [PATCH 13/17] Stylize --- .../manifold_interp/SnapshotInterpolator.cpp | 6 ++-- unit_tests/test_SnapshotInterpolator.cpp | 28 +++++++++---------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/SnapshotInterpolator.cpp index ef0e61d64..dd4bfb085 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/SnapshotInterpolator.cpp @@ -41,8 +41,10 @@ void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, CAROM_VERIFY(snapshot_ts.size() > 0); CAROM_VERIFY(output_ts.size() > 0); CAROM_VERIFY(output_snapshots.size() == 0); - CAROM_VERIFY(snapshot_ts[0]->getData()[0] - FLT_EPSILON <= output_ts[0]->getData()[0] && - output_ts[output_ts.size()-1]->getData()[0] <= snapshot_ts[snapshot_ts.size()-1]->getData()[0] + FLT_EPSILON); + CAROM_VERIFY(snapshot_ts[0]->getData()[0] - FLT_EPSILON <= + output_ts[0]->getData()[0] && + output_ts[output_ts.size()-1]->getData()[0] <= snapshot_ts[snapshot_ts.size() + -1]->getData()[0] + FLT_EPSILON); int n_out = output_ts.size(); int n_snap = snapshots.size(); diff --git a/unit_tests/test_SnapshotInterpolator.cpp b/unit_tests/test_SnapshotInterpolator.cpp index cb0de429e..377ce4610 100644 --- a/unit_tests/test_SnapshotInterpolator.cpp +++ b/unit_tests/test_SnapshotInterpolator.cpp @@ -38,27 +38,27 @@ TEST(InterpolationTest,test_accuracy) std::vector y{10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85}; //sin(x) std::vector y2{0, 0.909297426825682, 0.141120008059867, -0.958924274663138, -0.279415498198926, - 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, - 0.990607355694870, 0.650287840157117}; + 0.989358246623382, 0.412118485241757, -0.999990206550703, -0.536572918000435, + 0.990607355694870, 0.650287840157117}; //query points std::vector tq{0, 0.6000, 1.2000, 1.8000, 2.4000, 3.0000, 3.6000, 4.2000, 4.8000, 5.4000, 6.0000, 6.6000, - 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, - 13.2000, 13.8000, 14.4000, 15.000}; + 7.2000, 7.8000, 8.4000, 9.0000, 9.6000, 10.2000, 10.8000, 11.4000, 12.0000, 12.6000, + 13.2000, 13.8000, 14.4000, 15.000}; //pchip reference solution for original test function. std::vector yq1{10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000004, - 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, - 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, - 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, - 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, - 50.000000000000000, 53.832363636363638, 55.720727272727267, 58.433818181818161, - 67.055999999999969, 85.000000000000000}; + 10.000000000000000, 10.000000000000000, 10.000000000000002, 10.000000000000000, + 10.000000000000004, 10.000000000000000, 10.000000000000000, 10.000000000000002, + 10.000000000000000, 10.000000000000004, 10.102641509433962, 10.500000000000000, + 11.106230625292378, 12.213163262123807, 14.128630750038976, 27.078413223140512, + 50.000000000000000, 53.832363636363638, 55.720727272727267, 58.433818181818161, + 67.055999999999969, 85.000000000000000}; //pchip reference solution for sin(x) std::vector yq2{0, 0.569748887844739, 0.833039030477175, 0.906694689302138, 0.701592409322499, 0.141120008059867, - -0.288488198334352, -0.697095554949408, -0.939878053603591, -0.782971083698006, -0.279415498198926, - 0.188293444380395, 0.669217685146469, 0.965688937709033, 0.846474732609843, 0.412118485241757, - -0.077580693288345, -0.623537711025344, -0.971758328554163, -0.890773577229575, -0.536572918000435, - -0.041614069121016, 0.560852411851254, 0.957953731078007, 0.938810668593539, 0.650287840157117}; + -0.288488198334352, -0.697095554949408, -0.939878053603591, -0.782971083698006, -0.279415498198926, + 0.188293444380395, 0.669217685146469, 0.965688937709033, 0.846474732609843, 0.412118485241757, + -0.077580693288345, -0.623537711025344, -0.971758328554163, -0.890773577229575, -0.536572918000435, + -0.041614069121016, 0.560852411851254, 0.957953731078007, 0.938810668593539, 0.650287840157117}; std::vector snapshots; std::vector out_snapshots; From 86fa3656f3c7158ca4988255587f221dcf93ab22 Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 20 May 2024 22:53:11 -0700 Subject: [PATCH 14/17] removed old snapshotinterpolator unit test from CMakeLists --- CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a967d30f3..b5424af3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,8 +228,7 @@ if (ENABLE_EXAMPLES) weak_scaling random_test smoke_static - load_samples - snapshot_interpolator_test) + load_samples) if (USE_MFEM) set(regression_test_names From a1420d79e1a2b600fcd3de93e787de9b51ad34bf Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 3 Jun 2024 16:35:32 -0700 Subject: [PATCH 15/17] Renamed SnapshotInterpolator to PCHIPInterpolator. Addressed relevant PR Comments. --- CMakeLists.txt | 2 +- lib/CMakeLists.txt | 2 +- ...Interpolator.cpp => PCHIPInterpolator.cpp} | 41 ++++++++++--------- ...shotInterpolator.h => PCHIPInterpolator.h} | 29 +++++++------ ...polator.cpp => test_PCHIPInterpolator.cpp} | 4 +- 5 files changed, 39 insertions(+), 39 deletions(-) rename lib/algo/manifold_interp/{SnapshotInterpolator.cpp => PCHIPInterpolator.cpp} (82%) rename lib/algo/manifold_interp/{SnapshotInterpolator.h => PCHIPInterpolator.h} (73%) rename unit_tests/{test_SnapshotInterpolator.cpp => test_PCHIPInterpolator.cpp} (97%) diff --git a/CMakeLists.txt b/CMakeLists.txt index b5424af3f..24aaee036 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -282,7 +282,7 @@ if(GTEST_FOUND) GreedyCustomSampler NNLS basis_conversion - SnapshotInterpolator) + PCHIPInterpolator) foreach(stem IN LISTS unit_test_stems) add_executable(test_${stem} unit_tests/test_${stem}.cpp) target_link_libraries(test_${stem} PRIVATE ROM diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index ac3e66a8e..6b52a53fb 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -45,7 +45,7 @@ set(module_list algo/manifold_interp/Interpolator algo/manifold_interp/MatrixInterpolator algo/manifold_interp/VectorInterpolator - algo/manifold_interp/SnapshotInterpolator + algo/manifold_interp/PCHIPInterpolator hyperreduction/DEIM hyperreduction/GNAT hyperreduction/QDEIM diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.cpp b/lib/algo/manifold_interp/PCHIPInterpolator.cpp similarity index 82% rename from lib/algo/manifold_interp/SnapshotInterpolator.cpp rename to lib/algo/manifold_interp/PCHIPInterpolator.cpp index dd4bfb085..096565748 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.cpp +++ b/lib/algo/manifold_interp/PCHIPInterpolator.cpp @@ -14,7 +14,7 @@ * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ -#include "SnapshotInterpolator.h" +#include "PCHIPInterpolator.h" #include #include @@ -32,10 +32,10 @@ using namespace std; namespace CAROM { -void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, - std::vector& snapshots, - std::vector& output_ts, - std::vector& output_snapshots) +void PCHIPInterpolator::interpolate(std::vector& snapshot_ts, + std::vector& snapshots, + std::vector& output_ts, + std::vector& output_snapshots) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); @@ -124,15 +124,16 @@ void SnapshotInterpolator::interpolate(std::vector& snapshot_ts, d[n_snap-1]*computeH4(t,t_in[n_snap-2],t_in[n_snap-1]); counter += 1; } + CAROM_VERIFY(counter == n_out); } } -void SnapshotInterpolator::interpolate(std::vector& - snapshot_ts, - std::vector& snapshots, - int n_out, - std::vector& output_ts, - std::vector& output_snapshots) +void PCHIPInterpolator::interpolate(std::vector& + snapshot_ts, + std::vector& snapshots, + int n_out, + std::vector& output_ts, + std::vector& output_snapshots) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); @@ -159,9 +160,9 @@ void SnapshotInterpolator::interpolate(std::vector& interpolate(snapshot_ts,snapshots,output_ts, output_snapshots); } -const double SnapshotInterpolator::computeDerivative(double S1, double S2, +double PCHIPInterpolator::computeDerivative(double S1, double S2, double h1, - double h2) + double h2) const { double d = 0.0; double alpha = (h1 + 2*h2)/(3*(h1+h2)); @@ -176,41 +177,41 @@ const double SnapshotInterpolator::computeDerivative(double S1, double S2, return d; } -const double SnapshotInterpolator::computeH1(double x, double xl, double xr) +double PCHIPInterpolator::computeH1(double x, double xl, double xr) const { double h = xr - xl; return computePhi((xr-x)/h); } -const double SnapshotInterpolator::computeH2(double x, double xl, double xr) +double PCHIPInterpolator::computeH2(double x, double xl, double xr) const { double h = xr - xl; return computePhi((x-xl)/h); } -const double SnapshotInterpolator::computeH3(double x, double xl, double xr) +double PCHIPInterpolator::computeH3(double x, double xl, double xr) const { double h = xr-xl; return -h*computePsi((xr-x)/h); } -const double SnapshotInterpolator::computeH4(double x, double xl, double xr) +double PCHIPInterpolator::computeH4(double x, double xl, double xr) const { double h = xr-xl; return h*computePsi((x-xl)/h); } -const double SnapshotInterpolator::computePhi(double t) +double PCHIPInterpolator::computePhi(double t) const { return 3.*pow(t,2.) - 2*pow(t,3.); } -const double SnapshotInterpolator::computePsi(double t) +double PCHIPInterpolator::computePsi(double t) const { return pow(t,3.) - pow(t,2.); } -const int SnapshotInterpolator::sign(double a) +int PCHIPInterpolator::sign(double a) const { double TOL = 1e-15; if(abs(a) < TOL)return 0; diff --git a/lib/algo/manifold_interp/SnapshotInterpolator.h b/lib/algo/manifold_interp/PCHIPInterpolator.h similarity index 73% rename from lib/algo/manifold_interp/SnapshotInterpolator.h rename to lib/algo/manifold_interp/PCHIPInterpolator.h index ab272ab7a..4fd2e7014 100644 --- a/lib/algo/manifold_interp/SnapshotInterpolator.h +++ b/lib/algo/manifold_interp/PCHIPInterpolator.h @@ -1,7 +1,6 @@ -#ifndef included_SnapshotInterpolator_h -#define included_SnapshotInterpolator_h +#ifndef included_PCHIPInterpolator_h +#define included_PCHIPInterpolator_h -#include "Interpolator.h" #include #include @@ -11,19 +10,19 @@ class Matrix; class Vector; /** - * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE + * Implements PCHIP algorithm. Based on "A METHOD FOR UCTING LOCAL MONOTONE * PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) * */ -class SnapshotInterpolator +class PCHIPInterpolator { public: - SnapshotInterpolator() + PCHIPInterpolator() {} - ~SnapshotInterpolator() + ~PCHIPInterpolator() {} /** @@ -64,14 +63,14 @@ class SnapshotInterpolator private: - const double computeDerivative(double S1, double S2, double h1, double h2); - const double computeH1(double x, double xl, double xr); - const double computeH2(double x, double xl, double xr); - const double computeH3(double x, double xl, double xr); - const double computeH4(double x, double xl, double xr); - const double computePhi(double t); - const double computePsi(double t); - const int sign(double a); + double computeDerivative(double S1, double S2, double h1, double h2) const; + double computeH1(double x, double xl, double xr) const; + double computeH2(double x, double xl, double xr) const; + double computeH3(double x, double xl, double xr) const; + double computeH4(double x, double xl, double xr) const; + double computePhi(double t) const; + double computePsi(double t) const; + int sign(double a) const; }; } diff --git a/unit_tests/test_SnapshotInterpolator.cpp b/unit_tests/test_PCHIPInterpolator.cpp similarity index 97% rename from unit_tests/test_SnapshotInterpolator.cpp rename to unit_tests/test_PCHIPInterpolator.cpp index 377ce4610..2d2b2d15d 100644 --- a/unit_tests/test_SnapshotInterpolator.cpp +++ b/unit_tests/test_PCHIPInterpolator.cpp @@ -16,7 +16,7 @@ #ifdef CAROM_HAS_GTEST #include -#include "algo/manifold_interp/SnapshotInterpolator.h" +#include "algo/manifold_interp/PCHIPInterpolator.h" #include "linalg/Vector.h" #include @@ -87,7 +87,7 @@ TEST(InterpolationTest,test_accuracy) reference_snapshots.push_back(temp_y); } - CAROM::SnapshotInterpolator* interp = new CAROM::SnapshotInterpolator(); + CAROM::PCHIPInterpolator* interp = new CAROM::PCHIPInterpolator(); interp->interpolate(times,snapshots,out_times,out_snapshots); From f2a9a79ed5e230cd9712a218a34fe805a8c464e0 Mon Sep 17 00:00:00 2001 From: ptranq Date: Mon, 3 Jun 2024 16:36:51 -0700 Subject: [PATCH 16/17] Removed redundant else statement. --- lib/algo/manifold_interp/PCHIPInterpolator.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/lib/algo/manifold_interp/PCHIPInterpolator.cpp b/lib/algo/manifold_interp/PCHIPInterpolator.cpp index 096565748..d5ff5e0f5 100644 --- a/lib/algo/manifold_interp/PCHIPInterpolator.cpp +++ b/lib/algo/manifold_interp/PCHIPInterpolator.cpp @@ -170,10 +170,6 @@ double PCHIPInterpolator::computeDerivative(double S1, double S2, { d = S1*S2/(alpha*S2 + (1-alpha)*S1); } - else - { - d = 0.0; - } return d; } From 791b73f4f2c8b772f690ec34c52f012c46e6672f Mon Sep 17 00:00:00 2001 From: ptranq Date: Thu, 6 Jun 2024 21:42:45 -0700 Subject: [PATCH 17/17] Addressed Dylan's and Cole's comments. Added several CAROM_VERIFY statements to guarantee strictly increasing input / output ts, that there are >2 input ts, and that there are >1 output ts. --- .../manifold_interp/PCHIPInterpolator.cpp | 29 ++++++++++--------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/lib/algo/manifold_interp/PCHIPInterpolator.cpp b/lib/algo/manifold_interp/PCHIPInterpolator.cpp index d5ff5e0f5..dc0191a32 100644 --- a/lib/algo/manifold_interp/PCHIPInterpolator.cpp +++ b/lib/algo/manifold_interp/PCHIPInterpolator.cpp @@ -38,13 +38,19 @@ void PCHIPInterpolator::interpolate(std::vector& snapshot_ts, std::vector& output_snapshots) { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); - CAROM_VERIFY(snapshot_ts.size() > 0); - CAROM_VERIFY(output_ts.size() > 0); + CAROM_VERIFY(snapshot_ts.size() > 2); + CAROM_VERIFY(output_ts.size() > 1); CAROM_VERIFY(output_snapshots.size() == 0); CAROM_VERIFY(snapshot_ts[0]->getData()[0] - FLT_EPSILON <= output_ts[0]->getData()[0] && output_ts[output_ts.size()-1]->getData()[0] <= snapshot_ts[snapshot_ts.size() -1]->getData()[0] + FLT_EPSILON); + for(int i = 1; i < snapshot_ts.size(); ++i) + { + CAROM_VERIFY(snapshots[i-1]->dim() == snapshots[i]->dim()); + CAROM_VERIFY(snapshot_ts[i-1]->getData()[0] < snapshot_ts[i]->getData()[0]); + CAROM_VERIFY(output_ts[i-1]->getData()[0] < output_ts[i]->getData()[0]); + } int n_out = output_ts.size(); int n_snap = snapshots.size(); @@ -96,7 +102,7 @@ void PCHIPInterpolator::interpolate(std::vector& snapshot_ts, y_in[j+1]*computeH2(t,t_in[j],t_in[j+1]) + d[j]*computeH3(t,t_in[j],t_in[j+1]) + d[j+1]*computeH4(t,t_in[j],t_in[j+1]); - counter += 1; + counter++; } } @@ -122,7 +128,7 @@ void PCHIPInterpolator::interpolate(std::vector& snapshot_ts, y_in[n_snap-1]*computeH2(t,t_in[n_snap-2],t_in[n_snap-1]) + d[n_snap-2]*computeH3(t,t_in[n_snap-2],t_in[n_snap-1]) + d[n_snap-1]*computeH4(t,t_in[n_snap-2],t_in[n_snap-1]); - counter += 1; + counter++; } CAROM_VERIFY(counter == n_out); } @@ -137,15 +143,12 @@ void PCHIPInterpolator::interpolate(std::vector& { CAROM_VERIFY(snapshot_ts.size() == snapshots.size()); CAROM_VERIFY(snapshot_ts.size() > 0); - CAROM_VERIFY(n_out > 0); + CAROM_VERIFY(n_out > 2); CAROM_VERIFY(output_ts.size() == 0 && output_snapshots.size() == 0); - int n_snap = snapshots.size(); int n_dim = snapshots[0]->dim(); - - double t_min = snapshot_ts[0]->getData()[0]; double t_max = snapshot_ts[n_snap-1]->getData()[0]; double dt = (t_max-t_min)/(n_out-1); @@ -175,25 +178,25 @@ double PCHIPInterpolator::computeDerivative(double S1, double S2, double PCHIPInterpolator::computeH1(double x, double xl, double xr) const { - double h = xr - xl; + const double h = xr - xl; return computePhi((xr-x)/h); } double PCHIPInterpolator::computeH2(double x, double xl, double xr) const { - double h = xr - xl; + const double h = xr - xl; return computePhi((x-xl)/h); } double PCHIPInterpolator::computeH3(double x, double xl, double xr) const { - double h = xr-xl; + const double h = xr-xl; return -h*computePsi((xr-x)/h); } double PCHIPInterpolator::computeH4(double x, double xl, double xr) const { - double h = xr-xl; + const double h = xr-xl; return h*computePsi((x-xl)/h); } @@ -209,7 +212,7 @@ double PCHIPInterpolator::computePsi(double t) const int PCHIPInterpolator::sign(double a) const { - double TOL = 1e-15; + constexpr double TOL = 1e-15; if(abs(a) < TOL)return 0; else if(a > 0) return 1; else if (a < 0) return -1;