Add gmrf and precision builder functions as precursors for eventually adding DSEM - #1594
Add gmrf and precision builder functions as precursors for eventually adding DSEM#1594e-perl-NOAA wants to merge 32 commits into
Conversation
🎨 Chore: code formatting workflowOur automated workflows cannot run on forks because of permission issues, and thus, we ask that you run the following code locally and push any changes that are created to your feature branch. You will only be reminded of this once per PR. Thank you! Format C++ code
Format R code
styler::style_pkg() # Style R code
roxygen2::roxygenise() # Update documentation
styler::style_pkg() # Style R code again
roxygen2::roxygenise() # Update documentation again
usethis::use_tidy_description() # Style DESCRIPTION filePush changes
|
Wow 😮 , you did it! I haven't taken a close look at the code in those files yet, but for For the I'm not sure what your timeline is for working on the tests, but I'll be setting aside dedicated time each week to work on testing-related issues. Feel free to join those sessions (starting July 9 and listed on the FIMS calendar) if you have any questions or would like to work through any testing challenges together! |
|
@Bai-Li-NOAA Yeah, that sounds great! I won't have time to work on tests this week, but I'll definitely have questions for you when I start working on it next week! Also, I love that you are leaning into me naming my AI/agents Fernando 🤣 |
|
@nathanvaughan-NOAA can you review this PR? |
|
@kellijohnson-NOAA I need to add tests, I just haven't had time to work on them 😞 |
|
No worries @e-perl-NOAA, Fernando is a busy guy I'm sure he'll get to them eventually 😁 I'll just go through the functions. |
|
@e-perl-NOAA, there are two issues:
I started working on a comment that included a bunch of code - I realized it might be easier to set up a co-working session to go over these changes together. I think Fernando has gotten you to a nice spot, but at this point, writing more instructions to AI would probably take more time than implementing the needed changes directly! |
| } | ||
|
|
||
| // Register builder in Information Map (to be accessed by GMRF distribution) | ||
| info->dsem_builders[this->id] = builder; |
There was a problem hiding this comment.
dsem_builders does not exist in information. More code is needed to link the precision matrix builder to the GMRF both here and in information.
| }; | ||
|
|
||
| /** | ||
| * @brief Rcpp interface for DSEM growth. |
There was a problem hiding this comment.
remove beta_z from the next four files (rcpp_growth, rcpp_maturity, rcpp_recruitment, and rcpp_selectivity)
| } | ||
|
|
||
| // Centering: x - mu because TMB's GMRF expects input centered around a mean of 0. | ||
| vector<Type> x_centered(n_x); |
There was a problem hiding this comment.
I think this needs to be std::vector
There was a problem hiding this comment.
It would also be faster to have
vector x_centered;
x_centered.reserve(n_x);
for (size_t i = 0; i < n_x; ++i) {
x_centered.emplace_back(this->get_observed(i) - this->get_expected(i));
}
this avoids filling all the x_centered values with default first and overwriting them, emplace also avoids creating and then moving the centered value which push_back does.
There was a problem hiding this comment.
@nathanvaughan-NOAA
This is what implementing your suggestion gave me, does this look right?
// Centering: x - mu. TMB's GMRF expects input centered around a mean of 0.
// To improve performance, we build a std::vector using emplace_back to avoid
// default-constructing and then reassigning every element.
std::vector<Type> x_centered_std;
x_centered_std.reserve(n_x);
for (size_t i = 0; i < n_x; ++i) {
x_centered_std.emplace_back(this->get_observed(i) - this->get_expected(i));
}
// Evaluate TMB GMRF and multiply by -1 to convert from negative log-likelihood to log-likelihood.
this->lpdf = -1.0 * density::GMRF(*(this->precision_matrix_ptr))(Eigen::Map<const vector<Type>>(x_centered_std.data(), n_x));
|
|
||
| /** | ||
| * @brief Rcpp interface for DSEM growth. | ||
| */ |
There was a problem hiding this comment.
We should avoid setting this up to create module specific DSEM classes. Everything should be handled in the dsem interface so it can then point to arbitrary variables in any current/future modules without adding development overhead to them.
| int type = 0; /**< 1 = A path effect (Rho), 2 = Variance (Gamma) */ | ||
| int from = 0; /**< The variable the arrow starts from */ | ||
| int to = 0; /**< The variable the arrow points to */ | ||
| int beta_index = 0; /**< Which parameter in beta_z to use for this arrow */ |
There was a problem hiding this comment.
Am I right that using the index here is basically so you can map off multiple linkages to be estimated together? I'm not sure how common that would be but it may make more sense just to have the value in RAMPath and allow it to be estimable and/or mapped through TMB like other parameters? That basically replaces beta_index and start with beta_z that can be a fixed or estimated variable.
| } | ||
|
|
||
| // Poke the "strength" into the correct grid slot based on arrow type. | ||
| if (this->paths[r].type == 1) { |
There was a problem hiding this comment.
When I looked into SparseMatrix it looks like using .coeffRef is supposedly really slow due to memory moving issues and using Eigen::Triplet is much faster. Something like this apparently
int rows = 3;
int cols = 3;
Eigen::SparseMatrix<double> A(rows, cols);
std::vector<Eigen::Triplet<double>> tripletList;
tripletList.reserve(4);
tripletList.push_back(Eigen::Triplet<double>(0, 0, 1.5));
tripletList.push_back(Eigen::Triplet<double>(1, 2, 2.0));
tripletList.push_back(Eigen::Triplet<double>(2, 1, -3.5));
tripletList.push_back(Eigen::Triplet<double>(2, 2, 4.0));
A.setFromTriplets(tripletList.begin(), tripletList.end());
A.makeCompressed();
There was a problem hiding this comment.
@nathanvaughan-NOAA When Fernando and I try to implement your suggestion it is giving me the following. Does this look right to you?
// 2. Initialize sparse components
// We create three empty grids that are "sparse"
Eigen::SparseMatrix<Type> Rho_kk(static_cast<int>(n_k), static_cast<int>(n_k)); // Grid for causal paths (from -> to)
Eigen::SparseMatrix<Type> Gamma_kk(static_cast<int>(n_k), static_cast<int>(n_k)); // Grid for variances (<->)
Eigen::SparseMatrix<Type> I_kk(static_cast<int>(n_k), static_cast<int>(n_k)); // "Standard" grid (identity)
I_kk.setIdentity(); // Fill diagonal with 1s
// 3. Translate the "arrows" (RAMPath) into triplet lists for efficient sparse matrix construction.
std::vector<Eigen::Triplet<Type>> rho_triplets;
std::vector<Eigen::Triplet<Type>> gamma_triplets;
rho_triplets.reserve(this->paths.size());
gamma_triplets.reserve(this->paths.size());
for (size_t r = 0; r < this->paths.size(); ++r) {
// C++ starts counting at 0, but R starts at 1, so we subtract 1.
const int from = this->paths[r].from - 1;
const int to = this->paths[r].to - 1;
// Check if the user's arrow points to a slot that doesn't exist.
if (from < 0 || to < 0 || static_cast<size_t>(from) >= n_k ||
static_cast<size_t>(to) >= n_k) {
throw std::invalid_argument(
"DSEMPrecisionMatrixBuilder: RAM indices out of bounds.");
}
// Determine the "strength" of this arrow. If beta_index is 1 or
// more, it's a parameter the model is guessing, otherwise, it's a fixed number.
Type value = this->paths[r].start;
if (this->paths[r].beta_index >= 1) {
const size_t b_idx = static_cast<size_t>(this->paths[r].beta_index - 1);
if (b_idx >= this->beta_z.size()) {
throw std::invalid_argument(
"DSEMPrecisionMatrixBuilder: beta_index points past beta_z size.");
}
value = this->beta_z[b_idx];
}
// Add the "strength" to the correct triplet list based on arrow type.
if (this->paths[r].type == 1) {
rho_triplets.push_back(Eigen::Triplet<Type>(from, to, value));
} else if (this->paths[r].type == 2) {
gamma_triplets.push_back(Eigen::Triplet<Type>(from, to, value));
}
}
Rho_kk.setFromTriplets(rho_triplets.begin(), rho_triplets.end());
Gamma_kk.setFromTriplets(gamma_triplets.begin(), gamma_triplets.end());
There was a problem hiding this comment.
Yep this looks good to me, could probably use emplace_back here too to avoid the move but that is also a wider refactor question throughout FIMS.
| } | ||
|
|
||
| // Centering: x - mu because TMB's GMRF expects input centered around a mean of 0. | ||
| vector<Type> x_centered(n_x); |
There was a problem hiding this comment.
It would also be faster to have
vector x_centered;
x_centered.reserve(n_x);
for (size_t i = 0; i < n_x; ++i) {
x_centered.emplace_back(this->get_observed(i) - this->get_expected(i));
}
this avoids filling all the x_centered values with default first and overwriting them, emplace also avoids creating and then moving the centered value which push_back does.
This file implements the GMRF distribution functor, which evaluates the log-likelihood of Gaussian Markov Random Fields using a sparse precision matrix.
…::vector with emplace back in gmrf function
…of how parameters are defined
…e reworked anyways
|
@Andrea-Havron-NOAA Also, gemini is convinced that the following lines that we added yesterday in the rcpp_precision_builders.hpp need to be removed: because of the following: Does this seem right? |
|
@e-perl-NOAA should we mark this as a draft PR until you are ready for a final review? Or, is there something you currently need help with? |
|
@kellijohnson-NOAA I'm still waiting on some feedback from @Andrea-Havron-NOAA and/or @nathanvaughan-NOAA after my last two posts in the thread that contained some questions. |
|
Andrea is the expert on the GMRF stats but I’m good with explicitly filling out the symmetrical matrix for clarity.
|
This is a bit philosophical. The GMRF is a distribution that requires a precision matrix. Rather than initialize the GMRF and its precision matrix independently and add a lot of checks to make sure the user set up both correctly and linked them together correctly, I think they should be linked internally. I decided that the user specifies the precision matrix (in this case, dsem), and by default, this initializes GMRF. We could rename the interface to dsemGMRFInterface to be more explicit about what is happening under the hood. Given this, there are a couple of pieces missing as Gemini is right in that the gmrf is created locally and then destroyed. We can make the gmrf object persist by adding this code just below Also And then in the gmrf.hpp file, add: There is still a missing piece, we will need to link the distribution's I did want to mention, another design direction would be to create a gmrf interface which creates the DSEM builder instead given a precision type input from the user. I guess my rational for thinking about it the other way is that the parameters are associated with the precision matrix, so it is easier to set them up for registration if the primary class being initialized is the precision matrix which then calls the GMRF rather than the other way around. I will look over the suggested changes to precision_builders.hpp tomorrow! |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1594 +/- ##
==========================================
- Coverage 84.88% 83.73% -1.15%
==========================================
Files 105 56 -49
Lines 9473 2232 -7241
Branches 536 539 +3
==========================================
- Hits 8041 1869 -6172
+ Misses 1395 298 -1097
- Partials 37 65 +28 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
| std::shared_ptr<fims_distributions::PrecisionMatrixBuilderBase<Type>> precision_matrix_ptr = nullptr; | ||
|
|
||
| // IS THIS WHERE THIS GOES? | ||
| Eigen::SparseMatrix<Type> Q = precision_matrix_ptr->BuildPrecisionMatrixSparse(); |
There was a problem hiding this comment.
This needs to be inside the GMRF evaluate function
Co-authored-by: Andrea-Havron-NOAA <Andrea-Havron-NOAA@users.noreply.github.com>
5 participants
Disclaimer
I don't actually know what I'm doing here, but there was enthusiasm for progress being made on getting DSEM into FIMS so I embarked on this journey ⛵ with nothing but Copilot and a dream ☁️. Please keep that in mind as you make comments and I beg you to explain things in a way that you would see them in those "XXX for Dummies" books.
What is the feature?
Getting the behind the scenes functions of gmrf and a precision builder (that can be used with more than just dsem) that are precursors seemed like a good point to start merging things in since myself and copilot/notebookLM then all started to get confused on how to proceed from here (it started going a bit off the rails once we started getting to what needed to be added in the
information.hppfile) . I don't even think that these functions are "exposed" enough to write tests for? Maybe I'm wrong here and if so, @Bai-Li-NOAA I would DEFINITELY need your help.Once it seems like there is consensus that these are a good way forward and/or good to go, they can be merged in and I/we can incrementally work on getting DSEM into FIMS in subsequent PRs. Hopefully that method confuses copilot/notebookLM (and honestly, myself) less.
How have you implemented the solution?
Instructions for code reviewer
👋Hello reviewer👋, thank you for taking the time to review this PR!
nit:(for nitpicking) as the comment type. For example,nit:I prefer using adata.frame()instead of amatrixbecause ...This PR is now ready to be merged.Checklist