[Minuit2] Analytical Hessian indexing mismatch

[guitargeek]

When using analytic external Hessians in Minuit 2 (possible since 888a767), the index mapping between Minuit-internal floating parameters and all the "external" parameters of the function is not applied correctly.

Reproducer:

#include <Minuit2/FCNBase.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnPrint.h>

#include <cassert>
#include <vector>
#include <cmath>
#include <iostream>

using namespace ROOT::Minuit2;

class QuadraticFCN : public FCNBase {
public:
   double operator()(const std::vector<double>& p) const override {
      const double x = p[0];
      const double y = p[1];
      const double z = p[2];
      return x*x + 10*y*y + 100*z*z + 2*x*y + 4*x*z + 8*y*z;
   }

   bool HasGradient() const override { return true; }
   bool HasHessian()  const override { return true; }
   bool HasG2()       const override { return false; }

   std::vector<double> Gradient(const std::vector<double>& p) const override {
      return {
         2*p[0] + 2*p[1] + 4*p[2],
         20*p[1] + 2*p[0] + 8*p[2],
         200*p[2] + 4*p[0] + 8*p[1]
      };
   }

   std::vector<double> Hessian(const std::vector<double>&) const override {
      // Row-major 3x3
      return {
          2,  2,   4,
          2, 20,   8,
          4,  8, 200
      };
   }

   double Up() const override { return 1.0; }
};

void test_Hessian_external_indexing_bug()
{
   QuadraticFCN fcn;

   MnUserParameters upar;
   upar.Add("x", 1.0, 0.1);
   upar.Add("y", 2.0, 0.1);
   upar.Add("z", 3.0, 0.1);

   // FIX x to break the internal/external identity mapping
   upar.Fix("x");

   MnMigrad migrad(fcn, upar);
   FunctionMinimum min = migrad();

   const auto& h = min.Error().Hessian();

   // Internal parameters are (y, z)
   // Expected Hessian:
   // [20  8]
   // [ 8 200]
   //
   // Packed: [20, 8, 200]

   std::cout << "Internal Hessian:\n";
   std::cout << "h(0,0) = " << h(0,0) << "\n";
   std::cout << "h(0,1) = " << h(0,1) << "\n";
   std::cout << "h(1,1) = " << h(1,1) << "\n";

   std::cout << "Expected Hessian:\n";
   std::cout << "h(0,0) = 20\n";
   std::cout << "h(0,1) = 8\n";
   std::cout << "h(1,1) = 200\n";
}

Output:

Internal Hessian:
h(0,0) = 2
h(0,1) = 2
h(1,1) = 20
Expected Hessian:
h(0,0) = 20
h(0,1) = 8
h(1,1) = 200