From 6c6f432dc7574aad3dded05cbe52335c851ff89e Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 3 Jun 2025 21:32:55 -0700 Subject: [PATCH 01/11] Added initial version of matrix interpolator. --- .../manifold_interp/MatrixInterpolator.cpp | 91 ++++++++++++++++++- lib/algo/manifold_interp/MatrixInterpolator.h | 56 +++++++++++- 2 files changed, 143 insertions(+), 4 deletions(-) diff --git a/lib/algo/manifold_interp/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index 119ea1cf7..402f1bf05 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -45,7 +45,8 @@ 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, @@ -57,6 +58,7 @@ MatrixInterpolator::MatrixInterpolator(const std::vector & CAROM_VERIFY(matrix_type == "SPD" || matrix_type == "NS" || matrix_type == "R" || matrix_type == "B"); d_matrix_type = matrix_type; + d_compute_gradients = compute_gradients; // Rotate the reduced matrices for (int i = 0; i < reduced_matrices.size(); i++) @@ -102,10 +104,20 @@ std::shared_ptr MatrixInterpolator::interpolate(const Vector & point, std::shared_ptr interpolated_matrix; if (d_matrix_type == "SPD") { + if (d_compute_gradients) + { + std::cout << "Gradients are only implemented for B or G"; + CAROM_VERIFY(false); + } interpolated_matrix = interpolateSPDMatrix(point); } else if (d_matrix_type == "NS") { + if (d_compute_gradients) + { + std::cout << "Gradients are only implemented for B or G"; + CAROM_VERIFY(false); + } interpolated_matrix = interpolateNonSingularMatrix(point); } else @@ -467,6 +479,58 @@ std::shared_ptr MatrixInterpolator::interpolateNonSingularMatrix( return interpolated_matrix; } +double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, + const Vector & point2, const int index) +{ + Vector diff; + point1.minus(point2, diff); + double eps_norm_squared = epsilon * epsilon * diff.norm2(); + double res = 0.0; + + // Gaussian RBF + if (rbf == "G") + { + //res = std::exp(-eps_norm_squared); + res = -2*epsilon*epsilon*diff.item(index)*std::exp(-eps_norm_squared); + } + // Inverse quadratic RBF + else if (rbf == "IQ") + { + //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") + { + //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; +} + +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 + { + std::cout << "Interpolated gradients are only implemented for \"LS\" "; + CAROM_VERIFY(interp_method == "LS"); + } + return rbfs; +} + std::shared_ptr MatrixInterpolator::interpolateMatrix( const Vector & point) { @@ -503,9 +567,30 @@ 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 interpolated_matrix(obtainLogInterpolatedMatrix(rbf)); + + d_interpolation_gradient.push_back(interpolated_matrix); + } + } + else + { + std::cout << "Interpolated gradients are only implemented for \"LS\" "; + CAROM_VERIFY(d_interp_method == "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..b18e87a79 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.h +++ b/lib/algo/manifold_interp/MatrixInterpolator.h @@ -66,7 +66,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 +78,12 @@ 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: /** @@ -135,6 +142,41 @@ class MatrixInterpolator : public Interpolator */ std::shared_ptr interpolateMatrix(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 The matrix type/manifold */ @@ -154,6 +196,18 @@ class MatrixInterpolator : public Interpolator * @brief The reduced matrix of the reference point to the half power. */ std::shared_ptr d_x_half_power; + + /** + * @brief Flag that determines if a gradient with respect to the + * parameters should be computed. + */ + bool d_compute_gradients; + + /** + * @brief Flag that determines if a gradient with respect to the + * parameters should be computed. + */ + std::vector> d_interpolation_gradient; }; } From 330157e1dd553459014b9b49ad3f28ffe7f50657 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 3 Jun 2025 21:33:57 -0700 Subject: [PATCH 02/11] Added initial Matrix changes. --- lib/linalg/Matrix.cpp | 116 ++++++++++++++++++++++++++++++++++++++++++ lib/linalg/Matrix.h | 13 +++++ 2 files changed, 129 insertions(+) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index 559176bbc..c21dfa2ae 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -2103,6 +2103,122 @@ ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A) return eigenpair; } +std::vector NonSymmetricEigenSolve(const Matrix & A) +{ + char jobvl = 'V', jobrl = 'V'; + + int info; + int k = A.numColumns(); + int lwork = std::max(k*k, 10*k); + double* work = new double [lwork]; + double* e_real = new double [k]; + double* e_imaginary = new double [k]; + Matrix* ev_l = new Matrix(k, k, false); + Matrix* ev_r = new Matrix(k, k, false); + Matrix* A_copy = new Matrix(A); + + // A now in a row major representation. Put it + // into column major order. + for (int row = 0; row < k; ++row) { + for (int col = row+1; col < k; ++col) { + double tmp = A_copy->item(row, col); + A_copy->item(row, col) = A_copy->item(col, row); + A_copy->item(col, row) = tmp; + } + } + + // Now call lapack to do the eigensolve. + dgeev(&jobvl, &jobrl, &k, A_copy->getData(), &k, e_real, e_imaginary, ev_l->getData(), + &k, ev_r->getData(), &k, work, &lwork, &info); + + // Eigenvectors now in a column major representation. Put it + // into row major order. + for (int row = 0; row < k; ++row) { + for (int col = row+1; col < k; ++col) { + double tmp = ev_r->item(row, col); + ev_r->item(row, col) = ev_r->item(col, row); + ev_r->item(col, row) = tmp; + + tmp = ev_l->item(row, col); + ev_l->item(row, col) = ev_l->item(col,row); + ev_l->item(col,row) = tmp; + } + } + + std::vector eigenpairs; + + ComplexEigenPair eigenpair_right; + eigenpair_right.ev_real = new Matrix(k, k, false); + eigenpair_right.ev_imaginary = new Matrix(k, k, false); + + // Separate lapack eigenvector into real and imaginary parts + for (int i = 0; i < k; ++i) + { + for (int row = 0; row < k; ++row) { + eigenpair_right.ev_real->item(row, i) = ev_r->item(row, i); + } + if (e_imaginary[i] != 0) + { + for (int row = 0; row < k; ++row) { + eigenpair_right.ev_real->item(row, i + 1) = ev_r->item(row, i); + eigenpair_right.ev_imaginary->item(row, i) = ev_r->item(row, i + 1); + eigenpair_right.ev_imaginary->item(row, i + 1) = -ev_r->item(row, i + 1); + } + + // Skip the next eigenvalue since it'll be part of the complex + // conjugate pair. + ++i; + } + } + + for (int i = 0; i < k; i++) + { + eigenpair_right.eigs.push_back(std::complex(e_real[i], e_imaginary[i])); + } + + eigenpairs.push_back(eigenpair_right); + + ComplexEigenPair eigenpair_left; + eigenpair_left.ev_real = new Matrix(k, k, false); + eigenpair_left.ev_imaginary = new Matrix(k, k, false); + + // Separate lapack eigenvector into real and imaginary parts + for (int i = 0; i < k; ++i) + { + for (int row = 0; row < k; ++row) { + eigenpair_left.ev_real->item(row, i) = ev_l->item(row, i); + } + if (e_imaginary[i] != 0) + { + for (int row = 0; row < k; ++row) { + eigenpair_left.ev_real->item(row, i + 1) = ev_l->item(row, i); + eigenpair_left.ev_imaginary->item(row, i) = ev_l->item(row, i + 1); + eigenpair_left.ev_imaginary->item(row, i + 1) = -ev_l->item(row, i + 1); + } + + // Skip the next eigenvalue since it'll be part of the complex + // conjugate pair. + ++i; + } + } + + for (int i = 0; i < k; i++) + { + eigenpair_left.eigs.push_back(std::complex(e_real[i], e_imaginary[i])); + } + + eigenpairs.push_back(eigenpair_left); + + delete [] work; + delete [] e_real; + delete [] e_imaginary; + delete ev_r; + delete ev_l; + delete A_copy; + + return eigenpairs; +} + void SerialSVD(Matrix* A, Matrix* U, Vector* S, diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index d6e094f42..09ae4a5b7 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -1205,6 +1205,19 @@ EigenPair SymmetricRightEigenSolve(const Matrix & A); */ ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A); +/** + * @brief Computes the left and right eigenvectors/eigenvalues of an NxN + * real nonsymmetric matrix. + * + * @param[in] A The NxN real nonsymmetric matrix to be eigendecomposed. + * + * @return A standard vector of ComplexEigenPairs containing the + * right, then left, eigenvectors and eigenvalues of the eigensolve. + * The eigenvector matrices contained within the returning struct + * must be destroyed by the user. + */ +std::vector NonSymmetricEigenSolve(const Matrix & A); + std::unique_ptr SpaceTimeProduct(const CAROM::Matrix & As, const CAROM::Matrix & At, const CAROM::Matrix & Bs, const CAROM::Matrix & Bt, From 96fa182af56f91a93a4344e8431596a76552f64c Mon Sep 17 00:00:00 2001 From: ptranq Date: Wed, 4 Jun 2025 16:15:01 -0700 Subject: [PATCH 03/11] Added test and improved matrix interpolator gradients. --- .../manifold_interp/MatrixInterpolator.cpp | 7 +- unit_tests/CMakeLists.txt | 3 +- .../test_MatrixInterpolator_gradients.cpp | 201 ++++++++++++++++++ 3 files changed, 207 insertions(+), 4 deletions(-) create mode 100644 unit_tests/test_MatrixInterpolator_gradients.cpp diff --git a/lib/algo/manifold_interp/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index 402f1bf05..78a0e827c 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -479,7 +479,7 @@ std::shared_ptr MatrixInterpolator::interpolateNonSingularMatrix( return interpolated_matrix; } -double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, +double MatrixInterpolator::obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, const Vector & point2, const int index) { Vector diff; @@ -491,7 +491,7 @@ double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, if (rbf == "G") { //res = std::exp(-eps_norm_squared); - res = -2*epsilon*epsilon*diff.item(index)*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") @@ -510,7 +510,7 @@ double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, return res; } -std::vector obtainRBFGradientToTrainingPoints( +std::vector MatrixInterpolator::obtainRBFGradientToTrainingPoints( const std::vector & parameter_points, const std::string & interp_method, const std::string & rbf, double epsilon, const Vector & point, const int index) @@ -578,6 +578,7 @@ std::shared_ptr MatrixInterpolator::interpolateMatrix( d_rbf, d_epsilon, point, i); std::shared_ptr interpolated_matrix(obtainLogInterpolatedMatrix(rbf)); + // *interpolated_matrix += *d_rotated_reduced_matrices[d_ref_point]; d_interpolation_gradient.push_back(interpolated_matrix); } diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index 3532b126b..025e5621f 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -17,7 +17,8 @@ set(unit_test_names weak_scaling random_test smoke_static - load_samples) + load_samples + test_MatrixInterpolator_gradients) foreach(name IN LISTS unit_test_names) add_executable(${name} ${name}.cpp) diff --git a/unit_tests/test_MatrixInterpolator_gradients.cpp b/unit_tests/test_MatrixInterpolator_gradients.cpp new file mode 100644 index 000000000..5ce970e9d --- /dev/null +++ b/unit_tests/test_MatrixInterpolator_gradients.cpp @@ -0,0 +1,201 @@ +/****************************************************************************** + * + * 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) + * + *****************************************************************************/ + +#include "algo/manifold_interp/MatrixInterpolator.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include +#include +#include + +#ifndef _WIN32 +#include // mkdir +#else +#include // _mkdir +#define mkdir(dir, mode) _mkdir(dir) +#endif + +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::vector> createGradient(CAROM::Vector input) +{ + std::vector grad00 = dfunc00(input.getData()); + std::vector grad01 = dfunc01(input.getData()); + std::vector grad10 = dfunc10(input.getData()); + std::vector grad11 = dfunc11(input.getData()); + + std::shared_ptr gradient0(new Matrix(2,2,false)); + gradient0->item(0,0) = grad00[0]; gradient0->item(0,1) = grad01[0]; + gradient0->item(1,0) = grad10[0]; gradient0->item(1,1) = grad11[0]; + + std::shared_ptr gradient1(new Matrix(2,2,false)); + gradient1->item(0,0) = grad00[1]; gradient1->item(0,1) = grad01[1]; + gradient1->item(1,0) = grad10[1]; gradient1->item(1,1) = grad11[1]; + + std::vector> gradient{gradient0, gradient1}; + + return gradient; +} + +void printMatrix(std::shared_ptr input, int rows, int cols, string name) +{ + std::cout << name << "=" << std::endl; + for(int i = 0; i < rows; ++i) + { + for(int j = 0; j < cols; ++j) + { + std::cout << input->item(i,j) << ", "; + } + std::cout << std::endl; + } + std::cout << "-------------------" << std::endl; +} + +int main(int argc, char *argv[]) +{ + // CAROM::Vector input1(2, false); + // input1(0) = 0.0; input1(1) = 0.0; + // CAROM::Vector input2(2, false); + // input2(0) = 0.0; input2(1) = 0.1; + // CAROM::Vector input3(2, false); + // input3(0) = 0.1; input3(1) = 0.0; + // CAROM::Vector input4(2, false); + // input4(0) = 0.1; input4(1) = 0.1; + 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::vector> true_gradient = createGradient(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 = std::sqrt(std::numeric_limits::denorm_min()); + 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 fd_grad0(new Matrix(perturbation0->getData(), 2,2,false)); + *fd_grad0 -= *soln; + + std::shared_ptr fd_grad1(new Matrix(perturbation1->getData(), 2,2,false)); + *fd_grad1 -= *soln; + + for(int i = 0; i < 2; ++i) + { + for(int j = 0; j < 2; ++j) + { + fd_grad0->item(i,j) *= 1./epsilon; + fd_grad1->item(i,j) *= 1./epsilon; + } + } + + + + + + + + + printMatrix(soln, 2,2, "true solution"); + printMatrix(A,2,2,"interpolated Matrix"); + + printMatrix(true_gradient[0],2,2,"true_grad0"); + printMatrix(gradient[0],2,2,"grad0"); + printMatrix(fd_grad0, 2,2, "finite difference grad0"); + + printMatrix(true_gradient[1],2,2,"true_grad1"); + printMatrix(gradient[1],2,2,"grad1"); + printMatrix(fd_grad1, 2,2, "finite difference grad1"); + +} \ No newline at end of file From f0519b8cbc5461186f83e2777951165dc63abef6 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 19:45:12 -0700 Subject: [PATCH 04/11] Swapped the interpolation gradients to the Interpolator class. --- lib/algo/manifold_interp/Interpolator.cpp | 58 +++++++++++++++++- lib/algo/manifold_interp/Interpolator.h | 53 +++++++++++++++- .../manifold_interp/MatrixInterpolator.cpp | 60 ++----------------- lib/algo/manifold_interp/MatrixInterpolator.h | 45 +------------- .../test_MatrixInterpolator_gradients.cpp | 25 +++++--- 5 files changed, 130 insertions(+), 111 deletions(-) diff --git a/lib/algo/manifold_interp/Interpolator.cpp b/lib/algo/manifold_interp/Interpolator.cpp index e97fc38a7..69694a56a 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 + { + std::cout << "Interpolated gradients are only implemented for \"LS\" "; + CAROM_VERIFY(interp_method == "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); + double eps_norm_squared = epsilon * epsilon * diff.norm2(); + double res = 0.0; + + // Gaussian RBF + if (rbf == "G") + { + //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") + { + //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") + { + //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; diff --git a/lib/algo/manifold_interp/Interpolator.h b/lib/algo/manifold_interp/Interpolator.h index 45090ff7a..24a8ad8d5 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: /** @@ -144,9 +157,45 @@ class Interpolator * @param[in] point The unsampled parameter point. */ std::vector obtainRBFToTrainingPoints( +const std::vector & parameter_points, +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 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 78a0e827c..01478c04e 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -52,13 +52,13 @@ MatrixInterpolator::MatrixInterpolator(const std::vector & 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" || matrix_type == "B"); d_matrix_type = matrix_type; - d_compute_gradients = compute_gradients; // Rotate the reduced matrices for (int i = 0; i < reduced_matrices.size(); i++) @@ -479,58 +479,6 @@ std::shared_ptr MatrixInterpolator::interpolateNonSingularMatrix( return interpolated_matrix; } -double MatrixInterpolator::obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, - const Vector & point2, const int index) -{ - Vector diff; - point1.minus(point2, diff); - double eps_norm_squared = epsilon * epsilon * diff.norm2(); - double res = 0.0; - - // Gaussian RBF - if (rbf == "G") - { - //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") - { - //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") - { - //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; -} - -std::vector MatrixInterpolator::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 - { - std::cout << "Interpolated gradients are only implemented for \"LS\" "; - CAROM_VERIFY(interp_method == "LS"); - } - return rbfs; -} - std::shared_ptr MatrixInterpolator::interpolateMatrix( const Vector & point) { @@ -577,10 +525,10 @@ std::shared_ptr MatrixInterpolator::interpolateMatrix( d_interp_method, d_rbf, d_epsilon, point, i); - std::shared_ptr interpolated_matrix(obtainLogInterpolatedMatrix(rbf)); + std::shared_ptr gradient_matrix(obtainLogInterpolatedMatrix(rbf)); // *interpolated_matrix += *d_rotated_reduced_matrices[d_ref_point]; - d_interpolation_gradient.push_back(interpolated_matrix); + d_interpolation_gradient.push_back(gradient_matrix); } } else diff --git a/lib/algo/manifold_interp/MatrixInterpolator.h b/lib/algo/manifold_interp/MatrixInterpolator.h index b18e87a79..9fe98c5b3 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.h +++ b/lib/algo/manifold_interp/MatrixInterpolator.h @@ -142,41 +142,6 @@ class MatrixInterpolator : public Interpolator */ std::shared_ptr interpolateMatrix(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 The matrix type/manifold */ @@ -196,16 +161,10 @@ class MatrixInterpolator : public Interpolator * @brief The reduced matrix of the reference point to the half power. */ std::shared_ptr d_x_half_power; - - /** - * @brief Flag that determines if a gradient with respect to the - * parameters should be computed. - */ - bool d_compute_gradients; /** - * @brief Flag that determines if a gradient with respect to the - * parameters should be computed. + * @brief Gradient with respect to the parameters. Only exists after + * interpolate has been run */ std::vector> d_interpolation_gradient; }; diff --git a/unit_tests/test_MatrixInterpolator_gradients.cpp b/unit_tests/test_MatrixInterpolator_gradients.cpp index 5ce970e9d..5f72d5231 100644 --- a/unit_tests/test_MatrixInterpolator_gradients.cpp +++ b/unit_tests/test_MatrixInterpolator_gradients.cpp @@ -165,37 +165,44 @@ int main(int argc, char *argv[]) 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_grad0(new Matrix(perturbation0->getData(), 2,2,false)); *fd_grad0 -= *soln; std::shared_ptr fd_grad1(new Matrix(perturbation1->getData(), 2,2,false)); *fd_grad1 -= *soln; + 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_grad0->item(i,j) *= 1./epsilon; fd_grad1->item(i,j) *= 1./epsilon; + fd_grad_interp0->item(i,j) /= epsilon; + fd_grad_interp1->item(i,j) /= epsilon; } } - - - - - - - printMatrix(soln, 2,2, "true solution"); printMatrix(A,2,2,"interpolated Matrix"); printMatrix(true_gradient[0],2,2,"true_grad0"); - printMatrix(gradient[0],2,2,"grad0"); printMatrix(fd_grad0, 2,2, "finite difference grad0"); + printMatrix(gradient[0],2,2,"grad0"); + printMatrix(fd_grad_interp0, 2,2, "finite difference interpolation grad0"); + printMatrix(true_gradient[1],2,2,"true_grad1"); - printMatrix(gradient[1],2,2,"grad1"); printMatrix(fd_grad1, 2,2, "finite difference grad1"); + printMatrix(gradient[1],2,2,"grad1"); + printMatrix(fd_grad_interp1, 2,2, "finite difference interpolation grad0"); } \ No newline at end of file From f3ee8d4250b53f0d5b28dcb517b8975a39a255c3 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 20:56:33 -0700 Subject: [PATCH 05/11] Stylized and added new gtest unit test. --- lib/algo/manifold_interp/MatrixInterpolator.h | 1 + .../manifold_interp/VectorInterpolator.cpp | 25 +- lib/algo/manifold_interp/VectorInterpolator.h | 16 +- unit_tests/CMakeLists.txt | 3 +- .../test_MatrixInterpolator_gradients.cpp | 52 +- unit_tests/test_interpolation_gradients.cpp | 455 ++++++++++++++++++ 6 files changed, 521 insertions(+), 31 deletions(-) create mode 100644 unit_tests/test_interpolation_gradients.cpp diff --git a/lib/algo/manifold_interp/MatrixInterpolator.h b/lib/algo/manifold_interp/MatrixInterpolator.h index 9fe98c5b3..5717377c5 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, diff --git a/lib/algo/manifold_interp/VectorInterpolator.cpp b/lib/algo/manifold_interp/VectorInterpolator.cpp index 128fcf196..3cdaae1b2 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,25 @@ 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 + { + std::cout << "Interpolated gradients are only implemented for \"LS\" "; + CAROM_VERIFY(d_interp_method == "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..2eab3b291 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,12 @@ 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 +121,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 025e5621f..fcd174fdb 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -48,7 +48,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_MatrixInterpolator_gradients.cpp b/unit_tests/test_MatrixInterpolator_gradients.cpp index 5f72d5231..d2f0f4c2f 100644 --- a/unit_tests/test_MatrixInterpolator_gradients.cpp +++ b/unit_tests/test_MatrixInterpolator_gradients.cpp @@ -106,14 +106,7 @@ void printMatrix(std::shared_ptr input, int rows, int cols, strin int main(int argc, char *argv[]) { - // CAROM::Vector input1(2, false); - // input1(0) = 0.0; input1(1) = 0.0; - // CAROM::Vector input2(2, false); - // input2(0) = 0.0; input2(1) = 0.1; - // CAROM::Vector input3(2, false); - // input3(0) = 0.1; input3(1) = 0.0; - // CAROM::Vector input4(2, false); - // input4(0) = 0.1; input4(1) = 0.1; + bool SUCCESS = true; CAROM::Vector input1(2, false); input1(0) = 0.04; input1(1) = 0.04; CAROM::Vector input2(2, false); @@ -135,9 +128,6 @@ int main(int argc, char *argv[]) std::shared_ptr soln = createMatrix(target); - - std::vector> true_gradient = createGradient(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; @@ -146,14 +136,11 @@ int main(int argc, char *argv[]) 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); + 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 = std::sqrt(std::numeric_limits::denorm_min()); double epsilon = 1.0e-7; CAROM::Vector perturb0(2, false); perturb0(0) = epsilon; perturb0(1) = 0.0; @@ -168,12 +155,6 @@ int main(int argc, char *argv[]) std::shared_ptr interp_perturbation0 = interpolator.interpolate(perturb0); std::shared_ptr interp_perturbation1 = interpolator.interpolate(perturb1); - std::shared_ptr fd_grad0(new Matrix(perturbation0->getData(), 2,2,false)); - *fd_grad0 -= *soln; - - std::shared_ptr fd_grad1(new Matrix(perturbation1->getData(), 2,2,false)); - *fd_grad1 -= *soln; - std::shared_ptr fd_grad_interp0(new Matrix(*interp_perturbation0)); *fd_grad_interp0 -= *A; @@ -184,24 +165,41 @@ int main(int argc, char *argv[]) { for(int j = 0; j < 2; ++j) { - fd_grad0->item(i,j) *= 1./epsilon; - fd_grad1->item(i,j) *= 1./epsilon; 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-7) + { + SUCCESS = false; + break; + } + fd = fd_grad_interp0->item(i,j); + my = gradient[0]->item(i,j); + if(abs(fd - my) > 1.e-7) + { + SUCCESS = false; + break; + } + } + if(!SUCCESS) break; + } + // EXPECT_TRUE(SUCCESS); + printMatrix(soln, 2,2, "true solution"); printMatrix(A,2,2,"interpolated Matrix"); - printMatrix(true_gradient[0],2,2,"true_grad0"); - printMatrix(fd_grad0, 2,2, "finite difference grad0"); printMatrix(gradient[0],2,2,"grad0"); printMatrix(fd_grad_interp0, 2,2, "finite difference interpolation grad0"); - printMatrix(true_gradient[1],2,2,"true_grad1"); - printMatrix(fd_grad1, 2,2, "finite difference grad1"); printMatrix(gradient[1],2,2,"grad1"); printMatrix(fd_grad_interp1, 2,2, "finite difference interpolation grad0"); diff --git a/unit_tests/test_interpolation_gradients.cpp b/unit_tests/test_interpolation_gradients.cpp new file mode 100644 index 000000000..eaa0d89c0 --- /dev/null +++ b/unit_tests/test_interpolation_gradients.cpp @@ -0,0 +1,455 @@ +/****************************************************************************** + * + * 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/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 From f89a78520bb47357eae0a1bd7b544d3275cc28e5 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 20:56:55 -0700 Subject: [PATCH 06/11] Removed old test from cmakelists. --- unit_tests/CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index fcd174fdb..0410447c0 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -17,8 +17,7 @@ set(unit_test_names weak_scaling random_test smoke_static - load_samples - test_MatrixInterpolator_gradients) + load_samples) foreach(name IN LISTS unit_test_names) add_executable(${name} ${name}.cpp) From 41cc59abed16fca38baa1fa20d7819157fb1e97d Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 20:57:21 -0700 Subject: [PATCH 07/11] Removed old unit test --- .../test_MatrixInterpolator_gradients.cpp | 206 ------------------ 1 file changed, 206 deletions(-) delete mode 100644 unit_tests/test_MatrixInterpolator_gradients.cpp diff --git a/unit_tests/test_MatrixInterpolator_gradients.cpp b/unit_tests/test_MatrixInterpolator_gradients.cpp deleted file mode 100644 index d2f0f4c2f..000000000 --- a/unit_tests/test_MatrixInterpolator_gradients.cpp +++ /dev/null @@ -1,206 +0,0 @@ -/****************************************************************************** - * - * 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) - * - *****************************************************************************/ - -#include "algo/manifold_interp/MatrixInterpolator.h" -#include "linalg/Matrix.h" -#include "linalg/Vector.h" -#include -#include -#include - -#ifndef _WIN32 -#include // mkdir -#else -#include // _mkdir -#define mkdir(dir, mode) _mkdir(dir) -#endif - -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::vector> createGradient(CAROM::Vector input) -{ - std::vector grad00 = dfunc00(input.getData()); - std::vector grad01 = dfunc01(input.getData()); - std::vector grad10 = dfunc10(input.getData()); - std::vector grad11 = dfunc11(input.getData()); - - std::shared_ptr gradient0(new Matrix(2,2,false)); - gradient0->item(0,0) = grad00[0]; gradient0->item(0,1) = grad01[0]; - gradient0->item(1,0) = grad10[0]; gradient0->item(1,1) = grad11[0]; - - std::shared_ptr gradient1(new Matrix(2,2,false)); - gradient1->item(0,0) = grad00[1]; gradient1->item(0,1) = grad01[1]; - gradient1->item(1,0) = grad10[1]; gradient1->item(1,1) = grad11[1]; - - std::vector> gradient{gradient0, gradient1}; - - return gradient; -} - -void printMatrix(std::shared_ptr input, int rows, int cols, string name) -{ - std::cout << name << "=" << std::endl; - for(int i = 0; i < rows; ++i) - { - for(int j = 0; j < cols; ++j) - { - std::cout << input->item(i,j) << ", "; - } - std::cout << std::endl; - } - std::cout << "-------------------" << std::endl; -} - -int main(int argc, char *argv[]) -{ - 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-7) - { - SUCCESS = false; - break; - } - fd = fd_grad_interp0->item(i,j); - my = gradient[0]->item(i,j); - if(abs(fd - my) > 1.e-7) - { - SUCCESS = false; - break; - } - } - if(!SUCCESS) break; - } - // EXPECT_TRUE(SUCCESS); - - printMatrix(soln, 2,2, "true solution"); - printMatrix(A,2,2,"interpolated Matrix"); - - printMatrix(gradient[0],2,2,"grad0"); - printMatrix(fd_grad_interp0, 2,2, "finite difference interpolation grad0"); - - - printMatrix(gradient[1],2,2,"grad1"); - printMatrix(fd_grad_interp1, 2,2, "finite difference interpolation grad0"); - -} \ No newline at end of file From 087ee30291a8695dde3f4639a8cca7a2779d9ff2 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 21:07:14 -0700 Subject: [PATCH 08/11] Reverted changes to matrix class. This will come back later in DMD gradients. It is too soon for this now. --- lib/linalg/Matrix.cpp | 116 ------------------------------------------ lib/linalg/Matrix.h | 13 ----- 2 files changed, 129 deletions(-) diff --git a/lib/linalg/Matrix.cpp b/lib/linalg/Matrix.cpp index c21dfa2ae..559176bbc 100644 --- a/lib/linalg/Matrix.cpp +++ b/lib/linalg/Matrix.cpp @@ -2103,122 +2103,6 @@ ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A) return eigenpair; } -std::vector NonSymmetricEigenSolve(const Matrix & A) -{ - char jobvl = 'V', jobrl = 'V'; - - int info; - int k = A.numColumns(); - int lwork = std::max(k*k, 10*k); - double* work = new double [lwork]; - double* e_real = new double [k]; - double* e_imaginary = new double [k]; - Matrix* ev_l = new Matrix(k, k, false); - Matrix* ev_r = new Matrix(k, k, false); - Matrix* A_copy = new Matrix(A); - - // A now in a row major representation. Put it - // into column major order. - for (int row = 0; row < k; ++row) { - for (int col = row+1; col < k; ++col) { - double tmp = A_copy->item(row, col); - A_copy->item(row, col) = A_copy->item(col, row); - A_copy->item(col, row) = tmp; - } - } - - // Now call lapack to do the eigensolve. - dgeev(&jobvl, &jobrl, &k, A_copy->getData(), &k, e_real, e_imaginary, ev_l->getData(), - &k, ev_r->getData(), &k, work, &lwork, &info); - - // Eigenvectors now in a column major representation. Put it - // into row major order. - for (int row = 0; row < k; ++row) { - for (int col = row+1; col < k; ++col) { - double tmp = ev_r->item(row, col); - ev_r->item(row, col) = ev_r->item(col, row); - ev_r->item(col, row) = tmp; - - tmp = ev_l->item(row, col); - ev_l->item(row, col) = ev_l->item(col,row); - ev_l->item(col,row) = tmp; - } - } - - std::vector eigenpairs; - - ComplexEigenPair eigenpair_right; - eigenpair_right.ev_real = new Matrix(k, k, false); - eigenpair_right.ev_imaginary = new Matrix(k, k, false); - - // Separate lapack eigenvector into real and imaginary parts - for (int i = 0; i < k; ++i) - { - for (int row = 0; row < k; ++row) { - eigenpair_right.ev_real->item(row, i) = ev_r->item(row, i); - } - if (e_imaginary[i] != 0) - { - for (int row = 0; row < k; ++row) { - eigenpair_right.ev_real->item(row, i + 1) = ev_r->item(row, i); - eigenpair_right.ev_imaginary->item(row, i) = ev_r->item(row, i + 1); - eigenpair_right.ev_imaginary->item(row, i + 1) = -ev_r->item(row, i + 1); - } - - // Skip the next eigenvalue since it'll be part of the complex - // conjugate pair. - ++i; - } - } - - for (int i = 0; i < k; i++) - { - eigenpair_right.eigs.push_back(std::complex(e_real[i], e_imaginary[i])); - } - - eigenpairs.push_back(eigenpair_right); - - ComplexEigenPair eigenpair_left; - eigenpair_left.ev_real = new Matrix(k, k, false); - eigenpair_left.ev_imaginary = new Matrix(k, k, false); - - // Separate lapack eigenvector into real and imaginary parts - for (int i = 0; i < k; ++i) - { - for (int row = 0; row < k; ++row) { - eigenpair_left.ev_real->item(row, i) = ev_l->item(row, i); - } - if (e_imaginary[i] != 0) - { - for (int row = 0; row < k; ++row) { - eigenpair_left.ev_real->item(row, i + 1) = ev_l->item(row, i); - eigenpair_left.ev_imaginary->item(row, i) = ev_l->item(row, i + 1); - eigenpair_left.ev_imaginary->item(row, i + 1) = -ev_l->item(row, i + 1); - } - - // Skip the next eigenvalue since it'll be part of the complex - // conjugate pair. - ++i; - } - } - - for (int i = 0; i < k; i++) - { - eigenpair_left.eigs.push_back(std::complex(e_real[i], e_imaginary[i])); - } - - eigenpairs.push_back(eigenpair_left); - - delete [] work; - delete [] e_real; - delete [] e_imaginary; - delete ev_r; - delete ev_l; - delete A_copy; - - return eigenpairs; -} - void SerialSVD(Matrix* A, Matrix* U, Vector* S, diff --git a/lib/linalg/Matrix.h b/lib/linalg/Matrix.h index 09ae4a5b7..d6e094f42 100644 --- a/lib/linalg/Matrix.h +++ b/lib/linalg/Matrix.h @@ -1205,19 +1205,6 @@ EigenPair SymmetricRightEigenSolve(const Matrix & A); */ ComplexEigenPair NonSymmetricRightEigenSolve(const Matrix & A); -/** - * @brief Computes the left and right eigenvectors/eigenvalues of an NxN - * real nonsymmetric matrix. - * - * @param[in] A The NxN real nonsymmetric matrix to be eigendecomposed. - * - * @return A standard vector of ComplexEigenPairs containing the - * right, then left, eigenvectors and eigenvalues of the eigensolve. - * The eigenvector matrices contained within the returning struct - * must be destroyed by the user. - */ -std::vector NonSymmetricEigenSolve(const Matrix & A); - std::unique_ptr SpaceTimeProduct(const CAROM::Matrix & As, const CAROM::Matrix & At, const CAROM::Matrix & Bs, const CAROM::Matrix & Bt, From e9bb1852151fa79e8d7fd552bd05a172e3c08f3e Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 21:31:06 -0700 Subject: [PATCH 09/11] Ran stylize again. --- lib/algo/manifold_interp/Interpolator.cpp | 11 +++++++---- lib/algo/manifold_interp/MatrixInterpolator.cpp | 6 +++--- lib/algo/manifold_interp/VectorInterpolator.cpp | 2 +- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/lib/algo/manifold_interp/Interpolator.cpp b/lib/algo/manifold_interp/Interpolator.cpp index 69694a56a..b269ee44f 100644 --- a/lib/algo/manifold_interp/Interpolator.cpp +++ b/lib/algo/manifold_interp/Interpolator.cpp @@ -143,7 +143,8 @@ std::vector obtainRBFGradientToTrainingPoints( { for (int i = 0; i < parameter_points.size(); i++) { - rbfs.push_back(obtainRBFGradient(rbf, epsilon, point, parameter_points[i], index)); + rbfs.push_back(obtainRBFGradient(rbf, epsilon, point, parameter_points[i], + index)); } } else @@ -155,7 +156,7 @@ std::vector obtainRBFGradientToTrainingPoints( } double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, - const Vector & point2, const int index) + const Vector & point2, const int index) { Vector diff; point1.minus(point2, diff); @@ -172,14 +173,16 @@ double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, else if (rbf == "IQ") { //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)); + 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") { //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)); + res = -epsilon*epsilon*diff.item(index)/(std::sqrt(1.0+eps_norm_squared)* + (1.0+eps_norm_squared)); } return res; diff --git a/lib/algo/manifold_interp/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index 01478c04e..669dac1c7 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -522,8 +522,8 @@ std::shared_ptr MatrixInterpolator::interpolateMatrix( for (int i = 0; i < point.dim(); i++) { std::vector rbf = obtainRBFGradientToTrainingPoints(d_parameter_points, - d_interp_method, - d_rbf, d_epsilon, point, i); + d_interp_method, + d_rbf, d_epsilon, point, i); std::shared_ptr gradient_matrix(obtainLogInterpolatedMatrix(rbf)); // *interpolated_matrix += *d_rotated_reduced_matrices[d_ref_point]; @@ -531,7 +531,7 @@ std::shared_ptr MatrixInterpolator::interpolateMatrix( d_interpolation_gradient.push_back(gradient_matrix); } } - else + else { std::cout << "Interpolated gradients are only implemented for \"LS\" "; CAROM_VERIFY(d_interp_method == "LS"); diff --git a/lib/algo/manifold_interp/VectorInterpolator.cpp b/lib/algo/manifold_interp/VectorInterpolator.cpp index 3cdaae1b2..3afadb996 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.cpp +++ b/lib/algo/manifold_interp/VectorInterpolator.cpp @@ -129,7 +129,7 @@ std::shared_ptr VectorInterpolator::interpolate(const Vector & point) for (int i = 0; i < point.dim(); ++i) { std::vector rbf = obtainRBFGradientToTrainingPoints(d_parameter_points, - d_interp_method,d_rbf, d_epsilon, point, i); + d_interp_method,d_rbf, d_epsilon, point, i); std::shared_ptr gradient_vector(obtainLogInterpolatedVector(rbf)); d_interpolation_gradient.push_back(gradient_vector); } From a3a3471eb5b19e1b088da0cf5e66e387f83d5214 Mon Sep 17 00:00:00 2001 From: ptranq Date: Tue, 10 Jun 2025 21:34:07 -0700 Subject: [PATCH 10/11] Ran stylize again. But this time from within the scripts directory --- lib/algo/manifold_interp/Interpolator.h | 14 +- lib/algo/manifold_interp/MatrixInterpolator.h | 6 +- lib/algo/manifold_interp/VectorInterpolator.h | 4 +- unit_tests/test_interpolation_gradients.cpp | 250 +++++++++++------- 4 files changed, 170 insertions(+), 104 deletions(-) diff --git a/lib/algo/manifold_interp/Interpolator.h b/lib/algo/manifold_interp/Interpolator.h index 24a8ad8d5..e5adf837d 100644 --- a/lib/algo/manifold_interp/Interpolator.h +++ b/lib/algo/manifold_interp/Interpolator.h @@ -114,7 +114,7 @@ class Interpolator * parameters should be computed. */ bool d_compute_gradients; - + /** * @brief Gradient with respect to the parameters. Only exists after * interpolate has been run @@ -157,9 +157,9 @@ class Interpolator * @param[in] point The unsampled parameter point. */ std::vector obtainRBFToTrainingPoints( -const std::vector & parameter_points, -const std::string & interp_method, const std::string & rbf, double epsilon, -const Vector & point); + const std::vector & parameter_points, + 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 @@ -174,7 +174,7 @@ const Vector & point); * @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 + * @param[in] index The index of the parameter that we are differentiating against. */ std::vector obtainRBFGradientToTrainingPoints( @@ -191,11 +191,11 @@ std::vector obtainRBFGradientToTrainingPoints( influence. * @param[in] point1 The first point. * @param[in] point2 The second point. - * @param[in] index The index of the parameter that we are + * @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); + const Vector & point2, const int index); /** * @brief Compute the sum of the RBF weights. diff --git a/lib/algo/manifold_interp/MatrixInterpolator.h b/lib/algo/manifold_interp/MatrixInterpolator.h index 5717377c5..a983c60ee 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.h +++ b/lib/algo/manifold_interp/MatrixInterpolator.h @@ -83,7 +83,9 @@ class MatrixInterpolator : public Interpolator * @brief Returns the interpolated matrix's gradient. * */ - std::vector> getGradient(){return d_interpolation_gradient;} + std::vector> getGradient() { + return d_interpolation_gradient; + } private: @@ -162,7 +164,7 @@ 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 diff --git a/lib/algo/manifold_interp/VectorInterpolator.h b/lib/algo/manifold_interp/VectorInterpolator.h index 2eab3b291..b00ad3d50 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.h +++ b/lib/algo/manifold_interp/VectorInterpolator.h @@ -75,7 +75,9 @@ class VectorInterpolator : public Interpolator * @brief Returns the interpolated matrix's gradient. * */ - std::vector> getGradient(){return d_interpolation_gradient;} + std::vector> getGradient() { + return d_interpolation_gradient; + } private: diff --git a/unit_tests/test_interpolation_gradients.cpp b/unit_tests/test_interpolation_gradients.cpp index eaa0d89c0..666d6630c 100644 --- a/unit_tests/test_interpolation_gradients.cpp +++ b/unit_tests/test_interpolation_gradients.cpp @@ -90,57 +90,72 @@ TEST(InterpolationTest,GaussianMatrix) { bool SUCCESS = true; CAROM::Vector input1(2, false); - input1(0) = 0.04; input1(1) = 0.04; + input1(0) = 0.04; + input1(1) = 0.04; CAROM::Vector input2(2, false); - input2(0) = 0.04; input2(1) = 0.06; + input2(0) = 0.04; + input2(1) = 0.06; CAROM::Vector input3(2, false); - input3(0) = 0.06; input3(1) = 0.04; + input3(0) = 0.06; + input3(1) = 0.04; CAROM::Vector input4(2, false); - input4(0) = 0.06; input4(1) = 0.06; + input4(0) = 0.06; + input4(1) = 0.06; CAROM::Vector target(2,false); - target(0) = 0.05; target(1) = 0.05; + 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; + 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); + + 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(); + std::vector> gradient = + interpolator.getGradient(); double epsilon = 1.0e-7; CAROM::Vector perturb0(2, false); - perturb0(0) = epsilon; perturb0(1) = 0.0; + perturb0(0) = epsilon; + perturb0(1) = 0.0; perturb0 += target; CAROM::Vector perturb1(2, false); - perturb1(0) = 0.0; perturb1(1) = epsilon; + 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 interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); - std::shared_ptr fd_grad_interp0(new Matrix(*interp_perturbation0)); + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); *fd_grad_interp0 -= *A; - std::shared_ptr fd_grad_interp1(new Matrix(*interp_perturbation1)); + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); *fd_grad_interp1 -= *A; for(int i = 0; i < 2; ++i) @@ -161,17 +176,19 @@ TEST(InterpolationTest,GaussianMatrix) if(abs(fd - my) > 1.e-5) { SUCCESS = false; - std::cout << "fd = " << fd << ", my = " << my << "abs(fd-my) = " << abs(fd-my) << std::endl; + 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; + std::cout << "fd = " << fd << ", my = " << my << "abs(fd-my) = " << abs( + fd-my) << std::endl; break; - } + } } if(!SUCCESS) break; } @@ -181,57 +198,72 @@ TEST(InterpolationTest,InverseQuadraticMatrix) { bool SUCCESS = true; CAROM::Vector input1(2, false); - input1(0) = 0.04; input1(1) = 0.04; + input1(0) = 0.04; + input1(1) = 0.04; CAROM::Vector input2(2, false); - input2(0) = 0.04; input2(1) = 0.06; + input2(0) = 0.04; + input2(1) = 0.06; CAROM::Vector input3(2, false); - input3(0) = 0.06; input3(1) = 0.04; + input3(0) = 0.06; + input3(1) = 0.04; CAROM::Vector input4(2, false); - input4(0) = 0.06; input4(1) = 0.06; + input4(0) = 0.06; + input4(1) = 0.06; CAROM::Vector target(2,false); - target(0) = 0.05; target(1) = 0.05; + 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; + 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); + + 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(); + std::vector> gradient = + interpolator.getGradient(); double epsilon = 1.0e-7; CAROM::Vector perturb0(2, false); - perturb0(0) = epsilon; perturb0(1) = 0.0; + perturb0(0) = epsilon; + perturb0(1) = 0.0; perturb0 += target; CAROM::Vector perturb1(2, false); - perturb1(0) = 0.0; perturb1(1) = epsilon; + 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 interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); - std::shared_ptr fd_grad_interp0(new Matrix(*interp_perturbation0)); + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); *fd_grad_interp0 -= *A; - std::shared_ptr fd_grad_interp1(new Matrix(*interp_perturbation1)); + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); *fd_grad_interp1 -= *A; for(int i = 0; i < 2; ++i) @@ -253,14 +285,14 @@ TEST(InterpolationTest,InverseQuadraticMatrix) { 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; } @@ -271,57 +303,72 @@ TEST(InterpolationTest,InverseMultiQuadraticMatrix) { bool SUCCESS = true; CAROM::Vector input1(2, false); - input1(0) = 0.04; input1(1) = 0.04; + input1(0) = 0.04; + input1(1) = 0.04; CAROM::Vector input2(2, false); - input2(0) = 0.04; input2(1) = 0.06; + input2(0) = 0.04; + input2(1) = 0.06; CAROM::Vector input3(2, false); - input3(0) = 0.06; input3(1) = 0.04; + input3(0) = 0.06; + input3(1) = 0.04; CAROM::Vector input4(2, false); - input4(0) = 0.06; input4(1) = 0.06; + input4(0) = 0.06; + input4(1) = 0.06; CAROM::Vector target(2,false); - target(0) = 0.05; target(1) = 0.05; + 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; + 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); + + 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(); + std::vector> gradient = + interpolator.getGradient(); double epsilon = 1.0e-7; CAROM::Vector perturb0(2, false); - perturb0(0) = epsilon; perturb0(1) = 0.0; + perturb0(0) = epsilon; + perturb0(1) = 0.0; perturb0 += target; CAROM::Vector perturb1(2, false); - perturb1(0) = 0.0; perturb1(1) = epsilon; + 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 interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); - std::shared_ptr fd_grad_interp0(new Matrix(*interp_perturbation0)); + std::shared_ptr fd_grad_interp0(new Matrix( + *interp_perturbation0)); *fd_grad_interp0 -= *A; - std::shared_ptr fd_grad_interp1(new Matrix(*interp_perturbation1)); + std::shared_ptr fd_grad_interp1(new Matrix( + *interp_perturbation1)); *fd_grad_interp1 -= *A; for(int i = 0; i < 2; ++i) @@ -343,14 +390,14 @@ TEST(InterpolationTest,InverseMultiQuadraticMatrix) { 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; } @@ -361,56 +408,71 @@ TEST(InterpolationTest,Vector) { bool SUCCESS = true; CAROM::Vector input1(2, false); - input1(0) = 0.04; input1(1) = 0.04; + input1(0) = 0.04; + input1(1) = 0.04; CAROM::Vector input2(2, false); - input2(0) = 0.04; input2(1) = 0.06; + input2(0) = 0.04; + input2(1) = 0.06; CAROM::Vector input3(2, false); - input3(0) = 0.06; input3(1) = 0.04; + input3(0) = 0.06; + input3(1) = 0.04; CAROM::Vector input4(2, false); - input4(0) = 0.06; input4(1) = 0.06; + input4(0) = 0.06; + input4(1) = 0.06; CAROM::Vector target(2,false); - target(0) = 0.05; target(1) = 0.05; + 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; + 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); + + CAROM::VectorInterpolator interpolator(inputs, rotations, vectors, 0, "G", "LS", + 0.9, true); std::shared_ptr v = interpolator.interpolate(target); - std::vector> gradient = interpolator.getGradient(); + std::vector> gradient = + interpolator.getGradient(); double epsilon = 1.0e-7; CAROM::Vector perturb0(2, false); - perturb0(0) = epsilon; perturb0(1) = 0.0; + perturb0(0) = epsilon; + perturb0(1) = 0.0; perturb0 += target; CAROM::Vector perturb1(2, false); - perturb1(0) = 0.0; perturb1(1) = epsilon; + 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 interp_perturbation0 = interpolator.interpolate( + perturb0); + std::shared_ptr interp_perturbation1 = interpolator.interpolate( + perturb1); - std::shared_ptr fd_grad_interp0(new Vector(*interp_perturbation0)); + std::shared_ptr fd_grad_interp0(new Vector( + *interp_perturbation0)); *fd_grad_interp0 -= *v; - std::shared_ptr fd_grad_interp1(new Vector(*interp_perturbation1)); + std::shared_ptr fd_grad_interp1(new Vector( + *interp_perturbation1)); *fd_grad_interp1 -= *v; for(int i = 0; i < 2; ++i) @@ -421,20 +483,20 @@ TEST(InterpolationTest,Vector) 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; - } + 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); } From 75b79ba2aca8f13be760a6356b84e24c7d73add7 Mon Sep 17 00:00:00 2001 From: ptranq Date: Wed, 11 Jun 2025 21:58:25 -0700 Subject: [PATCH 11/11] Addressed Dylan's concerns on PR. --- lib/algo/manifold_interp/Interpolator.cpp | 15 ++++++--------- lib/algo/manifold_interp/MatrixInterpolator.cpp | 11 +++-------- lib/algo/manifold_interp/VectorInterpolator.cpp | 3 +-- unit_tests/test_interpolation_gradients.cpp | 2 +- 4 files changed, 11 insertions(+), 20 deletions(-) diff --git a/lib/algo/manifold_interp/Interpolator.cpp b/lib/algo/manifold_interp/Interpolator.cpp index b269ee44f..aa98bda91 100644 --- a/lib/algo/manifold_interp/Interpolator.cpp +++ b/lib/algo/manifold_interp/Interpolator.cpp @@ -149,8 +149,7 @@ std::vector obtainRBFGradientToTrainingPoints( } else { - std::cout << "Interpolated gradients are only implemented for \"LS\" "; - CAROM_VERIFY(interp_method == "LS"); + CAROM_ERROR("Interpolated gradients are only implemented for \"LS\""); } return rbfs; } @@ -160,19 +159,19 @@ double obtainRBFGradient(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 if (rbf == "G") { - //res = std::exp(-eps_norm_squared); + //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") { - //res = 1.0 / (1.0 + eps_norm_squared); + //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)); @@ -180,7 +179,7 @@ double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, // Inverse multiquadric RBF else if (rbf == "IMQ") { - //res = 1.0 / std::sqrt(1.0 + eps_norm_squared); + //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)); } @@ -188,8 +187,6 @@ double obtainRBFGradient(std::string rbf, double epsilon, const Vector & point1, return res; } - - double rbfWeightedSum(std::vector& rbf) { double sum = 0.0; @@ -205,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/MatrixInterpolator.cpp b/lib/algo/manifold_interp/MatrixInterpolator.cpp index 669dac1c7..37dd16e81 100644 --- a/lib/algo/manifold_interp/MatrixInterpolator.cpp +++ b/lib/algo/manifold_interp/MatrixInterpolator.cpp @@ -106,8 +106,7 @@ std::shared_ptr MatrixInterpolator::interpolate(const Vector & point, { if (d_compute_gradients) { - std::cout << "Gradients are only implemented for B or G"; - CAROM_VERIFY(false); + CAROM_ERROR("Gradients are only implemented for B or G"); } interpolated_matrix = interpolateSPDMatrix(point); } @@ -115,8 +114,7 @@ std::shared_ptr MatrixInterpolator::interpolate(const Vector & point, { if (d_compute_gradients) { - std::cout << "Gradients are only implemented for B or G"; - CAROM_VERIFY(false); + CAROM_ERROR("Gradients are only implemented for B or G"); } interpolated_matrix = interpolateNonSingularMatrix(point); } @@ -526,15 +524,12 @@ std::shared_ptr MatrixInterpolator::interpolateMatrix( d_rbf, d_epsilon, point, i); std::shared_ptr gradient_matrix(obtainLogInterpolatedMatrix(rbf)); - // *interpolated_matrix += *d_rotated_reduced_matrices[d_ref_point]; - d_interpolation_gradient.push_back(gradient_matrix); } } else { - std::cout << "Interpolated gradients are only implemented for \"LS\" "; - CAROM_VERIFY(d_interp_method == "LS"); + CAROM_ERROR("Interpolated gradients are only implemented for \"LS\""); } } diff --git a/lib/algo/manifold_interp/VectorInterpolator.cpp b/lib/algo/manifold_interp/VectorInterpolator.cpp index 3afadb996..516e5aafb 100644 --- a/lib/algo/manifold_interp/VectorInterpolator.cpp +++ b/lib/algo/manifold_interp/VectorInterpolator.cpp @@ -136,8 +136,7 @@ std::shared_ptr VectorInterpolator::interpolate(const Vector & point) } else { - std::cout << "Interpolated gradients are only implemented for \"LS\" "; - CAROM_VERIFY(d_interp_method == "LS"); + CAROM_ERROR("Interpolated gradients are only implemented for \"LS\""); } } diff --git a/unit_tests/test_interpolation_gradients.cpp b/unit_tests/test_interpolation_gradients.cpp index 666d6630c..3e14ac84a 100644 --- a/unit_tests/test_interpolation_gradients.cpp +++ b/unit_tests/test_interpolation_gradients.cpp @@ -1,6 +1,6 @@ /****************************************************************************** * - * Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC + * Copyright (c) 2013-2025, Lawrence Livermore National Security, LLC * and other libROM project developers. See the top-level COPYRIGHT * file for details. *