From 884d5ae8f1ad3073cbb266193755d63090024322 Mon Sep 17 00:00:00 2001 From: ptranq Date: Wed, 12 Jun 2024 08:56:32 -0700 Subject: [PATCH] Snapshot interpolation (#283) * 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. * Reduced scope to only the interpolator and a unit test. Will construct SnapshotDMD as a separat PR. * Fixed unit test compilation. * Stylize * Improved description of snapshot interpolator to include references * Stylize * Stylize * Addressed Dylan's Comments from PR#283. * Added snapshot interpolator unit test to CMakeLists * Added snapshot_interpolation unit test * Removed old snapshot interpolator test * Added improvements for floating point comparison * Stylize * removed old snapshotinterpolator unit test from CMakeLists * Renamed SnapshotInterpolator to PCHIPInterpolator. Addressed relevant PR Comments. * Removed redundant else statement. * 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. --- CMakeLists.txt | 3 +- lib/CMakeLists.txt | 1 + .../manifold_interp/PCHIPInterpolator.cpp | 222 ++++++++++++++++++ lib/algo/manifold_interp/PCHIPInterpolator.h | 77 ++++++ unit_tests/test_PCHIPInterpolator.cpp | 148 ++++++++++++ 5 files changed, 450 insertions(+), 1 deletion(-) create mode 100644 lib/algo/manifold_interp/PCHIPInterpolator.cpp create mode 100644 lib/algo/manifold_interp/PCHIPInterpolator.h create mode 100644 unit_tests/test_PCHIPInterpolator.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 6fe7260b6..e2506dd04 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -283,7 +283,8 @@ if(GTEST_FOUND) IncrementalSVDBrand GreedyCustomSampler NNLS - basis_conversion) + basis_conversion + 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 6a85038d2..23343006e 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/PCHIPInterpolator hyperreduction/DEIM hyperreduction/GNAT hyperreduction/QDEIM diff --git a/lib/algo/manifold_interp/PCHIPInterpolator.cpp b/lib/algo/manifold_interp/PCHIPInterpolator.cpp new file mode 100644 index 000000000..dc0191a32 --- /dev/null +++ b/lib/algo/manifold_interp/PCHIPInterpolator.cpp @@ -0,0 +1,222 @@ +/****************************************************************************** + * + * 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) + * + *****************************************************************************/ + +/** + * Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE + * PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE + * PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980) + * + */ +#include "PCHIPInterpolator.h" + +#include +#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 { + +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() > 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(); + int n_dim = snapshots[0]->dim(); + + 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++; + } + } + + 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] + FLT_EPSILON ) + { + 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++; + } + CAROM_VERIFY(counter == n_out); + } +} + +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); + 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); + 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); + } + interpolate(snapshot_ts,snapshots,output_ts, output_snapshots); +} + +double PCHIPInterpolator::computeDerivative(double S1, double S2, + double h1, + double h2) const +{ + double d = 0.0; + double alpha = (h1 + 2*h2)/(3*(h1+h2)); + if(S1*S2 > 0) + { + d = S1*S2/(alpha*S2 + (1-alpha)*S1); + } + return d; +} + +double PCHIPInterpolator::computeH1(double x, double xl, double xr) const +{ + const double h = xr - xl; + return computePhi((xr-x)/h); +} + +double PCHIPInterpolator::computeH2(double x, double xl, double xr) const +{ + const double h = xr - xl; + return computePhi((x-xl)/h); +} + +double PCHIPInterpolator::computeH3(double x, double xl, double xr) const +{ + const double h = xr-xl; + return -h*computePsi((xr-x)/h); +} + +double PCHIPInterpolator::computeH4(double x, double xl, double xr) const +{ + const double h = xr-xl; + return h*computePsi((x-xl)/h); +} + +double PCHIPInterpolator::computePhi(double t) const +{ + return 3.*pow(t,2.) - 2*pow(t,3.); +} + +double PCHIPInterpolator::computePsi(double t) const +{ + return pow(t,3.) - pow(t,2.); +} + +int PCHIPInterpolator::sign(double a) const +{ + constexpr 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/PCHIPInterpolator.h b/lib/algo/manifold_interp/PCHIPInterpolator.h new file mode 100644 index 000000000..4fd2e7014 --- /dev/null +++ b/lib/algo/manifold_interp/PCHIPInterpolator.h @@ -0,0 +1,77 @@ +#ifndef included_PCHIPInterpolator_h +#define included_PCHIPInterpolator_h + +#include +#include + +namespace CAROM { + +class Matrix; +class Vector; + +/** + * 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 PCHIPInterpolator +{ +public: + + + PCHIPInterpolator() + {} + ~PCHIPInterpolator() + {} + + /** + * @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) 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; +}; + +} +#endif \ No newline at end of file diff --git a/unit_tests/test_PCHIPInterpolator.cpp b/unit_tests/test_PCHIPInterpolator.cpp new file mode 100644 index 000000000..2d2b2d15d --- /dev/null +++ b/unit_tests/test_PCHIPInterpolator.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/PCHIPInterpolator.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::PCHIPInterpolator* interp = new CAROM::PCHIPInterpolator(); + + + 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