remove remaining eigen to gsl conversions #40331
Conversation
MialLewis
added
the
ISIS: Core
sf1919
moved this from In Progress
to Waiting for Review
in ISIS core workstream v6.15.0
MialLewis
moved this from Waiting for Review
to In Progress
in ISIS core workstream v6.15.0
rboston628
left a comment
rboston628
left a comment
There was a problem hiding this comment.
This simplifies a lot. I would just like better docs for the new method.
| // Replicates equivalent gsl function | ||
| Eigen::MatrixXd covar_from_jacobian(const map_type &J, double epsrel, const map_type &r) { | ||
| const Eigen::Index m = J.rows(), n = J.cols(); | ||
| const int dof = std::max<int>(1, static_cast<int>(m - n)); // avoid div by 0 for safety | ||
| const double s2 = r.squaredNorm() / dof; // residual variance | ||
|
|
||
| Eigen::JacobiSVD<Eigen::MatrixXd> svd(J, Eigen::ComputeThinU | Eigen::ComputeThinV); | ||
| const Eigen::VectorXd &s = svd.singularValues(); // length = min(m,n) | ||
| const Eigen::MatrixXd &V = svd.matrixV(); // n x n (thin: n x r) | ||
|
|
||
| // GSL-style relative threshold on singular values | ||
| const double smax = s.size() ? s.array().maxCoeff() : 0.0; | ||
| const double tol = epsrel * smax; | ||
|
|
||
| // Build diag(1/s_i^2) with cutoff | ||
| Eigen::VectorXd inv_s2 = s.unaryExpr([&](double si) { return (si > tol) ? 1.0 / (si * si) : 0.0; }); | ||
|
|
||
| // Covariance = s2 * V * diag(1/s_i^2) * V^T | ||
| Eigen::MatrixXd C = V * inv_s2.asDiagonal() * V.transpose(); | ||
| C *= s2; | ||
| return C; | ||
| } |
There was a problem hiding this comment.
In the GSL documentation there is no entry for gsl_multifit_covar. I did find some some documentation here, but this is a much simpler method than you have. It seems like you're following steps in GSL's overview of linear least-squares fitting.
Could you add a link to that section, and then add clarifying comments for s, s2, smax, V, etc., and also how they relate tot he math operations there?
There was a problem hiding this comment.
Thanks for this, I've moved this to draft as I wasn't able to get this working 100% as gsl_multifit_covar did. I'm planning to pick this back up later in the sprint once our "Must" issues are out of the way.
I was adding tests as I realised this!
There was a problem hiding this comment.
Hi @rboston628, I've reworked the new function to match up with the docs (with the help of AI, as this is a bit beyond my limited liner algebra ability). I've gone through it and added comments to document; I think it's correct.
If you could take another look it'd be much appreciated.
There was a problem hiding this comment.
The difference between reference files are very minimal as per the figure above; I'm going to adjust the results for this one.
MialLewis
force-pushed
the
0_remove_eigen_to_gsl_conversions
branch
from
25a5191 to
f9d4c88
Compare
Description of work
There were several points in the code where EigenMatrix/EigenVector objects were converted into gsl objects to use in gsl functions. This PR removes these cases, swapping out the gsl functions for eigen implementations.
Have enabled the system test which ref file I changed on windows and Mac, since the tests pass locally for me on Mac.
Here's a windows test: https://builds.mantidproject.org/job/build_branch/1446/
To test:
Reviewer
Your comments will be used as part of the gatekeeper process. Comment clearly on what you have checked and tested during your review. Provide an audit trail for any changes requested.
As per the review guidelines:
mantid-developersormantid-contributorsteams, add a review commentrerun cito authorize/rerun the CIGatekeeper
As per the gatekeeping guidelines: