OCCT

// Copyright (c) 2025 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <gtest/gtest.h>

// New API
#include <MathOpt_Powell.hxx>
#include <MathOpt_BFGS.hxx>

// Old API
#include <math_Powell.hxx>
#include <math_BFGS.hxx>
#include <math_MultipleVarFunction.hxx>
#include <math_MultipleVarFunctionWithGradient.hxx>
#include <math_Matrix.hxx>

#include <cmath>

namespace
{
constexpr double THE_TOLERANCE = 1.0e-6;

// ============================================================================
// Adapter classes for old API (inherit from math_MultipleVarFunction)
// ============================================================================

//! Quadratic function: f(x,y) = (x-1)^2 + (y-2)^2
//! Minimum at (1, 2) with f = 0
class QuadraticFuncOld : public math_MultipleVarFunction
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aDx = theX.At(0) - 1.0;
    const double aDy = theX.At(1) - 2.0;
    theF             = aDx * aDx + aDy * aDy;
    return true;
  }
};

//! Rosenbrock function: f(x,y) = 100*(y-x^2)^2 + (1-x)^2
//! Minimum at (1, 1) with f = 0
class RosenbrockFuncOld : public math_MultipleVarFunction
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aX  = theX(1);
    const double aY  = theX(2);
    const double aT1 = aY - aX * aX;
    const double aT2 = 1.0 - aX;
    theF             = 100.0 * aT1 * aT1 + aT2 * aT2;
    return true;
  }
};

//! Booth function: f(x,y) = (x + 2y - 7)^2 + (2x + y - 5)^2
//! Minimum at (1, 3) with f = 0
class BoothFuncOld : public math_MultipleVarFunction
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aX  = theX(1);
    const double aY  = theX(2);
    const double aT1 = aX + 2.0 * aY - 7.0;
    const double aT2 = 2.0 * aX + aY - 5.0;
    theF             = aT1 * aT1 + aT2 * aT2;
    return true;
  }
};

//! Sphere function in N dimensions: f(x) = sum(x_i^2)
//! Minimum at origin with f = 0
class SphereFuncOld : public math_MultipleVarFunction
{
private:
  int myN;

public:
  SphereFuncOld(int theN)
      : myN(theN)
  {
  }

  int NbVariables() const override { return myN; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    theF = 0.0;
    for (int i = theX.Lower(); i <= theX.Upper(); ++i)
    {
      theF += theX(i) * theX(i);
    }
    return true;
  }
};

// ============================================================================
// Adapter classes for old BFGS API (inherit from math_MultipleVarFunctionWithGradient)
// ============================================================================

//! Quadratic function with gradient: f(x,y) = (x-1)^2 + (y-2)^2
class QuadraticFuncGradOld : public math_MultipleVarFunctionWithGradient
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aDx = theX(1) - 1.0;
    const double aDy = theX(2) - 2.0;
    theF             = aDx * aDx + aDy * aDy;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theG) override
  {
    theG(1) = 2.0 * (theX(1) - 1.0);
    theG(2) = 2.0 * (theX(2) - 2.0);
    return true;
  }

  bool Values(const math_Vector& theX, double& theF, math_Vector& theG) override
  {
    return Value(theX, theF) && Gradient(theX, theG);
  }
};

//! Rosenbrock function with gradient
class RosenbrockFuncGradOld : public math_MultipleVarFunctionWithGradient
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aX  = theX(1);
    const double aY  = theX(2);
    const double aT1 = aY - aX * aX;
    const double aT2 = 1.0 - aX;
    theF             = 100.0 * aT1 * aT1 + aT2 * aT2;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theG) override
  {
    const double aX = theX(1);
    const double aY = theX(2);
    theG(1)         = -400.0 * aX * (aY - aX * aX) - 2.0 * (1.0 - aX);
    theG(2)         = 200.0 * (aY - aX * aX);
    return true;
  }

  bool Values(const math_Vector& theX, double& theF, math_Vector& theG) override
  {
    return Value(theX, theF) && Gradient(theX, theG);
  }
};

//! Booth function with gradient
class BoothFuncGradOld : public math_MultipleVarFunctionWithGradient
{
public:
  int NbVariables() const override { return 2; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    const double aX  = theX(1);
    const double aY  = theX(2);
    const double aT1 = aX + 2.0 * aY - 7.0;
    const double aT2 = 2.0 * aX + aY - 5.0;
    theF             = aT1 * aT1 + aT2 * aT2;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theG) override
  {
    const double aX  = theX(1);
    const double aY  = theX(2);
    const double aT1 = aX + 2.0 * aY - 7.0;
    const double aT2 = 2.0 * aX + aY - 5.0;
    theG(1)          = 2.0 * aT1 + 4.0 * aT2;
    theG(2)          = 4.0 * aT1 + 2.0 * aT2;
    return true;
  }

  bool Values(const math_Vector& theX, double& theF, math_Vector& theG) override
  {
    return Value(theX, theF) && Gradient(theX, theG);
  }
};

//! Sphere function with gradient
class SphereFuncGradOld : public math_MultipleVarFunctionWithGradient
{
private:
  int myN;

public:
  SphereFuncGradOld(int theN)
      : myN(theN)
  {
  }

  int NbVariables() const override { return myN; }

  bool Value(const math_Vector& theX, double& theF) override
  {
    theF = 0.0;
    for (int i = theX.Lower(); i <= theX.Upper(); ++i)
    {
      theF += theX(i) * theX(i);
    }
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theG) override
  {
    for (int i = theX.Lower(); i <= theX.Upper(); ++i)
    {
      theG(i) = 2.0 * theX(i);
    }
    return true;
  }

  bool Values(const math_Vector& theX, double& theF, math_Vector& theG) override
  {
    return Value(theX, theF) && Gradient(theX, theG);
  }
};

// ============================================================================
// Function classes for new API
// ============================================================================

struct QuadraticFuncNew
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aDx = theX.At(0) - 1.0;
    const double aDy = theX.At(1) - 2.0;
    theF             = aDx * aDx + aDy * aDy;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theGrad)
  {
    theGrad.ChangeAt(0) = 2.0 * (theX.At(0) - 1.0);
    theGrad.ChangeAt(1) = 2.0 * (theX.At(1) - 2.0);
    return true;
  }
};

struct RosenbrockFuncNew
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aX  = theX.At(0);
    const double aY  = theX.At(1);
    const double aT1 = aY - aX * aX;
    const double aT2 = 1.0 - aX;
    theF             = 100.0 * aT1 * aT1 + aT2 * aT2;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theGrad)
  {
    const double aX     = theX.At(0);
    const double aY     = theX.At(1);
    theGrad.ChangeAt(0) = -400.0 * aX * (aY - aX * aX) - 2.0 * (1.0 - aX);
    theGrad.ChangeAt(1) = 200.0 * (aY - aX * aX);
    return true;
  }
};

struct BoothFuncNew
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aX  = theX.At(0);
    const double aY  = theX.At(1);
    const double aT1 = aX + 2.0 * aY - 7.0;
    const double aT2 = 2.0 * aX + aY - 5.0;
    theF             = aT1 * aT1 + aT2 * aT2;
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theGrad)
  {
    const double aX     = theX.At(0);
    const double aY     = theX.At(1);
    const double aT1    = aX + 2.0 * aY - 7.0;
    const double aT2    = 2.0 * aX + aY - 5.0;
    theGrad.ChangeAt(0) = 2.0 * aT1 + 4.0 * aT2;
    theGrad.ChangeAt(1) = 4.0 * aT1 + 2.0 * aT2;
    return true;
  }
};

struct SphereFuncNew
{
  bool Value(const math_Vector& theX, double& theF)
  {
    theF = 0.0;
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theF += theX.At(i) * theX.At(i);
    }
    return true;
  }

  bool Gradient(const math_Vector& theX, math_Vector& theGrad)
  {
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theGrad.ChangeAt(i) = 2.0 * theX.At(i);
    }
    return true;
  }
};

//! Create identity matrix for Powell directions.
math_Matrix CreateIdentityMatrix(int theN)
{
  math_Matrix aMat(1, theN, 1, theN, 0.0);
  for (int i = 1; i <= theN; ++i)
  {
    aMat(i, i) = 1.0;
  }
  return aMat;
}

} // namespace

// ============================================================================
// math_Powell vs MathOpt::Powell comparison tests
// ============================================================================

TEST(MathOpt_NDim_ComparisonTest, Powell_Quadratic)
{
  QuadraticFuncOld anOldFunc;
  QuadraticFuncNew aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = 5.0;
  anOldStart(2) = 7.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = 5.0;
  aNewStart.ChangeAt(1) = 7.0;

  // Old API
  math_Powell anOldSolver(anOldFunc, 1.0e-8, 200);
  math_Matrix aDirections = CreateIdentityMatrix(2);
  anOldSolver.Perform(anOldFunc, anOldStart, aDirections);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 200;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::Powell(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);

  // Verify correctness
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aNewResult.Solution->At(1), 2.0, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_NDim_ComparisonTest, Powell_Booth)
{
  BoothFuncOld anOldFunc;
  BoothFuncNew aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = 0.0;
  anOldStart(2) = 0.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = 0.0;
  aNewStart.ChangeAt(1) = 0.0;

  // Old API
  math_Powell anOldSolver(anOldFunc, 1.0e-8, 200);
  math_Matrix aDirections = CreateIdentityMatrix(2);
  anOldSolver.Perform(anOldFunc, anOldStart, aDirections);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 200;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::Powell(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);

  // Verify correctness
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aNewResult.Solution->At(1), 3.0, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_NDim_ComparisonTest, Powell_Rosenbrock)
{
  RosenbrockFuncOld anOldFunc;
  RosenbrockFuncNew aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = -1.0;
  anOldStart(2) = 1.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = -1.0;
  aNewStart.ChangeAt(1) = 1.0;

  // Old API - needs more iterations for Rosenbrock
  math_Powell anOldSolver(anOldFunc, 1.0e-6, 1000);
  math_Matrix aDirections = CreateIdentityMatrix(2);
  anOldSolver.Perform(anOldFunc, anOldStart, aDirections);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 1000;
  aConfig.XTolerance    = 1.0e-6;
  aConfig.FTolerance    = 1.0e-8;
  auto aNewResult       = MathOpt::Powell(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the minimum (relaxed tolerance for Rosenbrock)
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), 1.0e-3);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), 1.0e-3);

  // Verify correctness (Rosenbrock minimum at (1, 1))
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, 1.0e-3);
  EXPECT_NEAR(aNewResult.Solution->At(1), 1.0, 1.0e-3);
}

TEST(MathOpt_NDim_ComparisonTest, Powell_Sphere3D)
{
  SphereFuncOld anOldFunc(3);
  SphereFuncNew aNewFunc;

  math_Vector anOldStart(1, 3);
  anOldStart(1) = 3.0;
  anOldStart(2) = -2.0;
  anOldStart(3) = 4.0;

  math_Vector aNewStart(3);
  aNewStart.ChangeAt(0) = 3.0;
  aNewStart.ChangeAt(1) = -2.0;
  aNewStart.ChangeAt(2) = 4.0;

  // Old API
  math_Powell anOldSolver(anOldFunc, 1.0e-8, 200);
  math_Matrix aDirections = CreateIdentityMatrix(3);
  anOldSolver.Perform(anOldFunc, anOldStart, aDirections);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 200;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::Powell(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the origin
  for (size_t i = 0; i < aNewResult.Solution->Size(); ++i)
  {
    EXPECT_NEAR(anOldSolver.Location().At(i), aNewResult.Solution->At(i), THE_TOLERANCE);
    EXPECT_NEAR(aNewResult.Solution->At(i), 0.0, THE_TOLERANCE);
  }
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

// ============================================================================
// math_BFGS vs MathOpt::BFGS comparison tests
// ============================================================================

TEST(MathOpt_NDim_ComparisonTest, BFGS_Quadratic)
{
  QuadraticFuncGradOld anOldFunc;
  QuadraticFuncNew     aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = 5.0;
  anOldStart(2) = 7.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = 5.0;
  aNewStart.ChangeAt(1) = 7.0;

  // Old API
  math_BFGS anOldSolver(2, 1.0e-8, 100);
  anOldSolver.Perform(anOldFunc, anOldStart);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 100;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);

  // Verify correctness
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aNewResult.Solution->At(1), 2.0, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_NDim_ComparisonTest, BFGS_Booth)
{
  BoothFuncGradOld anOldFunc;
  BoothFuncNew     aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = 0.0;
  anOldStart(2) = 0.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = 0.0;
  aNewStart.ChangeAt(1) = 0.0;

  // Old API
  math_BFGS anOldSolver(2, 1.0e-8, 100);
  anOldSolver.Perform(anOldFunc, anOldStart);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 100;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), THE_TOLERANCE);
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);

  // Verify correctness
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aNewResult.Solution->At(1), 3.0, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_NDim_ComparisonTest, BFGS_Rosenbrock)
{
  RosenbrockFuncGradOld anOldFunc;
  RosenbrockFuncNew     aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = -1.0;
  anOldStart(2) = 1.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = -1.0;
  aNewStart.ChangeAt(1) = 1.0;

  // Old API - needs more iterations for Rosenbrock
  math_BFGS anOldSolver(2, 1.0e-6, 1000);
  anOldSolver.Perform(anOldFunc, anOldStart);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 1000;
  aConfig.XTolerance    = 1.0e-6;
  aConfig.FTolerance    = 1.0e-8;
  auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the minimum (relaxed tolerance for Rosenbrock)
  EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), 1.0e-3);
  EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), 1.0e-3);

  // Verify correctness (Rosenbrock minimum at (1, 1))
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, 1.0e-3);
  EXPECT_NEAR(aNewResult.Solution->At(1), 1.0, 1.0e-3);
}

TEST(MathOpt_NDim_ComparisonTest, BFGS_Sphere3D)
{
  SphereFuncGradOld anOldFunc(3);
  SphereFuncNew     aNewFunc;

  math_Vector anOldStart(1, 3);
  anOldStart(1) = 3.0;
  anOldStart(2) = -2.0;
  anOldStart(3) = 4.0;

  math_Vector aNewStart(3);
  aNewStart.ChangeAt(0) = 3.0;
  aNewStart.ChangeAt(1) = -2.0;
  aNewStart.ChangeAt(2) = 4.0;

  // Old API
  math_BFGS anOldSolver(3, 1.0e-8, 100);
  anOldSolver.Perform(anOldFunc, anOldStart);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 100;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;
  auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should find the origin
  for (size_t i = 0; i < aNewResult.Solution->Size(); ++i)
  {
    EXPECT_NEAR(anOldSolver.Location().At(i), aNewResult.Solution->At(i), THE_TOLERANCE);
    EXPECT_NEAR(aNewResult.Solution->At(i), 0.0, THE_TOLERANCE);
  }
  EXPECT_NEAR(anOldSolver.Minimum(), *aNewResult.Value, THE_TOLERANCE);
  EXPECT_NEAR(*aNewResult.Value, 0.0, THE_TOLERANCE);
}

// ============================================================================
// Iteration count comparison tests
// ============================================================================

TEST(MathOpt_NDim_ComparisonTest, IterationCount_BFGSVsPowell)
{
  QuadraticFuncNew aNewFunc;

  math_Vector aStart(2);
  aStart.ChangeAt(0) = 10.0;
  aStart.ChangeAt(1) = 10.0;

  MathOpt::Config aConfig;
  aConfig.MaxIterations = 500;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;

  auto aPowellResult = MathOpt::Powell(aNewFunc, aStart, aConfig);
  auto aBFGSResult   = MathOpt::BFGS(aNewFunc, aStart, aConfig);

  ASSERT_TRUE(aPowellResult.IsDone());
  ASSERT_TRUE(aBFGSResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(*aPowellResult.Value, *aBFGSResult.Value, THE_TOLERANCE);

  // BFGS should use fewer iterations on quadratic functions (due to superlinear convergence)
  EXPECT_LT(aBFGSResult.NbIterations, aPowellResult.NbIterations);
}

TEST(MathOpt_NDim_ComparisonTest, IterationCount_LBFGSVsBFGS)
{
  SphereFuncNew aNewFunc;

  math_Vector aStart(5);
  for (size_t i = 0; i < aStart.Size(); ++i)
  {
    aStart.ChangeAt(i) = static_cast<double>(i + 1);
  }

  MathOpt::Config aConfig;
  aConfig.MaxIterations = 100;
  aConfig.XTolerance    = 1.0e-8;
  aConfig.FTolerance    = 1.0e-10;

  auto aBFGSResult  = MathOpt::BFGS(aNewFunc, aStart, aConfig);
  auto aLBFGSResult = MathOpt::LBFGS(aNewFunc, aStart, 10, aConfig);

  ASSERT_TRUE(aBFGSResult.IsDone());
  ASSERT_TRUE(aLBFGSResult.IsDone());

  // Both should find the same minimum
  EXPECT_NEAR(*aBFGSResult.Value, *aLBFGSResult.Value, THE_TOLERANCE);

  // Both should find the origin
  for (size_t i = 0; i < aBFGSResult.Solution->Size(); ++i)
  {
    EXPECT_NEAR(aBFGSResult.Solution->At(i), 0.0, THE_TOLERANCE);
    EXPECT_NEAR(aLBFGSResult.Solution->At(i), 0.0, THE_TOLERANCE);
  }
}

// ============================================================================
// Accuracy comparison tests
// ============================================================================

TEST(MathOpt_NDim_ComparisonTest, Accuracy_BFGSVsOld)
{
  QuadraticFuncGradOld anOldFunc;
  QuadraticFuncNew     aNewFunc;

  math_Vector anOldStart(1, 2);
  anOldStart(1) = 5.0;
  anOldStart(2) = 7.0;

  math_Vector aNewStart(2);
  aNewStart.ChangeAt(0) = 5.0;
  aNewStart.ChangeAt(1) = 7.0;

  // Tighter tolerance
  const double aTightTol = 1.0e-12;

  // Old API
  math_BFGS anOldSolver(2, aTightTol, 200);
  anOldSolver.Perform(anOldFunc, anOldStart);

  // New API
  MathOpt::Config aConfig;
  aConfig.MaxIterations = 200;
  aConfig.XTolerance    = aTightTol;
  aConfig.FTolerance    = 1.0e-14;
  auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

  ASSERT_TRUE(anOldSolver.IsDone());
  ASSERT_TRUE(aNewResult.IsDone());

  // Both should achieve good accuracy
  EXPECT_NEAR(anOldSolver.Location()(1), 1.0, 1.0e-8);
  EXPECT_NEAR(anOldSolver.Location()(2), 2.0, 1.0e-8);
  EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, 1.0e-8);
  EXPECT_NEAR(aNewResult.Solution->At(1), 2.0, 1.0e-8);
}

// ============================================================================
// Different starting points comparison tests
// ============================================================================

TEST(MathOpt_NDim_ComparisonTest, DifferentStartingPoints)
{
  QuadraticFuncGradOld anOldFunc;
  QuadraticFuncNew     aNewFunc;

  // Test multiple starting points
  std::vector<std::pair<double, double>> aStartPoints = {{0.0, 0.0},
                                                         {10.0, 10.0},
                                                         {-5.0, 5.0},
                                                         {100.0, -50.0}};

  for (const auto& aPoint : aStartPoints)
  {
    math_Vector anOldStart(1, 2);
    anOldStart(1) = aPoint.first;
    anOldStart(2) = aPoint.second;

    math_Vector aNewStart(2);
    aNewStart.ChangeAt(0) = aPoint.first;
    aNewStart.ChangeAt(1) = aPoint.second;

    // Old API
    math_BFGS anOldSolver(2, 1.0e-8, 200);
    anOldSolver.Perform(anOldFunc, anOldStart);

    // New API
    MathOpt::Config aConfig;
    aConfig.MaxIterations = 200;
    aConfig.XTolerance    = 1.0e-8;
    aConfig.FTolerance    = 1.0e-10;
    auto aNewResult       = MathOpt::BFGS(aNewFunc, aNewStart, aConfig);

    ASSERT_TRUE(anOldSolver.IsDone())
      << "Old API failed for start point (" << aPoint.first << ", " << aPoint.second << ")";
    ASSERT_TRUE(aNewResult.IsDone())
      << "New API failed for start point (" << aPoint.first << ", " << aPoint.second << ")";

    // Both should find the same minimum
    EXPECT_NEAR(anOldSolver.Location().At(0), aNewResult.Solution->At(0), THE_TOLERANCE);
    EXPECT_NEAR(anOldSolver.Location().At(1), aNewResult.Solution->At(1), THE_TOLERANCE);
    EXPECT_NEAR(aNewResult.Solution->At(0), 1.0, THE_TOLERANCE);
    EXPECT_NEAR(aNewResult.Solution->At(1), 2.0, THE_TOLERANCE);
  }
}