AlignAndFocusPowderSlim for NOMAD and POWGEN (#40664)

This pulls the following into `ornl-next`
* #40642
* #40644 
* #40645 
* #40643 
* #40657
* #40331 
* #40633

---------

Co-authored-by: mantid-builder <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Ross Whitfield <[email protected]>
Co-authored-by: MialLewis <[email protected]>

Author: 5 people

.pre-commit-config.yaml

@@ -68,7 +68,7 @@ repos:
           )$
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.14.11
+    rev: v0.14.13
     # ruff must appear before black in the list of hooks
     hooks:
       - id: ruff-check

Framework/CurveFitting/CMakeLists.txt

@@ -41,6 +41,7 @@ set(SRC_FILES
     src/FitMW.cpp
     src/EigenComplexMatrix.cpp
     src/EigenComplexVector.cpp
+    src/EigenFunctions.cpp
     src/EigenMatrix.cpp
     src/EigenMatrixView.cpp
     src/EigenVector.cpp
@@ -347,6 +348,7 @@ set(INC_FILES
     inc/MantidCurveFitting/Functions/VesuvioResolution.h
     inc/MantidCurveFitting/Functions/Voigt.h
     inc/MantidCurveFitting/Functions/ConvTempCorrection.h
+    inc/MantidCurveFitting/EigenFunctions.h
     inc/MantidCurveFitting/GSLFunctions.h
     inc/MantidCurveFitting/GeneralDomainCreator.h
     inc/MantidCurveFitting/HalfComplex.h
@@ -409,9 +411,9 @@ set(TEST_FILES
     EigenComplexVectorTest.h
     EigenFortranMatrixTest.h
     EigenFortranVectorTest.h
+    EigenFunctionsTest.h
     EigenJacobianTest.h
     EigenMatrixTest.h
-    EigenToGSLTest.h
     EigenVectorTest.h
     EigenViewTest.h
     FitMWTest.h

Framework/CurveFitting/inc/MantidCurveFitting/EigenFunctions.h

@@ -0,0 +1,15 @@
+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2026 ISIS Rutherford Appleton Laboratory UKRI,
+//   NScD Oak Ridge National Laboratory, European Spallation Source,
+//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+#pragma once
+
+#include "MantidCurveFitting/EigenMatrix.h"
+
+namespace Mantid::CurveFitting {
+
+MANTID_CURVEFITTING_DLL Eigen::MatrixXd covar_from_jacobian(const map_type &J, double epsrel);
+
+} // namespace Mantid::CurveFitting

Framework/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h

@@ -50,19 +50,5 @@ int gsl_f(const gsl_vector *x, void *params, gsl_vector *f);
 int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J);
 int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J);
 
-/// take data from Eigen Vector and take a gsl view
-inline gsl_vector_view getGSLVectorView(vec_map_type &v) { return gsl_vector_view_array(v.data(), v.size()); }
-/// take data from an Eigen Matrix and return a transposed a gsl view.
-inline gsl_matrix_view getGSLMatrixView(map_type &tr) { return gsl_matrix_view_array(tr.data(), tr.cols(), tr.rows()); }
-
-/// take const data from Eigen Vector and take a gsl view
-inline gsl_vector_const_view const getGSLVectorView_const(const vec_map_type v) {
-  return gsl_vector_const_view_array(v.data(), v.size());
-}
-/// take data from a constEigen Matrix and return a transposed gsl view.
-inline gsl_matrix_const_view const getGSLMatrixView_const(const map_type m) {
-  return gsl_matrix_const_view_array(m.data(), m.cols(), m.rows());
-}
-
 } // namespace CurveFitting
 } // namespace Mantid

Framework/CurveFitting/src/CostFunctions/CostFuncFitting.cpp

@@ -9,12 +9,11 @@
 //----------------------------------------------------------------------
 #include "MantidCurveFitting/CostFunctions/CostFuncFitting.h"
 #include "MantidAPI/IConstraint.h"
+#include "MantidCurveFitting/EigenFunctions.h"
 #include "MantidCurveFitting/EigenJacobian.h"
