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