Interpolation gradients

lib/algo/manifold_interp/Interpolator.cpp

@@ -34,7 +34,8 @@ Interpolator::Interpolator(const std::vector<Vector> & parameter_points,
                            int ref_point,
                            std::string rbf,
                            std::string interp_method,
-                           double closest_rbf_val)
+                           double closest_rbf_val,
+                           bool compute_gradients)
 {
     CAROM_VERIFY(parameter_points.size() == rotation_matrices.size());
     CAROM_VERIFY(parameter_points.size() > 1);
@@ -53,6 +54,7 @@ Interpolator::Interpolator(const std::vector<Vector> & parameter_points,
     MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
     MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs);
 
+    d_compute_gradients = compute_gradients;
     d_parameter_points = parameter_points;
     d_rotation_matrices = rotation_matrices;
     d_ref_point = ref_point;
@@ -131,6 +133,60 @@ std::vector<double> obtainRBFToTrainingPoints(
     return rbfs;
 }
 
+std::vector<double> obtainRBFGradientToTrainingPoints(
+    const std::vector<Vector> & parameter_points,
+    const std::string & interp_method, const std::string & rbf, double epsilon,
+    const Vector & point, const int index)
+{
+    std::vector<double> rbfs;
+    if (interp_method == "LS")
+    {
+        for (int i = 0; i < parameter_points.size(); i++)
+        {
+            rbfs.push_back(obtainRBFGradient(rbf, epsilon, point, parameter_points[i],
+                                             index));
+        }
+    }
+    else
+    {
+        CAROM_ERROR("Interpolated gradients are only implemented for \"LS\"");
+    }
+    return rbfs;
+}
+
+double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1,
+                         const Vector & point2, const int index)
+{
+    Vector diff;
+    point1.minus(point2, diff);
+    const double eps_norm_squared = epsilon * epsilon * diff.norm2();
+    double res = 0.0;
+
+    // Gaussian RBF
+    if (rbf == "G")
+    {
+        //Derivative of Gaussian RBF, res = std::exp(-eps_norm_squared);
+        res = -2.0*epsilon*epsilon*diff.item(index)*std::exp(-eps_norm_squared);
+    }
+    // Inverse quadratic RBF
+    else if (rbf == "IQ")
+    {
+        //Derivative of Inverse Quadratic RBF, res = 1.0 / (1.0 + eps_norm_squared);
+        res = -2.0*epsilon*epsilon*diff.item(index)/((1.0+eps_norm_squared)*
+                (1.0+eps_norm_squared));
+
+    }
+    // Inverse multiquadric RBF
+    else if (rbf == "IMQ")
+    {
+        //Derivative of Inverse multiquadritic RBF, res = 1.0 / std::sqrt(1.0 + eps_norm_squared);
+        res = -epsilon*epsilon*diff.item(index)/(std::sqrt(1.0+eps_norm_squared)*
+                (1.0+eps_norm_squared));
+    }
+
+    return res;
+}
+
 double rbfWeightedSum(std::vector<double>& rbf)
 {
     double sum = 0.0;
@@ -146,7 +202,7 @@ double obtainRBF(std::string rbf, double epsilon, const Vector & point1,
 {
     Vector diff;
     point1.minus(point2, diff);
-    double eps_norm_squared = epsilon * epsilon * diff.norm2();
+    const double eps_norm_squared = epsilon * epsilon * diff.norm2();
     double res = 0.0;
 
     // Gaussian RBF

lib/algo/manifold_interp/Interpolator.h

@@ -56,7 +56,8 @@ class Interpolator
                  int ref_point,
                  std::string rbf,
                  std::string interp_method,
-                 double closest_rbf_val = 0.9);
+                 double closest_rbf_val = 0.9,
+                 bool compute_gradients = false);
 
     /**
      * @brief The rank of the process this object belongs to.
@@ -108,6 +109,18 @@ class Interpolator
      */
     std::unique_ptr<Matrix> d_lambda_T;
 
+    /**
+     * @brief Flag that determines if a gradient with respect to the
+     *        parameters should be computed.
+     */
+    bool d_compute_gradients;
+
+    /**
+     * @brief Gradient with respect to the parameters.  Only exists after
+     *          interpolate has been run
+     */
+    std::vector<std::shared_ptr<Vector>> d_interpolation_gradient;
+
 private:
 
     /**
@@ -148,6 +161,42 @@ std::vector<double> obtainRBFToTrainingPoints(
     const std::string & interp_method, const std::string & rbf, double epsilon,
     const Vector & point);
 
+/**
+ * @brief Compute the RBF gradient from the parameter points with the
+ *        unsampled parameter point.
+ *
+ * @param[in] parameter_points The parameter points.
+ * @param[in] interp_method  The interpolation method type
+ *                           ("LS" == linear solve,
+ *                           "IDW" == inverse distance weighting,
+ *                           "LP" == lagrangian polynomials)
+ * @param[in] rbf Which RBF to compute.
+ * @param[in] epsilon   The RBF parameter that determines the width of
+                         influence.
+ * @param[in] point The unsampled parameter point.
+ * @param[in] index The index of the parameter that we are
+                     differentiating against.
+ */
+std::vector<double> obtainRBFGradientToTrainingPoints(
+    const std::vector<Vector> & parameter_points,
+    const std::string & interp_method, const std::string & rbf, double epsilon,
+    const Vector & point, const int index);
+
+
+/**
+ * @brief Compute the RBF gradient between two points.
+ *
+ * @param[in] rbf Which RBF to compute.
+ * @param[in] epsilon   The RBF parameter that determines the width of
+                         influence.
+ * @param[in] point1 The first point.
+ * @param[in] point2 The second point.
+ * @param[in] index The index of the parameter that we are
+                     differentiating against.
+ */
+double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1,
+                         const Vector & point2, const int index);
+
 /**
  * @brief Compute the sum of the RBF weights.
  *

lib/algo/manifold_interp/MatrixInterpolator.cpp

@@ -45,13 +45,15 @@ MatrixInterpolator::MatrixInterpolator(const std::vector<Vector> &
                                        std::string matrix_type,
                                        std::string rbf,
                                        std::string interp_method,
-                                       double closest_rbf_val) :
+                                       double closest_rbf_val,
+                                       bool compute_gradients) :
     Interpolator(parameter_points,
                  rotation_matrices,
                  ref_point,
                  rbf,
                  interp_method,
-                 closest_rbf_val)
+                 closest_rbf_val,
+                 compute_gradients)
 {
     CAROM_VERIFY(reduced_matrices.size() == rotation_matrices.size());
     CAROM_VERIFY(matrix_type == "SPD" || matrix_type == "NS" || matrix_type == "R"
@@ -102,10 +104,18 @@ std::shared_ptr<Matrix> MatrixInterpolator::interpolate(const Vector & point,
     std::shared_ptr<Matrix> interpolated_matrix;
     if (d_matrix_type == "SPD")
     {
+        if (d_compute_gradients)
+        {
+            CAROM_ERROR("Gradients are only implemented for B or G");
+        }
         interpolated_matrix = interpolateSPDMatrix(point);
     }
     else if (d_matrix_type == "NS")
     {
+        if (d_compute_gradients)
+        {
+            CAROM_ERROR("Gradients are only implemented for B or G");
+        }
         interpolated_matrix = interpolateNonSingularMatrix(point);
     }
     else
@@ -503,9 +513,28 @@ std::shared_ptr<Matrix> MatrixInterpolator::interpolateMatrix(
     // Interpolate gammas to get gamma for new point
     std::shared_ptr<Matrix> interpolated_matrix(obtainLogInterpolatedMatrix(rbf));
 
+    if (d_compute_gradients)
+    {
+        if(d_interp_method == "LS")
+        {
+            for (int i = 0; i < point.dim(); i++)
+            {
+                std::vector<double> rbf = obtainRBFGradientToTrainingPoints(d_parameter_points,
+                                          d_interp_method,
+                                          d_rbf, d_epsilon, point, i);
+
+                std::shared_ptr<Matrix> gradient_matrix(obtainLogInterpolatedMatrix(rbf));
+                d_interpolation_gradient.push_back(gradient_matrix);
+            }
+        }
+        else
+        {
+            CAROM_ERROR("Interpolated gradients are only implemented for \"LS\"");
+        }
+    }
+
     // The exp mapping is X + the interpolated gamma
     *interpolated_matrix += *d_rotated_reduced_matrices[d_ref_point];
     return interpolated_matrix;
 }
-
-}
+}

lib/algo/manifold_interp/MatrixInterpolator.h

@@ -58,6 +58,7 @@ class MatrixInterpolator : public Interpolator
      *                              "LP" == lagrangian polynomials)
      * @param[in] closest_rbf_val   The RBF parameter determines the width of influence.
      *                              Set the RBF value of the nearest two parameter points to a value between 0.0 to 1.0
+     * @param[in] compute_gradients Choose whether or not to compute the gradient along with the interpolation.
      */
     MatrixInterpolator(const std::vector<Vector> & parameter_points,
                        const std::vector<std::shared_ptr<Matrix>> & rotation_matrices,
@@ -66,7 +67,8 @@ class MatrixInterpolator : public Interpolator
                        std::string matrix_type,
                        std::string rbf = "G",
                        std::string interp_method = "LS",
-                       double closest_rbf_val = 0.9);
+                       double closest_rbf_val = 0.9,
+                       bool compute_gradients = false);
 
     /**
      * @brief Obtain the interpolated reduced matrix of the unsampled parameter point.
@@ -77,6 +79,14 @@ class MatrixInterpolator : public Interpolator
     std::shared_ptr<Matrix> interpolate(const Vector & point,
                                         bool orthogonalize = false);
 
+    /**
+    * @brief Returns the interpolated matrix's gradient.
+    *
+    */
+    std::vector<std::shared_ptr<Matrix>> getGradient() {
+        return d_interpolation_gradient;
+    }
+
 private:
 
     /**
@@ -154,6 +164,12 @@ class MatrixInterpolator : public Interpolator
      * @brief The reduced matrix of the reference point to the half power.
      */
     std::shared_ptr<Matrix> d_x_half_power;