-#include "MantidCurveFitting/GSLFunctions.h"
 #include "MantidCurveFitting/SeqDomain.h"
 #include "MantidKernel/Exception.h"
 
-#include <gsl/gsl_multifit_nlin.h>
 #include <limits>
 #include <utility>
 
@@ -132,15 +131,8 @@ void CostFuncFitting::calActiveCovarianceMatrix(EigenMatrix &covar, double epsre
   // calculate the derivatives
   m_function->functionDeriv(*m_domain, J);
 
-  // let the GSL to compute the covariance matrix
-  EigenMatrix covarTr(covar.tr());
-  EigenMatrix tempJTr;
-  tempJTr = j.transpose();
-  gsl_matrix_const_view tempJTrGSL = getGSLMatrixView_const(tempJTr.inspector());
-  gsl_matrix_view covarTrGSL = getGSLMatrixView(covarTr.mutator());
-
-  gsl_multifit_covar(&tempJTrGSL.matrix, epsrel, &covarTrGSL.matrix);
-  covar = covarTr.tr();
+  // let Eigen compute the covariance matrix
+  covar = covar_from_jacobian(j, epsrel);
 }
 
 /** Calculates covariance matrix

Framework/CurveFitting/src/EigenFunctions.cpp

@@ -0,0 +1,70 @@
+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2026 ISIS Rutherford Appleton Laboratory UKRI,
+//   NScD Oak Ridge National Laboratory, European Spallation Source,
+//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidCurveFitting/EigenFunctions.h"
+
+namespace Mantid::CurveFitting {
+
+/** Mimics gsl_multifit_covar(J, epsrel, covar)
+ * @param J :: Jacobian
+ * @param epsrel :: relative threshold to decide rank deficiency
+ * @return covariance matrix, rows/cols set to zero for dependent params
+ */
+Eigen::MatrixXd covar_from_jacobian(const map_type &J, double epsrel) {
+  if (epsrel < 0.0) {
+    throw std::invalid_argument("epsrel must be non-negative");
+  }
+
+  const Eigen::Index nr = J.rows(); // nr = num rows
+  const Eigen::Index nc = J.cols(); // nc = num cols
+  if ((nc == 0) || (nr == 0))
+    return Eigen::MatrixXd{};
+
+  // Pivoted QR Decomposition: JP = QR
+  Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(J);
+
+  // R is (nc x np) in general; the useful part for J^T J is (nc x nc) block
+  Eigen::MatrixXd R = qr.matrixR().topLeftCorner(nc, nc).template triangularView<Eigen::Upper>().toDenseMatrix();
+
+  // Determine rank using the same rule as GSL: compare diag entries to |R_11|
+  const double r11 = std::abs(R(0, 0));
+  Eigen::Index rank = 0;
+  if (r11 > 0.0) {
+    for (Eigen::Index k = 0; k < nc; ++k) {
+      if (std::abs(R(k, k)) > epsrel * r11) { // column considered linerally independant
+        ++rank;
+      } else {
+        break; // with pivoting, following first dependent subsequent cols should also be dependant
+      }
+    }
+  } else {
+    // R11 == 0, everything dependent
+    rank = 0;
+  }
+
+  // Build covariance in the pivoted parameter order
+  Eigen::MatrixXd cov_pivot = Eigen::MatrixXd::Zero(nc, nc);
+
+  if (rank > 0) {
+    // cov = (R^T R)^{-1} for the independent cols = R^{-1} R^{-T}
+    const Eigen::MatrixXd R1 = R.topLeftCorner(rank, rank).template triangularView<Eigen::Upper>();
+    Eigen::MatrixXd invR1 = R1.inverse();
+    Eigen::MatrixXd cov1 = invR1 * invR1.transpose();
+    cov_pivot.topLeftCorner(rank, rank) = cov1;
+  }
+
+  // Unpivot back to original parameter order:
+  // cov = (J^T J)^{-1} = P (R^T R)^{-1} P^T
+  const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> P = qr.colsPermutation();
+  Eigen::MatrixXd cov = P * cov_pivot * P.transpose();
+
+  return cov;
+}
+
+} // namespace Mantid::CurveFitting

Framework/CurveFitting/src/FuncMinimizers/DampedGaussNewtonMinimizer.cpp

@@ -18,7 +18,6 @@
 #include "MantidKernel/Logger.h"
 
 #include <cmath>
-#include <gsl/gsl_blas.h>
 
 namespace Mantid::CurveFitting::FuncMinimisers {
 
@@ -121,9 +120,7 @@ bool DampedGaussNewtonMinimizer::iterate(size_t /*iteration*/) {
   // Try the stop condition
   EigenVector p(n);
   m_leastSquares->getParameters(p);
-  gsl_vector_view dx_gsl = getGSLVectorView(dx.mutator());
-  double dx_norm = gsl_blas_dnrm2(&dx_gsl.vector);
-  return dx_norm >= m_relTol;
+  return dx.norm() >= m_relTol;
 }
 
 /// Return current value of the cost function

Framework/CurveFitting/src/FuncMinimizers/LevenbergMarquardtMDMinimizer.cpp

@@ -18,7 +18,6 @@
 #include "MantidKernel/Logger.h"
 
 #include <cmath>
-#include <gsl/gsl_blas.h>
 
 namespace Mantid::CurveFitting::FuncMinimisers {
 namespace {
@@ -194,15 +193,11 @@ bool LevenbergMarquardtMDMinimizer::iterate(size_t /*iteration*/) {
   double dL;
   // der -> - der - 0.5 * hessian * dx
 
-  EigenMatrix tempHessianTr = m_costFunction->getHessian().tr();
-  const gsl_matrix_const_view tempHessianTrGSL = getGSLMatrixView_const(tempHessianTr.inspector());
-  const gsl_vector_const_view dxGSL = getGSLVectorView_const(dx.inspector());
-  gsl_vector_view ddGSL = getGSLVectorView(dd.mutator());
-
-  gsl_blas_dgemv(CblasNoTrans, -0.5, &tempHessianTrGSL.matrix, &dxGSL.vector, 1., &ddGSL.vector);
+  EigenVector Trdx = m_costFunction->getHessian().tr() * dx;
+  Trdx *= -0.5;
+  dd += Trdx;
   // calculate the linear part of the change in cost function
-  // dL = - der * dx - 0.5 * dx * hessian * dx
-  gsl_blas_ddot(&ddGSL.vector, &dxGSL.vector, &dL);
+  dL = dd.dot(dx);
 
   double F1 = m_costFunction->val();
   if (verbose) {
@@ -216,7 +211,7 @@ bool LevenbergMarquardtMDMinimizer::iterate(size_t /*iteration*/) {
   if (m_rho >= 0) {
     EigenVector p(n);
     m_costFunction->getParameters(p);
-    double dx_norm = gsl_blas_dnrm2(&dxGSL.vector);
+    double dx_norm = dx.norm();
     if (dx_norm < absError) {
       if (verbose) {
         g_log.warning() << "Successful fit, parameters changed by less than " << absError << '\n';

Framework/CurveFitting/src/Functions/ChebfunBase.cpp

@@ -11,11 +11,11 @@
 #include "MantidCurveFitting/GSLFunctions.h"
 #include "MantidCurveFitting/HalfComplex.h"
 
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_errno.h>
 #include <gsl/gsl_fft_halfcomplex.h>
 #include <gsl/gsl_fft_real.h>
 
+#include <Eigen/Eigenvalues>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -684,21 +684,16 @@ std::vector<double> ChebfunBase::roots(const std::vector<double> &a) const {
     C.set(N + i, lasti, tmp);
   }
 
-  gsl_vector_complex *eigenval = gsl_vector_complex_alloc(N2);
-  auto workspace = gsl_eigen_nonsymm_alloc(N2);
-
-  EigenMatrix C_tr(C.tr());
-  gsl_matrix_view C_tr_gsl = getGSLMatrixView(C_tr.mutator());
-  gsl_eigen_nonsymm(&C_tr_gsl.matrix, eigenval, workspace);
-  gsl_eigen_nonsymm_free(workspace);
+  Eigen::EigenSolver<Eigen::MatrixXd> es(C.mutator());
+  const Eigen::VectorXcd evals = es.eigenvalues();
 
   const double Dx = endX() - startX();
   bool isFirst = true;
   double firstIm = 0;
   for (size_t i = 0; i < N2; ++i) {
-    auto val = gsl_vector_complex_get(eigenval, i);
-    double re = GSL_REAL(val);
-    double im = GSL_IMAG(val);
+    const std::complex<double> val = evals[static_cast<Eigen::Index>(i)];
+    double re = val.real();
+    double im = val.imag();
     double ab = re * re + im * im;
     if (fabs(ab - 1.0) > 1e-2) {
       isFirst = true;
@@ -715,8 +710,6 @@ std::vector<double> ChebfunBase::roots(const std::vector<double> &a) const {
       isFirst = true;
     }
   }
-  gsl_vector_complex_free(eigenval);
-
   return r;
 }
 
@@ -824,25 +817,24 @@ std::vector<double> ChebfunBase::smooth(const std::vector<double> &xvalues, cons
       return y;
     }
 
-    // high frequency filter values: smooth decreasing function
-    auto ri0f = static_cast<double>(i0 + 1);
-    double xm = (1.0 + ri0f) / 2;
-    ym /= ri0f;
-    double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
-    double b1 = ym - a1 * xm;
-
-    // std::cerr << "(a1,b1) = (" << a1 << ',' << b1 << ')' << '\n';
-
     // calculate coeffs of a cubic c3*i^3 + c2*i^2 + c1*i + c0
     // which will replace the linear a1*i + b1 in building the
     // second part of the filter
     double c0, c1, c2, c3;
     {
+      // high frequency filter values: smooth decreasing function
+      auto ri0f = static_cast<double>(i0 + 1);
+      double xm = (1.0 + ri0f) / 2;
+      ym /= ri0f;
+
       auto x0 = double(i0 + 1);
       auto x1 = double(n + 1);
       double sigma = g_tolerance / noise / 10;
       double s = sigma / (1.0 - sigma);
       double m1 = log(s);
+      double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
+      double b1 = ym - a1 * xm;
+      // std::cerr << "(a1,b1) = (" << a1 << ',' << b1 << ')' << '\n';
       double m0 = a1 * x0 + b1;
       if (a1 < 0.0) {
         c3 = 0.0;

Framework/CurveFitting/test/CostFunctions/LeastSquaresTest.h

@@ -22,7 +22,6 @@
 #include "MantidCurveFitting/Functions/UserFunction.h"
 #include "MantidCurveFitting/GSLFunctions.h"
 
-#include <gsl/gsl_blas.h>
 #include <sstream>
 
 using namespace Mantid;
@@ -201,15 +200,11 @@ class LeastSquaresTest : public CxxTest::TestSuite {
     dx.set(1, -0.2);
 
     double L; // = d*dx + 0.5 * dx * H * dx
+    EigenVector Trdx = H.tr() * dx;
+    Trdx *= 0.5;
+    g += Trdx;
+    L = g.dot(dx);
 
-    EigenMatrix temp_H_tr = H.tr();
-    const gsl_matrix_const_view temp_H_tr_gsl = getGSLMatrixView_const(temp_H_tr.inspector());
-    const gsl_vector_const_view dx_gsl = getGSLVectorView_const(dx.inspector());
-    gsl_vector_view g_gsl = getGSLVectorView(g.mutator());
-
-    gsl_blas_dgemv(CblasNoTrans, 0.5, &temp_H_tr_gsl.matrix, &dx_gsl.vector, 1., &g_gsl.vector);
-
-    gsl_blas_ddot(&g_gsl.vector, &dx_gsl.vector, &L);
     TS_ASSERT_DELTA(L, -0.145, 1e-10); // L + costFun->val() == 0
   }