diff --git a/lib/algo/manifold_interp/Interpolator.cpp b/lib/algo/manifold_interp/Interpolator.cpp index e97fc38a7..aa98bda91 100644 --- a/lib/algo/manifold_interp/Interpolator.cpp +++ b/lib/algo/manifold_interp/Interpolator.cpp @@ -34,7 +34,8 @@ Interpolator::Interpolator(const std::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 & 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 obtainRBFToTrainingPoints( return rbfs; } +std::vector obtainRBFGradientToTrainingPoints( + const std::vector & parameter_points, + const std::string & interp_method, const std::string & rbf, double epsilon, + const Vector & point, const int index) +{ + std::vector 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& 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 diff --git a/lib/algo/manifold_interp/Interpolator.h b/lib/algo/manifold_interp/Interpolator.h index 45090ff7a..e5adf837d 100644 --- a/lib/algo/manifold_interp/Interpolator.h +++ b/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 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> d_interpolation_gradient; + private: /** @@ -148,6 +161,42 @@ std::vector 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 obtainRBFGradientToTrainingPoints( + const std::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. * diff --git a/lib/algo/manifold_interp/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index 119ea1cf7..37dd16e81 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -45,13 +45,15 @@ MatrixInterpolator::MatrixInterpolator(const std::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 MatrixInterpolator::interpolate(const Vector & point, std::shared_ptr 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 MatrixInterpolator::interpolateMatrix( // Interpolate gammas to get gamma for new point std::shared_ptr interpolated_matrix(obtainLogInterpolatedMatrix(rbf)); + if (d_compute_gradients) + { + if(d_interp_method == "LS") + { + for (int i = 0; i < point.dim(); i++) + { + std::vector rbf = obtainRBFGradientToTrainingPoints(d_parameter_points, + d_interp_method, + d_rbf, d_epsilon, point, i); + + std::shared_ptr 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; } - -} +} \ No newline at end of file diff --git a/lib/algo/manifold_interp/MatrixInterpolator.h b/lib/algo/manifold_interp/MatrixInterpolator.h index 67a966e5f..a983c60ee 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.h +++ b/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 & parameter_points, const std::vector> & 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 interpolate(const Vector & point, bool orthogonalize = false); + /** + * @brief Returns the interpolated matrix's gradient. + * + */ + std::vector> 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 d_x_half_power; + + /** + * @brief Gradient with respect to the parameters. Only exists after + * interpolate has been run + */ + std::vector> d_interpolation_gradient; }; } diff --git a/lib/algo/manifold_interp/VectorInterpolator.cpp b/lib/algo/manifold_interp/VectorInterpolator.cpp index 128fcf196..516e5aafb 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.cpp +++ b/lib/algo/manifold_interp/VectorInterpolator.cpp @@ -45,13 +45,15 @@ VectorInterpolator::VectorInterpolator(const std::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 VectorInterpolator::interpolate(const Vector & point) std::unique_ptr log_interpolated_vector = obtainLogInterpolatedVector( rbf); + if (d_compute_gradients) + { + if(d_interp_method == "LS") + { + for (int i = 0; i < point.dim(); ++i) + { + std::vector rbf = obtainRBFGradientToTrainingPoints(d_parameter_points, + d_interp_method,d_rbf, d_epsilon, point, i); + std::shared_ptr 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 interpolated_vector = d_rotated_reduced_vectors[d_ref_point]->plus( diff --git a/lib/algo/manifold_interp/VectorInterpolator.h b/lib/algo/manifold_interp/VectorInterpolator.h index 930ffd3ba..b00ad3d50 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.h +++ b/lib/algo/manifold_interp/VectorInterpolator.h @@ -53,6 +53,7 @@ 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 & parameter_points, const std::vector> & rotation_matrices, @@ -60,7 +61,8 @@ class VectorInterpolator : public Interpolator 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 interpolate(const Vector & point); + /** + * @brief Returns the interpolated matrix's gradient. + * + */ + std::vector> getGradient() { + return d_interpolation_gradient; + } + private: /** @@ -113,6 +123,12 @@ class VectorInterpolator : public Interpolator * @brief The reduced elements in tangential space. */ std::vector> d_gammas; + + /** + * @brief Gradient with respect to the parameters. Only exists after + * interpolate has been run + */ + std::vector> d_interpolation_gradient; }; std::unique_ptr obtainInterpolatedVector( diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index 3532b126b..0410447c0 100644 --- a/unit_tests/CMakeLists.txt +++ b/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) diff --git a/unit_tests/test_interpolation_gradients.cpp b/unit_tests/test_interpolation_gradients.cpp new file mode 100644 index 000000000..3e14ac84a --- /dev/null +++ b/unit_tests/test_interpolation_gradients.cpp @@ -0,0 +1,517 @@ +/****************************************************************************** + * + * Copyright (c) 2013-2025, 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/MatrixInterpolator.h" +#include "algo/manifold_interp/VectorInterpolator.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include + +using namespace std; +using namespace CAROM; + +double func00(double* input) +{ + return 4*input[0] + 3*input[1]; +} +std::vector dfunc00(double* input) +{ + std::vector output{4.0, 3.0}; + return output; +} +double func01(double* input) +{ + return 2*input[0]*input[1] + 2*input[0]*input[0]; +} +std::vector dfunc01(double* input) +{ + std::vector output{2*input[1]+4*input[0], 2*input[0]}; + return output; +} +double func10(double* input) +{ + return input[1]*input[1] + 2.0; +} +std::vector dfunc10(double* input) +{ + std::vector output{0.0, 2*input[1]}; + return output; +} +double func11(double* input) +{ + return std::cos(3*input[0]); +} +std::vector dfunc11(double* input) +{ + std::vector output{-3*std::sin(3*input[0]), 0.0}; + return output; +} +std::shared_ptr createMatrix(CAROM::Vector input) +{ + std::shared_ptr output(new Matrix(2,2,false)); + output->item(0,0) = func00(input.getData()); + output->item(0,1) = func01(input.getData()); + output->item(1,0) = func10(input.getData()); + output->item(1,1) = func11(input.getData()); + return output; +} + +std::shared_ptr createVector(CAROM::Vector input) +{ + std::shared_ptr output(new Vector(2,false)); + output->item(0) = func00(input.getData()); + output->item(1) = func01(input.getData()); + return output; +} + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +TEST(InterpolationTest,GaussianMatrix) +{ + bool SUCCESS = true; + CAROM::Vector input1(2, false); + input1(0) = 0.04; + input1(1) = 0.04; + CAROM::Vector input2(2, false); + input2(0) = 0.04; + input2(1) = 0.06; + CAROM::Vector input3(2, false); + input3(0) = 0.06; + input3(1) = 0.04; + CAROM::Vector input4(2, false); + input4(0) = 0.06; + input4(1) = 0.06; + + CAROM::Vector target(2,false); + target(0) = 0.05; + target(1) = 0.05; + + //Interpolate a 2x2 matrix. We should have 4 different functions. + + std::shared_ptr matrix1 = createMatrix(input1); + std::shared_ptr matrix2 = createMatrix(input2); + std::shared_ptr matrix3 = createMatrix(input3); + std::shared_ptr matrix4 = createMatrix(input4); + + std::shared_ptr soln = createMatrix(target); + + std::shared_ptr eye(new Matrix(2,2,false)); + eye->item(0,0) = 1.0; + eye->item(0,1) = 0.0; + eye->item(0,1) = 0.0; + eye->item(1,1) = 1.0; + + std::vector inputs{input1, input2, input3, input4}; + std::vector> rotations{eye,eye,eye,eye}; + std::vector> matrices{matrix1, matrix2, matrix3, matrix4}; + + CAROM::MatrixInterpolator interpolator(inputs, rotations, matrices, 0, "B", "G", + "LS", 0.9, true); + + std::shared_ptr A = interpolator.interpolate(target); + std::vector> gradient = + interpolator.getGradient(); + + double epsilon = 1.0e-7; + CAROM::Vector perturb0(2, false); + perturb0(0) = epsilon; + perturb0(1) = 0.0; + perturb0 += target; + CAROM::Vector perturb1(2, false); + perturb1(0) = 0.0; + perturb1(1) = epsilon; + perturb1 += target; + + std::shared_ptr perturbation0 = createMatrix(perturb0); + std::shared_ptr perturbation1 = createMatrix(perturb1); + + std::shared_ptr interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); + + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); + *fd_grad_interp0 -= *A; + + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); + *fd_grad_interp1 -= *A; + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + fd_grad_interp0->item(i,j) /= epsilon; + fd_grad_interp1->item(i,j) /= epsilon; + } + } + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + double fd = fd_grad_interp0->item(i,j); + double my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + std::cout << "fd = " << fd << ", my = " << my << "abs(fd-my) = " << abs( + fd-my) << std::endl; + break; + } + fd = fd_grad_interp1->item(i,j); + my = gradient[1]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + std::cout << "fd = " << fd << ", my = " << my << "abs(fd-my) = " << abs( + fd-my) << std::endl; + break; + } + } + if(!SUCCESS) break; + } + EXPECT_TRUE(SUCCESS); +} +TEST(InterpolationTest,InverseQuadraticMatrix) +{ + bool SUCCESS = true; + CAROM::Vector input1(2, false); + input1(0) = 0.04; + input1(1) = 0.04; + CAROM::Vector input2(2, false); + input2(0) = 0.04; + input2(1) = 0.06; + CAROM::Vector input3(2, false); + input3(0) = 0.06; + input3(1) = 0.04; + CAROM::Vector input4(2, false); + input4(0) = 0.06; + input4(1) = 0.06; + + CAROM::Vector target(2,false); + target(0) = 0.05; + target(1) = 0.05; + + //Interpolate a 2x2 matrix. We should have 4 different functions. + + std::shared_ptr matrix1 = createMatrix(input1); + std::shared_ptr matrix2 = createMatrix(input2); + std::shared_ptr matrix3 = createMatrix(input3); + std::shared_ptr matrix4 = createMatrix(input4); + + std::shared_ptr soln = createMatrix(target); + + std::shared_ptr eye(new Matrix(2,2,false)); + eye->item(0,0) = 1.0; + eye->item(0,1) = 0.0; + eye->item(0,1) = 0.0; + eye->item(1,1) = 1.0; + + std::vector inputs{input1, input2, input3, input4}; + std::vector> rotations{eye,eye,eye,eye}; + std::vector> matrices{matrix1, matrix2, matrix3, matrix4}; + + CAROM::MatrixInterpolator interpolator(inputs, rotations, matrices, 0, "B", + "IQ", "LS", 0.9, true); + + std::shared_ptr A = interpolator.interpolate(target); + std::vector> gradient = + interpolator.getGradient(); + + double epsilon = 1.0e-7; + CAROM::Vector perturb0(2, false); + perturb0(0) = epsilon; + perturb0(1) = 0.0; + perturb0 += target; + CAROM::Vector perturb1(2, false); + perturb1(0) = 0.0; + perturb1(1) = epsilon; + perturb1 += target; + + std::shared_ptr perturbation0 = createMatrix(perturb0); + std::shared_ptr perturbation1 = createMatrix(perturb1); + + std::shared_ptr interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); + + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); + *fd_grad_interp0 -= *A; + + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); + *fd_grad_interp1 -= *A; + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + fd_grad_interp0->item(i,j) /= epsilon; + fd_grad_interp1->item(i,j) /= epsilon; + } + } + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + double fd = fd_grad_interp0->item(i,j); + double my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + break; + } + fd = fd_grad_interp0->item(i,j); + my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + break; + } + } + if(!SUCCESS) break; + } + EXPECT_TRUE(SUCCESS); +} + +TEST(InterpolationTest,InverseMultiQuadraticMatrix) +{ + bool SUCCESS = true; + CAROM::Vector input1(2, false); + input1(0) = 0.04; + input1(1) = 0.04; + CAROM::Vector input2(2, false); + input2(0) = 0.04; + input2(1) = 0.06; + CAROM::Vector input3(2, false); + input3(0) = 0.06; + input3(1) = 0.04; + CAROM::Vector input4(2, false); + input4(0) = 0.06; + input4(1) = 0.06; + + CAROM::Vector target(2,false); + target(0) = 0.05; + target(1) = 0.05; + + //Interpolate a 2x2 matrix. We should have 4 different functions. + + std::shared_ptr matrix1 = createMatrix(input1); + std::shared_ptr matrix2 = createMatrix(input2); + std::shared_ptr matrix3 = createMatrix(input3); + std::shared_ptr matrix4 = createMatrix(input4); + + std::shared_ptr soln = createMatrix(target); + + std::shared_ptr eye(new Matrix(2,2,false)); + eye->item(0,0) = 1.0; + eye->item(0,1) = 0.0; + eye->item(0,1) = 0.0; + eye->item(1,1) = 1.0; + + std::vector inputs{input1, input2, input3, input4}; + std::vector> rotations{eye,eye,eye,eye}; + std::vector> matrices{matrix1, matrix2, matrix3, matrix4}; + + CAROM::MatrixInterpolator interpolator(inputs, rotations, matrices, 0, "B", + "IMQ", "LS", 0.9, true); + + std::shared_ptr A = interpolator.interpolate(target); + std::vector> gradient = + interpolator.getGradient(); + + double epsilon = 1.0e-7; + CAROM::Vector perturb0(2, false); + perturb0(0) = epsilon; + perturb0(1) = 0.0; + perturb0 += target; + CAROM::Vector perturb1(2, false); + perturb1(0) = 0.0; + perturb1(1) = epsilon; + perturb1 += target; + + std::shared_ptr perturbation0 = createMatrix(perturb0); + std::shared_ptr perturbation1 = createMatrix(perturb1); + + std::shared_ptr interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); + + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); + *fd_grad_interp0 -= *A; + + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); + *fd_grad_interp1 -= *A; + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + fd_grad_interp0->item(i,j) /= epsilon; + fd_grad_interp1->item(i,j) /= epsilon; + } + } + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + double fd = fd_grad_interp0->item(i,j); + double my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + break; + } + fd = fd_grad_interp0->item(i,j); + my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + break; + } + } + if(!SUCCESS) break; + } + EXPECT_TRUE(SUCCESS); +} + +TEST(InterpolationTest,Vector) +{ + bool SUCCESS = true; + CAROM::Vector input1(2, false); + input1(0) = 0.04; + input1(1) = 0.04; + CAROM::Vector input2(2, false); + input2(0) = 0.04; + input2(1) = 0.06; + CAROM::Vector input3(2, false); + input3(0) = 0.06; + input3(1) = 0.04; + CAROM::Vector input4(2, false); + input4(0) = 0.06; + input4(1) = 0.06; + + CAROM::Vector target(2,false); + target(0) = 0.05; + target(1) = 0.05; + + + std::shared_ptr matrix1 = createVector(input1); + std::shared_ptr matrix2 = createVector(input2); + std::shared_ptr matrix3 = createVector(input3); + std::shared_ptr matrix4 = createVector(input4); + + std::shared_ptr soln = createVector(target); + + std::shared_ptr eye(new Matrix(2,2,false)); + eye->item(0,0) = 1.0; + eye->item(0,1) = 0.0; + eye->item(0,1) = 0.0; + eye->item(1,1) = 1.0; + + std::vector inputs{input1, input2, input3, input4}; + std::vector> rotations{eye,eye,eye,eye}; + std::vector> vectors{matrix1, matrix2, matrix3, matrix4}; + + CAROM::VectorInterpolator interpolator(inputs, rotations, vectors, 0, "G", "LS", + 0.9, true); + + std::shared_ptr v = interpolator.interpolate(target); + std::vector> gradient = + interpolator.getGradient(); + + double epsilon = 1.0e-7; + CAROM::Vector perturb0(2, false); + perturb0(0) = epsilon; + perturb0(1) = 0.0; + perturb0 += target; + CAROM::Vector perturb1(2, false); + perturb1(0) = 0.0; + perturb1(1) = epsilon; + perturb1 += target; + + std::shared_ptr perturbation0 = createVector(perturb0); + std::shared_ptr perturbation1 = createVector(perturb1); + + std::shared_ptr interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); + + std::shared_ptr fd_grad_interp0(new Vector( + *interp_perturbation0)); + *fd_grad_interp0 -= *v; + + std::shared_ptr fd_grad_interp1(new Vector( + *interp_perturbation1)); + *fd_grad_interp1 -= *v; + + for(int i = 0; i < 2; ++i) + { + fd_grad_interp0->item(i) /= epsilon; + fd_grad_interp1->item(i) /= epsilon; + } + + for(int i = 0; i < 2; ++i) + { + double fd = fd_grad_interp0->item(i); + double my = gradient[0]->item(i); + if(abs(fd - my) > 1.e-5) + { + SUCCESS = false; + break; + } + fd = fd_grad_interp0->item(i); + my = gradient[0]->item(i); + if(abs(fd - my) > 1.e-5) + { + 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