feat(edm): implement EDM prediction framework (Simplex, S-Map, GP-EDM)

[aman-raj-srivastva]

Summary

This PR implements the core EDM prediction framework for Empirical Dynamic Modeling (EDM) in FIMS as part of the GSoC 2026 project:
"Add Empirical Dynamic Modeling to FIMS."

This work builds on the delay embedding infrastructure introduced in PR #1528 and implements the prediction algorithms described in Issue #1487.

Target Branch: main-edm

This PR completes the standalone prediction framework for EDM. Integration of these prediction methods into the FIMS likelihood/TMB pipeline will be addressed in a subsequent PR.

---

What this PR adds

Generic EDM Prediction Framework

Implemented a reusable and extensible prediction framework for EDM algorithms.

Features include:

• Common prediction interface
• Shared prediction result structures
• Extensible architecture for additional EDM prediction methods
• Configurable forecast horizon support (shared across all prediction methods)
• Modular design for future integration with the FIMS objective function

---

Simplex Projection

Implemented the Standard Simplex Projection algorithm (Sugihara & May, 1990).

Features include:

• Nearest-neighbor based prediction
• Distance-weighted forecasting
• Integration with the delay embedding infrastructure

---

Distance Weighting Infrastructure

Implemented reusable distance weighting utilities shared across prediction algorithms.

Features include:

• Squared Euclidean distance calculation utilities
• Modular exponential and Gaussian weighting functions
• Shared infrastructure for Simplex, S-Map, and GP-EDM

Note on distance metric: This implementation uses squared Euclidean distances in all kernel functions (rather than Euclidean distances as in rEDM). This is an intentional design decision to maintain compatibility with CppAD automatic differentiation, which avoids the sqrt() operation. Predictions are numerically verified against pure-R squared-distance reference implementations.

---

Standard S-Map

Implemented the Standard S-Map prediction algorithm (Sugihara, 1994).

Features include:

• Weighted local linear regression
• Configurable exponential and Gaussian kernel functions
• Shared linear algebra utilities (Gaussian elimination with partial pivoting)
• Integration with delay embeddings

---

GP-EDM

Implemented the GP-EDM prediction algorithm following the reference GPEDM implementation (Munch et al. 2017; Rogers 2023).

Features include:

• Gaussian Process regression on delay embeddings
• Automatic Relevance Determination (ARD) length-scale parameters (phi)
• ARD squared-exponential kernel with half-Normal priors on phi
• Rprop optimizer for marginal likelihood maximization
• Shared covariance kernel and linear algebra utilities
• Rcpp interfaces for both hyperparameter fitting and prediction

Bug fix: Corrected a missing 0.5 factor on the log-determinant term in the GP log-marginal likelihood. This was identified during numerical verification against the GPEDM reference package.

Known intentional differences from the GPEDM reference package:

Aspect	GPEDM Reference	This Implementation	Reason
Matrix inversion	Cholesky (chol2inv)	Gaussian elimination	Simpler for v1; both are correct for SPD matrices
rho parameter	Estimated (multi-population)	Not estimated	FIMS single-population case only
Input scaling	global/local/none	Not implemented (user responsibility)	Scope reduction for v1

---

Rcpp Interfaces

Implemented R-callable interfaces for all prediction methods via Rcpp modules.

Includes:

• SimplexProjectionInterface — fields: embedding_dimension, n_neighbors
• SMapProjectionInterface — fields: embedding_dimension, theta, kernel
• GPEdmProjectionInterface — fields: embedding_dimension, phi, sigma2, ve; methods: fit(), predict()
• to_matrix() helper bridging DelayEmbeddingInterface storage to the C++ prediction framework
• All interfaces registered with EDMInterfaceBase::live_objects for Rcpp memory safety

All three interfaces are callable from R via methods::new().

---

Validation

Numerical verification was performed for all three prediction methods:

Algorithm	Reference	Max Difference
Simplex Projection	Pure-R squared-distance implementation	2.22e-16 (machine ε)
S-Map	Pure-R WLS squared-distance implementation	1.30e-14
GP-EDM (fixed params)	GPEDM R package (predict.GP)	2.53e-10

For GP-EDM, the sigma2 and ve gradients were additionally verified to match the reference package exactly (these do not depend on the per-dimension distance matrices), confirming the covariance matrix math is correct.

---

Testing

This PR includes comprehensive testing for the EDM prediction framework.

Current coverage:

• 114 GoogleTests (C++)
• 87 testthat tests (R)

Testing covers:

• Prediction framework behavior
• Simplex Projection (IO correctness, numerical agreement, edge cases, error handling)
• Standard S-Map (IO correctness, numerical agreement, kernel comparison, error handling)
• GP-EDM (IO correctness, fit → predict chain, hyperparameter localization, error handling)
• Rcpp interfaces (construction, field access, invalid ID handling)
• Numerical agreement with pure-R reference implementations

---

Future Work

This PR completes the standalone EDM prediction framework. Future work can focus on:

• Integrating the prediction framework into the FIMS likelihood/TMB pipeline
• CppAD/TMB automatic differentiation support throughout the EDM module
• Observation uncertainty propagation through the delay embedding
• Additional statistical model integration

---

Related Issues

• Addresses Implement Generic EDM Prediction Functor #1487
• Builds on Implement Delay Embedding Generator for Time Series #1484
• Builds on Add EDM delay embedding infrastructure and tests #1528

[aman-raj-srivastva]

Reverted the S-Map commit temporarily for further review before re-pushing

[aman-raj-srivastva]

@stevemunch @nathanvaughan-NOAA while implementing the S-Map predictor, I came across a few design questions where I'd appreciate feedback before I push the implementation.

I have currently implemented the standard S-Map formulation described in Section 2.1 of Esguerra & Munch (2024), i.e., the local linear model that serves as the baseline/M-step component of the full HMS-map.

For each local neighborhood, the implementation constructs a design matrix:

$$X_i = [1, x_t, x_{t-\tau}, \ldots, x_{t-(E-1)\tau}]$$

and solves the weighted least-squares system:

$$(X^T W X)\beta = X^T W y$$

using a CppAD-safe Gaussian elimination routine.

Before finalizing the implementation, I wanted to get your thoughts on a couple of design choices.

1. Weighting Kernel / Distance Metric

Classic S-Map implementations (e.g., Sugihara 1994 and rEDM) typically use an exponential kernel based on Euclidean distance:

$$w_i = \exp\left(-\theta \frac{d_i}{\bar d}\right)$$

However, Section 2.1 of Esguerra & Munch (2024) defines the weighting function as:

$$w_i = \exp\left(-\theta^2 \left(\frac{d_i}{D}\right)^2\right)$$

which is effectively a Gaussian kernel based on squared Euclidean distance.

For the generic FIMS S-Map predictor, would you prefer:

• the formulation used in the HMS-map paper, or
• compatibility with the more traditional rEDM-style weighting?

Using squared distances has the additional benefit of avoiding a sqrt() inside the AD tape.

2. Target Alignment in Delay Embedding

The paper formulates the delay embedding map as a one-step-ahead forecast:

$$x_{t+1} = f(x_t, \ldots, x_{t-E+1})$$

While wiring the implementation, I noticed that the current DelayEmbeddingMatrix assigns target_values[row] and embedded_values[row][0] to the same observation, which appears to correspond to a 0-step-ahead target.

Am I interpreting this correctly, or is the intended target shift handled elsewhere in the workflow?

If not, would aligning the embedding construction with the standard one-step-ahead formulation be the preferred approach?

3. Scope of This PR

My understanding is that this PR should focus on the standard S-Map predictor only:

• local weighted linear regression
• prediction interface integration
• testing and validation infrastructure

The EM iteration and state-estimation machinery required for the full HMS-map would then be a separate follow-up effort.

Just wanted to confirm that this matches the intended project scope before I proceed further.

Thanks! I currently have the S-Map implementation and associated tests working locally, and I wanted to verify these design details before pushing the next commit.

[stevemunch]

Aman-In my experience using distance v. squared distance in the weighting kernel doesn’t make a huge difference in performance except in rare cases.  I’d use whatever was easier/faster.  In the state-estimation step of the HMS algorithm, we are updating all of the states based on the current estimates of the S-map coefficients and the observations.  There’s no time lag in this step.  In the S-map step, the model is fit one step ahead as usual.  Regarding HMS map being a separate project - I will let others chime in.  For my two cents - that’s the part that requires the most care in implementation. The rest is pretty straightforward. Critically ‘standard EDM’ does not handle observation error explicitly.  The statistically minded folks at NMFS generally care about observation uncertainty, so it would be most valuable to if you completed the HMS map code or some other algorithm for handling observation uncertainty.Sent from my iPadOn Jun 24, 2026, at 12:04 PM, Aman Raj ***@***.***> wrote:aman-raj-srivastva left a comment (NOAA-FIMS/FIMS#1561) @stevemunch @nathanvaughan-NOAA while implementing the S-Map predictor, I came across a few design questions where I'd appreciate feedback before I push the implementation. I have currently implemented the standard S-Map formulation described in Section 2.1 of Esguerra & Munch (2024), i.e., the local linear model that serves as the baseline/M-step component of the full HMS-map. For each local neighborhood, the implementation constructs a design matrix: $$X_i = [1, x_t, x_{t-\tau}, \ldots, x_{t-(E-1)\tau}]$$ and solves the weighted least-squares system: $$(X^T W X)\beta = X^T W y$$ using a CppAD-safe Gaussian elimination routine. Before finalizing the implementation, I wanted to get your thoughts on a couple of design choices. 1. Weighting Kernel / Distance Metric Classic S-Map implementations (e.g., Sugihara 1994 and rEDM) typically use an exponential kernel based on Euclidean distance: $$w_i = \exp\left(-\theta \frac{d_i}{\bar d}\right)$$ However, Section 2.1 of Esguerra & Munch (2024) defines the weighting function as: $$w_i = \exp\left(-\theta^2 \left(\frac{d_i}{D}\right)^2\right)$$ which is effectively a Gaussian kernel based on squared Euclidean distance. For the generic FIMS S-Map predictor, would you prefer: the formulation used in the HMS-map paper, or compatibility with the more traditional rEDM-style weighting? Using squared distances has the additional benefit of avoiding a sqrt() inside the AD tape. 2. Target Alignment in Delay Embedding The paper formulates the delay embedding map as a one-step-ahead forecast: $$x_{t+1} = f(x_t, \ldots, x_{t-E+1})$$ While wiring the implementation, I noticed that the current DelayEmbeddingMatrix assigns target_values[row] and embedded_values[row][0] to the same observation, which appears to correspond to a 0-step-ahead target. Am I interpreting this correctly, or is the intended target shift handled elsewhere in the workflow? If not, would aligning the embedding construction with the standard one-step-ahead formulation be the preferred approach? 3. Scope of This PR My understanding is that this PR should focus on the standard S-Map predictor only: local weighted linear regression prediction interface integration testing and validation infrastructure The EM iteration and state-estimation machinery required for the full HMS-map would then be a separate follow-up effort. Just wanted to confirm that this matches the intended project scope before I proceed further. Thanks! I currently have the S-Map implementation and associated tests working locally, and I wanted to verify these design details before pushing the next commit. —Reply to this email directly, view it on GitHub, or unsubscribe.Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

[nathanvaughan-NOAA]

Hey @aman-raj-srivastva, great job on your progress so far. For these questions I would suggest

1. Can you add the option for the user to choose the kernal function? That seems like it would offer the most utility for comparing with other approaches and leave options open for different functions to be added in the future. @stevemunch should be able to explain the rational for the change in kernal. I think the main difference is just that the guassian has fatter tails than the exponential but I thought both methods generally estimated a scale factor that may minimize the impact of the kernal choice. I just saw @stevemunch above say basically that same thing.

2. Good catch on the indexing I missed that the first time. I think it's best to change that to the 1 step ahead index. From @stevemunch 's comment above the state-estimation step would be updating the base values (why it's so beneficial to set this all up as linked pointers) and then the S-MAP step would use the embedded and target objects to make predictions.

3. @stevemunch is correct that the HMS uncertainty integration is critical to practically applying this and an important feature, but the mechanics will be complex and I think fairly isolated from this initial structure development so I think it's fine to compartmentalize for now and keep this pull request focused.

[aman-raj-srivastva]

Hello everyone!

I've pushed new commits that completes this implementation phase of my GSoC proposal.

@stevemunch @nathanvaughan-NOAA I had a few questions before moving forward:

1. During validation, I noticed what appears to be an issue in the reference GPEDM R package that causes the length-scale gradients to remain zero during fitting. Would it be helpful to document this discrepancy in the PR description or code comments for future reference, or would it be better to open a separate issue against the GPEDM package?

2. Is the current to_matrix() helper (used to bridge DelayEmbeddingInterface to DelayEmbeddingMatrix) the preferred pattern for FIMS, or would you recommend a different approach for future TMB integration?

3. For GP-EDM::fit(), the optimized hyperparameters are currently written back to the interface object (phi, sigma2, ve). Would you prefer this behaviour, or should fit() simply return the optimized values without modifying the object state?

Thank you! I'd really appreciate any feedback or suggestions before moving on to the next phase.

[stevemunch]

Aman-Some more details on the gradient issue would be really helpful.  SteveSent from my iPadOn Jul 7, 2026, at 5:04 AM, Aman Raj ***@***.***> wrote:aman-raj-srivastva left a comment (NOAA-FIMS/FIMS#1561) Hello everyone! I've pushed new commits that completes this implementation phase of my GSoC proposal. @stevemunch @nathanvaughan-NOAA I had a few questions before moving forward: During validation, I noticed what appears to be an issue in the reference GPEDM R package that causes the length-scale gradients to remain zero during fitting. Would it be helpful to document this discrepancy in the PR description or code comments for future reference, or would it be better to open a separate issue against the GPEDM package? Is the current to_matrix() helper (used to bridge DelayEmbeddingInterface to DelayEmbeddingMatrix) the preferred pattern for FIMS, or would you recommend a different approach for future TMB integration? For GP-EDM::fit(), the optimized hyperparameters are currently written back to the interface object (phi, sigma2, ve). Would you prefer this behaviour, or should fit() simply return the optimized values without modifying the object state? Thank you! I'd really appreciate any feedback or suggestions before moving on to the next phase. —Reply to this email directly, view it on GitHub, or unsubscribe.Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

[aman-raj-srivastva]

Hi @stevemunch! Happy to elaborate.

While validating the GP-EDM implementation against the reference package, I traced the optimization routine for the ARD length-scale parameters (phi). These gradients depend on the per-dimension squared-distance matrices, where for each embedding dimension (d),

$D_d(i,j) = \left(x_{id} - x_{jd}\right)^2$

In the current GPEDM R implementation, these matrices are constructed inside fmingrad_Rprop() as:

for (i in 1:ncol(X)) {
  D[[i]] <- laGP::distance(X[, i])
}

From what I observed, X[, i] is simplified to a vector (the default behavior of [ in R). When this vector is passed to laGP::distance(), it appears to be interpreted as a single observation rather than an N × 1 matrix, resulting in a 1 × 1 distance matrix instead of the expected pairwise N × N distance matrix for that embedding dimension.

Since the ARD gradient is computed from these distance matrices,

dC <- -D[[i]] * Cd
dl[i] <- vQ %*% as.vector(dC)

the likelihood contribution to the gradient appears to become zero, leaving the optimization of the length-scale parameters largely driven by the prior.

In the FIMS implementation, I compute the pairwise squared-distance matrices explicitly for each embedding dimension in C++, so each ARD parameter contributes correctly to the covariance matrix and its gradient.

As a sanity check, the gradients for sigma2 and ve (which do not depend on these distance matrices) match the reference implementation, while only the phi gradients differ. After correcting this in FIMS, the optimized length scales differ from the reference package and produce a slightly better log-marginal likelihood on the logistic-map example that I used for validation.

If it would be helpful, I can also prepare a minimal reproducible example demonstrating the behavior, or open an issue against the GPEDM package if we agree this is an issue in the reference implementation.

[stevemunch]

Aman,I don’t really do R, so am not 100% sure of the problem.  But if the only thing that matters is the prior, the inverse length scale parameters should all go to zero.  Similarly, if the gradient is a scalar (and applied to all the length scales), then all length scales should end up the same, since they have the same prior. So, I don’t think I understand the issue still. Sorry.You could certainly post this as an issue to GPEDM and see what Tanya says.  Sent from my iPhoneOn Jul 7, 2026, at 10:16 AM, Aman Raj ***@***.***> wrote:aman-raj-srivastva left a comment (NOAA-FIMS/FIMS#1561) Hi @stevemunch! Happy to elaborate. While validating the GP-EDM implementation against the reference package, I traced the optimization routine for the ARD length-scale parameters (phi). These gradients depend on the per-dimension squared-distance matrices, where for each embedding dimension (d), $D_d(i,j) = \left(x_{id} - x_{jd}\right)^2$ In the current GPEDM R implementation, these matrices are constructed inside fmingrad_Rprop() as: for (i in 1:ncol(X)) { D[[i]] <- laGP::distance(X[, i]) } From what I observed, X[, i] is simplified to a vector (the default behavior of [ in R). When this vector is passed to laGP::distance(), it appears to be interpreted as a single observation rather than an N × 1 matrix, resulting in a 1 × 1 distance matrix instead of the expected pairwise N × N distance matrix for that embedding dimension. Since the ARD gradient is computed from these distance matrices, dC <- -D[[i]] * Cd dl[i] <- vQ %*% as.vector(dC) the likelihood contribution to the gradient appears to become zero, leaving the optimization of the length-scale parameters largely driven by the prior. In the FIMS implementation, I compute the pairwise squared-distance matrices explicitly for each embedding dimension in C++, so each ARD parameter contributes correctly to the covariance matrix and its gradient. As a sanity check, the gradients for sigma2 and ve (which do not depend on these distance matrices) match the reference implementation, while only the phi gradients differ. After correcting this in FIMS, the optimized length scales differ from the reference package and produce a slightly better log-marginal likelihood on the logistic-map example that I used for validation. If it would be helpful, I can also prepare a minimal reproducible example demonstrating the behavior, or open an issue against the GPEDM package if we agree this is an issue in the reference implementation. —Reply to this email directly, view it on GitHub, or unsubscribe.Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you were mentioned.Message ID: ***@***.***>

[aman-raj-srivastva]

Aman,I don’t really do R, so am not 100% sure of the problem. But if the only thing that matters is the prior, the inverse length scale parameters should all go to zero. Similarly, if the gradient is a scalar (and applied to all the length scales), then all length scales should end up the same, since they have the same prior. So, I don’t think I understand the issue still. Sorry.You could certainly post this as an issue to GPEDM and see what Tanya says.

Thanks for taking a closer look. Your observation helped me identify where my interpretation went wrong.

Looking back, the zero phi gradients I observed were actually an artifact of my verification code rather than the GPEDM optimizer itself. Specifically, fitGP() does not retain the per-dimension distance matrices (gp_fit_ref$inputs$D is NULL), so when I later called getlikegrad(..., D = NULL) manually, the phi gradients evaluated to zero. That behavior was not occurring during the actual optimization.

As you pointed out, if the optimizer were truly receiving zero (or identical) gradients for all length-scale parameters, we would expect the fitted phi values to collapse toward the prior or converge to very similar values, which is not what the package produces. So I no longer believe there is an issue in the reference implementation.

Based on this, I believe the FIMS implementation is behaving as intended. Unless there are any additional suggestions or changes you'd like to see for this PR, I'll plan to move on to the next stage of the project, focusing on integrating the EDM prediction framework into the FIMS likelihood/TMB pipeline while continuing to address any review feedback on this PR.

Thanks again!

[mollystevens-noaa]

Hey @aman-raj-srivastva ! @nathanvaughan-NOAA and I rebased main-edm to main and resolved merge conflicts on your pull request. There are currently tests failing, but it could have been due to our rebase. We created a safety branch called 'main-edm-backup' before rebasing just in case. In the future, we can pull your changes in first (open to workflow suggestions here). I can rebase to main on a schedule (e.g. every Monday morning) so you know when to expect them. Hopefully this didn't break anything for you! We're holding off on merging the pull request until this passes all the tests, please let us know if you'd like a working meeting to unravel everything.

[aman-raj-srivastva]

Hi @mollystevens-noaa and @nathanvaughan-NOAA,
From what I understood so far is that main-edm will periodically rebased onto main (e.g. on Mondays) so it stays consistent with the latest changes.
Then, for any active feature branch (such as my current integration work), I will first update my local main-edm to match the latest upstream/main-edm, and then rebase my feature branch to the updated main-edm before continuing development.
With the recent main-edm rebase, rebasing my feature/edm-fims-integration branch seems to replay the earlier delay embedding commits. Would you recommend creating a new branch from the updated main-edm or continuing with the current branch?
If that sounds right, just let me know. If not, I'm happy to discuss it in a working meeting.

[mollystevens-noaa]

Hey @aman-raj-srivastva ! Please continue working in the current branch to get things working properly and passing all tests. Once this is completed, I'll merge the pull request. I wasn't expecting the number of merge conflicts (apologies if something was unintentionally kicked out), but hopefully the weekly rebasing and additional communication will mitigate these issues! Also, if you develop a preference on the order of rebasing/committing/etc to make your life easier, please let me know!

[aman-raj-srivastva]

Hey @aman-raj-srivastva ! Please continue working in the current branch to get things working properly and passing all tests. Once this is completed, I'll merge the pull request. I wasn't expecting the number of merge conflicts (apologies if something was unintentionally kicked out), but hopefully the weekly rebasing and additional communication will mitigate these issues! Also, if you develop a preference on the order of rebasing/committing/etc to make your life easier, please let me know!

Thanks, Molly! That makes sense. I'll continue working on the current branch, resolve the conflicts from the rebase, make sure everything is passing locally, and then push the updated changes for review.

[mollystevens-noaa]

Hey @aman-raj-srivastva ! I met with @kellijohnson-NOAA and @Andrea-Havron-NOAA this afternoon, and we have an updated plan forward that will require much less overhead. The rebase to main conducted here will be the only full rebase we'll do until the end of the project, where we'll leave a few days for you to handle any additional merge conflicts / breaks that occur. For the remainder of the project, you can just work off this main-edm branch, and @Andrea-Havron-NOAA will let us know if we need to cherry pick any changes from the main branch that will affect your work. Bringing in @nathanvaughan-NOAA for awareness. Hope your resolution of conflicts here is going well, and please let us know if you have any questions or concerns!

[mollystevens-noaa]

Thanks @aman-raj-srivastva ! I see this version is passing all checks, and looks good to me. Before I merge this pull request, I wanted to circle back with @nathanvaughan-NOAA and @Andrea-Havron-NOAA to make sure that it fits as expected within the FIMS framework.

[nathanvaughan-NOAA]

Thanks @aman-raj-srivastva for all the great work and @mollystevens-noaa for testing out the new functions. I had a discussion with @Andrea-Havron-NOAA who will have a more detailed comment but we are concerned that at the moment the EDM module is still isolated from the rest of the FIMS codebase. Specifically there are a lot of calculations added inside the GP-EDM for estimating scale that should be a part of the FIMS distributions module. In it's final state the EDM module should only be creating the delay embedding map, calculating distance weights for each of the embedded vectors, and calculating predicted values from the embedding map and weights which are all deterministic calculations. The multivariate normal embedding function should be a part of the distributions functor so it can be used in EDM and for other relationships in FIMS, and the scale parameters should be specified as fixed or estimable similarly to all the other parameters in FIMS. @e-perl-NOAA is working on a similar issue in pull request #1594 to add a GMRF likelihood function.

[Andrea-Havron-NOAA]

It is more efficient to calculate this using the Cholesky decomposition than Gaussian Elimination, which can be evaluated using Eigen's ldlt, for example:
Given a matrix Sigma(n,n):

// Perform LLT (Cholesky) decomposition using Eigen via TMB
Eigen::LDLT<Eigen::Matrix<scalartype,Dynamic,Dynamic> > ldlt(Sigma);

// 1. Solve Sigma * alpha = y  -->  L L^T alpha = y
vector<Type> alpha = ldlt.solve(y);

For reference, see TMB's MVNORM function here

[Andrea-Havron-NOAA]

@msupernaw, should this use CppAD::sign (uses an atomic) instead so it is differentiable?

[Andrea-Havron-NOAA]

It's great to the see the multivariate normal function in FIMS! I would prefer you put it, however, in the distribution functors so that this likelihood can be used by the entire code base. You can use TMB's MVNORM function. If you use an explicit definition, update the inverse calculation of Sigma to use the cholesky decomposition instead of the Gaussian Elimination as mentioned above.

[Andrea-Havron-NOAA]

use fims_math::exp() throughout

[Andrea-Havron-NOAA]

see TMB example for efficiently calculating the log determinant from a cholesky decomposition

[Andrea-Havron-NOAA]

as noted above for the multivariate normal distribution, please add the beta distribution to the distribution functor folder and expose using the FIMS distribution pathway so these distribtions are available throughout the codebase.

[Andrea-Havron-NOAA]

it is safe to use a sqrt() inside the AD tape, we use fims_math::sqrt(). Generally, when taking the sqrt of a parameter, we estimate the parameter in log space and take the exponential to keep the value positive. If you can't use this trick, however, when evaluating the s-map projections, please feel free to leave the kGaussian approach as is.

[Andrea-Havron-NOAA]

use cholesky factorization as noted in the gp_edm_projection file.

[Andrea-Havron-NOAA]

Thanks for this PR and all your hard work! As @nathanvaughan-NOAA mentioned above, I think there are some changes that can be made to better integrate this code into the larger FIMS code base. Specifically,

1. distributions should all be defined in the inst/include/distributions/functors folder following the existing pattern for FIMS distributions. It is ok to use TMB macros here as long as they are wrapped in a #ifdef TMB_MODEL/#endif statement
2. go through the codebase and update basic functions to pull from fims_math instead of std::. You also defined a couple of functions that are already defined in fims::math. This just helps cut down redundancy and helps keep all our basic functions in one place so they are easy to access and test.

I noticed in your next PR, you have more code that helps integrate the edm piece into FIMS, I am looking forward to checking it out!