Skip to content

Conversation

@joehellmersNOAA
Copy link
Collaborator

What type of PR is this? (check all applicable)

  • Refactor
  • Feature
  • Bug Fix
  • Optimization
  • Example
  • Documentation

Description

This pull request is needed to fix the C++ weights, and does so by separating them into two new classes that are not part of the associated operators.

Related Issues & Documents

QA Instructions, Screenshots, Recordings

Please replace this line with instructions on how to test your changes, a note
on the devices and browsers this has been tested on, as well as any relevant
images for UI changes.

Added/updated tests?

_We encourage you to test all code included with MOLE, including examples.

  • Yes
  • No, and this is why: please replace this line with details on why tests
    have not been included
  • I need help with writing tests

Read Contributing Guide and Code of Conduct

[optional] Are there any post deployment tasks we need to perform?

No

[optional] What gif best describes this PR or how it makes you feel?

Copy link
Collaborator

@jbrzensk jbrzensk left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah @joehellmersNOAA !!
I put some comments, mostly have to do with spacing, using the arma:: before Armadillo things, and const before the constant declaration in the test.

Good for 1D, I think that 2D and 3D are not as trivial as we think, hence the 2D in MATLAB, and no 3D implementations! I would leave the 2D and 3D until there is a better understanding, and a need, for them.

b.at(0) = -1.0;
b.at(m+1) = 1.0;


Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nit: extra lineIDK

#include "gradient.h"
#include "weightsP.h"


Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nit: extra line


};


Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nit; extra blank line

#include "mole.h"
#include <gtest/gtest.h>

vec Q_baseline[3];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

arma::vec


vec Q_baseline[3];
vec P_baseline[3];
Real tolerance = 1.0e-08;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

constexpr

// Baseline weights - Generated from MatLab/Octave
Q_baseline[0] = { 1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000};
Q_baseline[1] = { 1.000000000,1.125064293,0.751414447,1.162097097,0.962852897,0.997142531,0.962852897,
1.162097097,0.751414447,1.125064293,1.000000000 };
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nit; If you break the line, line up the numbers. And stay consistent. Some of the lines line up with the parentheses, others line up with the brackets. I vote to line up the numbers, it looks nicer.

P_baseline[2] = { 0.315722926,1.390677390,0.629532984,1.234237262,0.919144328,1.009804636,
1.000880474,1.000880474,1.009804636,0.919144328,1.234237262,0.629532984,
1.390677390,0.315722926 };
for (int k : {2,4,6}) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

put a blank line before the for, for consistency

{
Gradient G(k, m, dx);

vec b(m+2);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

arma::vec

b.at(m+1) = 1.0;


sp_mat Gtranspose = G.t();
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

arma::sp_mat

for (int k : {2,4,6}) {
run_Qtest(k);
run_Ptest(k);
}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For consistency, you put spaces between the closing loops and function in the previous Qweights and Pweights. Do the same here.

Copy link
Collaborator

@aboada aboada left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, still working on the review. I recommend applying clang-format (see: ClangFormat Documentation) to fix any formatting issues.

#include <cassert>

#ifdef EIGEN
//#ifdef EIGEN
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are we commenting out this line? Is Eigen mandatory now?

}
#endif

vec Utils::spsolve_eigenQR(const sp_mat &A, const vec &b) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you elaborate on why QR factorization is needed here as opposed to LU? I’m assuming this is for computing P and Q. Are we expecting any ill-conditioned matrices?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Moreover, the only difference between the LU method and the QR one is:

1c1
< vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) {
---
> vec Utils::spsolve_eigenQR(const sp_mat &A, const vec &b) {
4c4
<   Eigen::SparseLU<Eigen::SparseMatrix<Real>, Eigen::COLAMDOrdering<int>> solver;
---
>   Eigen::SparseQR<Eigen::SparseMatrix<Real>, Eigen::COLAMDOrdering<int>> solver;

We should consider refactoring to avoid code duplication.

Copy link
Collaborator

@aboada aboada Dec 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we give this test a more descriptive name? Since we may add more functionality in the future, how will we know that test6 is associated with P and Q?
Changing the test name may incur additional updates to the instructions in the CMakeList. If that's too much, we can do it in a separate PR.


for (int k : {2,4,6}) {
run_Qtest(k);
run_Ptest(k);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unit tests should be independent, isolated, and ideally testing only one thing (focusing on a single responsibility). Could we split the tests for P and Q?


TEST(WeightTests, Accuracy) {

// Baseline weights - Generated from MatLab/Octave
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of hardcoding these values from Matlab/Octave, why don't we check that $G^Tp = b_{m+2}$ (and the same for $Q$), as shown in Corbino and Castillo?

We can also have a separate test analogous to examples/matlab_octave/integration1D.m.

1.000880474,1.009804636,0.919144328,1.234237262,0.629532984,1.390677390,0.315722926
};

for (int k : {2,4,6}) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To ensure robustness, we should also add tests that check for unexpected, unsupported, corner cases, etc.

For instance, what happens with negative k when calling the WeightsQ constructor (k,m,dx). What about odd values? Negative m or dx, etc.

@joehellmersNOAA
Copy link
Collaborator Author

Thanks, still working on the review. I recommend applying clang-format (see: ClangFormat Documentation) to fix any formatting issues.

I have no problem passing this C++ code through some sort of linter, but others should weigh in if this is the best one or not. @jbrzensk @cpaolini

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants