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

@@ -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

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

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 <cfloat>
+#include <limits.h>
+#include <cmath>
+#include "linalg/Vector.h"
+#include "mpi.h"
+
+/* Use C++11 built-in shared pointers if available; else fallback to Boost. */
+#if __cplusplus >= 201103L
+#include <memory>
+#else
+#include <boost/shared_ptr.hpp>
+#endif
+using namespace std;
+
+namespace CAROM {
+
+void PCHIPInterpolator::interpolate(std::vector<Vector*>& snapshot_ts,
+                                    std::vector<Vector*>& snapshots,
+                                    std::vector<Vector*>& output_ts,
+                                    std::vector<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<double> 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<Vector*>&
+                                    snapshot_ts,
+                                    std::vector<Vector*>& snapshots,
+                                    int n_out,
+                                    std::vector<Vector*>& output_ts,
+                                    std::vector<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;
+}
+
+}

lib/algo/manifold_interp/PCHIPInterpolator.h

@@ -0,0 +1,77 @@
+#ifndef included_PCHIPInterpolator_h
+#define included_PCHIPInterpolator_h
+
+#include <vector>
+#include <string>
+
+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<Vector*>& snapshot_ts,
+                     std::vector<Vector*>& snapshots,
+                     std::vector<Vector*>& output_ts,
+                     std::vector<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<Vector*>& snapshot_ts,
+                     std::vector<Vector*>& snapshots,
+                     int n_out,
+                     std::vector<Vector*>& output_ts,
+                     std::vector<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