+
+    /**
+     * @brief Gradient with respect to the parameters.  Only exists after
+     *          interpolate has been run
+     */
+    std::vector<std::shared_ptr<Matrix>> d_interpolation_gradient;
 };
 
 }

lib/algo/manifold_interp/VectorInterpolator.cpp

@@ -45,13 +45,15 @@ VectorInterpolator::VectorInterpolator(const std::vector<Vector> &
                                        int ref_point,
                                        std::string rbf,
                                        std::string interp_method,
-                                       double closest_rbf_val) :
+                                       double closest_rbf_val,
+                                       bool compute_gradients) :
     Interpolator(parameter_points,
                  rotation_matrices,
                  ref_point,
                  rbf,
                  interp_method,
-                 closest_rbf_val)
+                 closest_rbf_val,
+                 compute_gradients)
 {
     CAROM_VERIFY(reduced_vectors.size() == rotation_matrices.size());
 
@@ -120,6 +122,24 @@ std::shared_ptr<Vector> VectorInterpolator::interpolate(const Vector & point)
     std::unique_ptr<Vector> log_interpolated_vector = obtainLogInterpolatedVector(
                 rbf);
 
+    if (d_compute_gradients)
+    {
+        if(d_interp_method == "LS")
+        {
+            for (int i = 0; i < point.dim(); ++i)
+            {
+                std::vector<double> rbf = obtainRBFGradientToTrainingPoints(d_parameter_points,
+                                          d_interp_method,d_rbf, d_epsilon, point, i);
+                std::shared_ptr<Vector> gradient_vector(obtainLogInterpolatedVector(rbf));
+                d_interpolation_gradient.push_back(gradient_vector);
+            }
+        }
+        else
+        {
+            CAROM_ERROR("Interpolated gradients are only implemented for \"LS\"");
+        }
+    }
+
     // The exp mapping is X + the interpolated gamma
     std::shared_ptr<Vector> interpolated_vector =
         d_rotated_reduced_vectors[d_ref_point]->plus(

lib/algo/manifold_interp/VectorInterpolator.h

@@ -53,14 +53,16 @@ class VectorInterpolator : public Interpolator
      *                              "LP" == lagrangian polynomials)
      * @param[in] closest_rbf_val   The RBF parameter determines the width of influence.
      *                              Set the RBF value of the nearest two parameter points to a value between 0.0 to 1.0
+     * @param[in] compute_gradients Choose whether or not to compute the gradient along with the interpolation.
      */
     VectorInterpolator(const std::vector<Vector> & parameter_points,
                        const std::vector<std::shared_ptr<Matrix>> & rotation_matrices,
                        const std::vector<std::shared_ptr<Vector>> & reduced_vectors,
                        int ref_point,
                        std::string rbf = "G",
                        std::string interp_method = "LS",
-                       double closest_rbf_val = 0.9);
+                       double closest_rbf_val = 0.9,
+                       bool compute_gradients = false);
 
     /**
      * @brief Obtain the interpolated reduced vector of the unsampled parameter point.
@@ -69,6 +71,14 @@ class VectorInterpolator : public Interpolator
      */
     std::shared_ptr<Vector> interpolate(const Vector & point);
 
+    /**
+    * @brief Returns the interpolated matrix's gradient.
+    *
+    */
+    std::vector<std::shared_ptr<Vector>> getGradient() {
+        return d_interpolation_gradient;
+    }
+
 private:
 
     /**
@@ -113,6 +123,12 @@ class VectorInterpolator : public Interpolator
      * @brief The reduced elements in tangential space.
      */
     std::vector<std::shared_ptr<Vector>> d_gammas;
+
+    /**
+     * @brief Gradient with respect to the parameters.  Only exists after
+     *          interpolate has been run
+     */
+    std::vector<std::shared_ptr<Vector>> d_interpolation_gradient;
 };
 
 std::unique_ptr<Vector> obtainInterpolatedVector(

unit_tests/CMakeLists.txt

@@ -47,7 +47,8 @@ set(unit_test_stems
     GreedyCustomSampler
     NNLS
     basis_conversion
-    PCHIPInterpolator)
+    PCHIPInterpolator
+    interpolation_gradients)
 
 foreach(stem IN LISTS unit_test_stems)
     add_executable(test_${stem} test_${stem}.cpp)