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>

#include <MathOpt_GlobOpt.hxx>
#include <math_GlobOptMin.hxx>
#include <math_MultipleVarFunction.hxx>
#include <math_Vector.hxx>

#include <cmath>

namespace
{
constexpr double THE_TOLERANCE            = 1.0e-6;
constexpr double THE_STOCHASTIC_TOLERANCE = 1.0e-2;

// ============================================================================
// Test function classes
// ============================================================================

//! Sphere function: f(x) = sum(x_i^2)
struct SphereFunc
{
  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;
  }
};

//! Rastrigin function - highly multimodal
struct RastriginFunc
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const size_t aN  = theX.Size();
    const double aPI = M_PI;
    theF             = 10.0 * aN;
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theF += theX.At(i) * theX.At(i) - 10.0 * std::cos(2.0 * aPI * theX.At(i));
    }
    return true;
  }
};

//! Rosenbrock function
struct RosenbrockFunc
{
  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;
  }
};

//! Booth function
struct BoothFunc
{
  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;
  }
};

//! Box-local minimum at 0.25 and a lower unbounded minimum outside the box.
struct OutsideMinimumFunc
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aX = theX.At(0);
    if (aX < 0.0 || aX > 1.0)
    {
      ++NbOutsideCalls;
      const double anOutsideDelta = aX - 2.0;
      theF                        = anOutsideDelta * anOutsideDelta - 1.0;
      return true;
    }
    const double aDelta = aX - 0.25;
    theF                = aDelta * aDelta;
    return true;
  }

  int NbOutsideCalls = 0;
};

//! Objective with an implicit feasible half-box.
struct PartiallyDefinedFunc
{
  bool Value(const math_Vector& theX, double& theF)
  {
    if (theX.At(0) < 0.0)
    {
      ++NbRejectedCalls;
      return false;
    }
    const double aDelta = theX.At(0) - 0.25;
    theF                = aDelta * aDelta;
    return true;
  }

  int NbRejectedCalls = 0;
};

struct RejectAllFunc
{
  bool Value(const math_Vector&, double&) { return false; }
};

//! Himmelblau function - has 4 identical local minima
//! f(x,y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2
struct HimmelblauFunc
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aX  = theX.At(0);
    const double aY  = theX.At(1);
    const double aT1 = aX * aX + aY - 11.0;
    const double aT2 = aX + aY * aY - 7.0;
    theF             = aT1 * aT1 + aT2 * aT2;
    return true;
  }
};

//! Goldstein-Price function - complex multimodal
struct GoldsteinPriceFunc
{
  bool Value(const math_Vector& theX, double& theF)
  {
    const double aX = theX.At(0);
    const double aY = theX.At(1);

    const double aA =
      1.0
      + (aX + aY + 1.0) * (aX + aY + 1.0)
          * (19.0 - 14.0 * aX + 3.0 * aX * aX - 14.0 * aY + 6.0 * aX * aY + 3.0 * aY * aY);

    const double aB =
      30.0
      + (2.0 * aX - 3.0 * aY) * (2.0 * aX - 3.0 * aY)
          * (18.0 - 32.0 * aX + 12.0 * aX * aX + 48.0 * aY - 36.0 * aX * aY + 27.0 * aY * aY);

    theF = aA * aB;
    return true;
  }
};

// Old API adapter
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;
  }
};

} // namespace

// ============================================================================
// GlobalMinimum unified interface tests
// ============================================================================

TEST(MathOpt_GlobOptTest, GlobalMinimum_PSO_Sphere)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::PSO, 100);
  aConfig.NbPopulation = 40;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_GlobOptTest, GlobalMinimum_MultiStart_Sphere)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::MultiStart, 100);
  aConfig.NbStarts = 20;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_GlobOptTest, GlobalMinimum_PSOHybrid_Booth)
{
  BoothFunc aFunc;

  math_Vector aLower(size_t{2}, -10.0);
  math_Vector aUpper(size_t{2}, 10.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::PSOHybrid, 100);
  aConfig.NbPopulation = 40;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 3.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_GlobOptTest, GlobalMinimum_DE_Sphere)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 100);
  aConfig.NbPopulation  = 40;
  aConfig.MutationScale = 0.8;
  aConfig.CrossoverProb = 0.9;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

// ============================================================================
// Differential Evolution tests
// ============================================================================

TEST(MathOpt_GlobOptTest, DifferentialEvolution_Sphere)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.MaxIterations = 100;
  aConfig.NbPopulation  = 40;
  aConfig.MutationScale = 0.8;
  aConfig.CrossoverProb = 0.9;

  auto aResult = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_GlobOptTest, DifferentialEvolution_Rosenbrock)
{
  RosenbrockFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.MaxIterations = 200;
  aConfig.NbPopulation  = 50;
  aConfig.MutationScale = 0.8;
  aConfig.CrossoverProb = 0.9;

  auto aResult = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 1.0, 0.1);
  EXPECT_NEAR(aResult.Solution->At(1), 1.0, 0.1);
  EXPECT_LT(*aResult.Value, 0.1);
}

TEST(MathOpt_GlobOptTest, DifferentialEvolution_Rastrigin)
{
  RastriginFunc aFunc;

  math_Vector aLower(size_t{2}, -5.12);
  math_Vector aUpper(size_t{2}, 5.12);

  MathOpt::GlobalConfig aConfig;
  aConfig.MaxIterations = 200;
  aConfig.NbPopulation  = 60;
  aConfig.MutationScale = 0.5;
  aConfig.CrossoverProb = 0.9;

  auto aResult = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  // Should find global minimum at origin
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, 0.5);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, 0.5);
  EXPECT_LT(*aResult.Value, 1.0);
}

TEST(MathOpt_GlobOptTest, DifferentialEvolutionRejectsInvalidCandidatesWithoutPoisoningBest)
{
  PartiallyDefinedFunc  aFunc;
  math_Vector           aLower(size_t{1}, -1.0);
  math_Vector           anUpper(size_t{1}, 1.0);
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 100);
  aConfig.NbPopulation       = 20;
  aConfig.PolishBudgetPerDim = 0;

  const MathUtils::VectorResult aResult =
    MathOpt::DifferentialEvolution(aFunc, aLower, anUpper, aConfig);

  EXPECT_GT(aFunc.NbRejectedCalls, 0);
  EXPECT_NE(aResult.Status, MathUtils::Status::CallbackError);
  EXPECT_NE(aResult.Status, MathUtils::Status::NumericalError);
  ASSERT_TRUE(aResult.Solution.has_value());
  ASSERT_TRUE(aResult.Value.has_value());
}

TEST(MathOpt_GlobOptTest, DifferentialEvolutionBudgetExhaustionRetainsBest)
{
  SphereFunc            aFunc;
  math_Vector           aLower(size_t{1}, -1.0);
  math_Vector           anUpper(size_t{1}, 1.0);
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 1);
  aConfig.NbPopulation       = 4;
  aConfig.Tolerance          = 1.0e-30;
  aConfig.PolishBudgetPerDim = 0;

  const MathUtils::VectorResult aResult =
    MathOpt::DifferentialEvolution(aFunc, aLower, anUpper, aConfig);

  EXPECT_EQ(aResult.Status, MathUtils::Status::MaxIterations);
  EXPECT_TRUE(aResult.Solution.has_value());
  EXPECT_TRUE(aResult.Value.has_value());
}

TEST(MathOpt_GlobOptTest, DifferentialEvolutionPreservesAllRejectedCallbackStatus)
{
  RejectAllFunc         aFunc;
  math_Vector           aLower(size_t{1}, -1.0);
  math_Vector           anUpper(size_t{1}, 1.0);
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 10);
  aConfig.NbPopulation = 4;

  EXPECT_EQ(MathOpt::DifferentialEvolution(aFunc, aLower, anUpper, aConfig).Status,
            MathUtils::Status::CallbackError);
}

// ============================================================================
// MultiStart tests
// ============================================================================

TEST(MathOpt_GlobOptTest, MultiStart_Sphere)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.NbStarts = 20;

  auto aResult = MathOpt::MultiStart(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
}

TEST(MathOpt_GlobOptTest, MultiStart_Booth)
{
  BoothFunc aFunc;

  math_Vector aLower(size_t{2}, -10.0);
  math_Vector aUpper(size_t{2}, 10.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.NbStarts = 30;

  auto aResult = MathOpt::MultiStart(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 1.0, 0.5);
  EXPECT_NEAR(aResult.Solution->At(1), 3.0, 0.5);
  EXPECT_LT(*aResult.Value, 1.0);
}

TEST(MathOpt_GlobOptTest, MultiStartEvaluatesAndSolvesInsideBounds)
{
  OutsideMinimumFunc    aFunc;
  math_Vector           aLower(size_t{1}, 0.0);
  math_Vector           anUpper(size_t{1}, 1.0);
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::MultiStart, 100);
  aConfig.NbStarts = 10;

  const MathUtils::VectorResult aResult = MathOpt::MultiStart(aFunc, aLower, anUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  ASSERT_TRUE(aResult.Solution.has_value());
  EXPECT_EQ(aFunc.NbOutsideCalls, 0);
  EXPECT_NEAR(aResult.Solution->At(0), 0.25, 1.0e-4);
  EXPECT_NEAR(*aResult.Value, 0.0, 1.0e-8);
}

// ============================================================================
// Himmelblau function (multiple global minima)
// ============================================================================

TEST(MathOpt_GlobOptTest, Himmelblau_FindsOneMinimum)
{
  HimmelblauFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 200);
  aConfig.NbPopulation = 50;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());

  // Himmelblau has 4 global minima, all with f = 0
  // The algorithm should find one of them
  EXPECT_NEAR(*aResult.Value, 0.0, THE_STOCHASTIC_TOLERANCE);
}

// ============================================================================
// Goldstein-Price function
// ============================================================================

TEST(MathOpt_GlobOptTest, GoldsteinPrice)
{
  GoldsteinPriceFunc aFunc;

  math_Vector aLower(size_t{2}, -2.0);
  math_Vector aUpper(size_t{2}, 2.0);

  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::DifferentialEvolution, 200);
  aConfig.NbPopulation = 60;

  auto aResult = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());

  // Global minimum at (0, -1) with f = 3
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, 0.1);
  EXPECT_NEAR(aResult.Solution->At(1), -1.0, 0.1);
  EXPECT_NEAR(*aResult.Value, 3.0, 0.5);
}

// ============================================================================
// Higher dimensional tests
// ============================================================================

TEST(MathOpt_GlobOptTest, DifferentialEvolution_Sphere5D)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{5}, -5.0);
  math_Vector aUpper(size_t{5}, 5.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.MaxIterations = 300;
  aConfig.NbPopulation  = 80;
  aConfig.MutationScale = 0.7;
  aConfig.CrossoverProb = 0.9;

  auto aResult = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult.IsDone());
  for (int i = 1; i <= 5; ++i)
  {
    EXPECT_NEAR(aResult.Solution->At(static_cast<size_t>(i - 1)), 0.0, 0.1);
  }
  EXPECT_LT(*aResult.Value, 0.1);
}

// ============================================================================
// Strategy comparison tests
// ============================================================================

TEST(MathOpt_GlobOptTest, StrategyComparison_Rosenbrock)
{
  RosenbrockFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  // Test PSO
  MathOpt::GlobalConfig aConfigPSO(MathOpt::GlobalStrategy::PSO, 200);
  aConfigPSO.NbPopulation = 50;

  auto aResultPSO = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfigPSO);

  // Test DE
  MathOpt::GlobalConfig aConfigDE(MathOpt::GlobalStrategy::DifferentialEvolution, 200);
  aConfigDE.NbPopulation = 50;

  auto aResultDE = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfigDE);

  ASSERT_TRUE(aResultPSO.IsDone());
  ASSERT_TRUE(aResultDE.IsDone());

  // Both should find reasonable solutions
  EXPECT_LT(*aResultPSO.Value, 1.0);
  EXPECT_LT(*aResultDE.Value, 1.0);
}

// ============================================================================
// Reproducibility test
// ============================================================================

TEST(MathOpt_GlobOptTest, Reproducibility_DE)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aConfig;
  aConfig.MaxIterations = 100;
  aConfig.NbPopulation  = 40;

  auto aResult1 = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);
  auto aResult2 = MathOpt::DifferentialEvolution(aFunc, aLower, aUpper, aConfig);

  ASSERT_TRUE(aResult1.IsDone());
  ASSERT_TRUE(aResult2.IsDone());

  // With same seed, results should be identical
  EXPECT_DOUBLE_EQ(*aResult1.Value, *aResult2.Value);
}

// ============================================================================
// Comparison with old API
// ============================================================================

TEST(MathOpt_GlobOptTest, CompareWithOldAPI_Sphere)
{
  SphereFuncOld anOldFunc(2);
  SphereFunc    aNewFunc;

  math_Vector anOldLower(1, 2, -5.0);
  math_Vector anOldUpper(1, 2, 5.0);

  math_Vector aNewLower(size_t{2}, -5.0);
  math_Vector aNewUpper(size_t{2}, 5.0);

  // Old API
  math_GlobOptMin anOldSolver(&anOldFunc, anOldLower, anOldUpper, 0.01);
  anOldSolver.Perform();
  ASSERT_TRUE(anOldSolver.isDone());
  ASSERT_GT(anOldSolver.NbExtrema(), 0);
  math_Vector anOldSol(1, 2);
  anOldSolver.Points(1, anOldSol);

  // New API
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::PSOHybrid, 100);
  aConfig.NbPopulation = 40;

  auto aNewResult = MathOpt::GlobalMinimum(aNewFunc, aNewLower, aNewUpper, aConfig);

  ASSERT_TRUE(aNewResult.IsDone());

  EXPECT_LT(anOldSolver.GetF(), 0.1);
  EXPECT_LT(*aNewResult.Value, 0.1);
  EXPECT_NEAR(anOldSolver.GetF(), *aNewResult.Value, 0.1);
  EXPECT_NEAR(anOldSol.At(0), aNewResult.Solution->At(0), 0.1);
  EXPECT_NEAR(anOldSol.At(1), aNewResult.Solution->At(1), 0.1);
}

// ============================================================================
// New GlobalMinimum overload tests
// ============================================================================

TEST(MathOpt_GlobOptTest, GlobalMinimum_PSO_CustomPSOConfig)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  MathOpt::GlobalConfig aGlobConfig(MathOpt::GlobalStrategy::PSO, 200);
  aGlobConfig.NbPopulation = 40;

  MathOpt::PSOConfig aPSOConfig;
  aPSOConfig.NbParticles   = 60;
  aPSOConfig.MaxIterations = 150;
  aPSOConfig.Tolerance     = 1.0e-6;

  aPSOConfig.InertiaSchedule = MathOpt::PSOInertiaSchedule::LinearDecay;
  aPSOConfig.OmegaMin        = 0.3;

  MathOpt::PSOStats aStats;
  auto              aResult =
    MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aGlobConfig, &aPSOConfig, nullptr, &aStats);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 0.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
  EXPECT_GT(aStats.NbFunctionEvals, 0);
  EXPECT_GT(aStats.NbIterations, 0);
}

TEST(MathOpt_GlobOptTest, GlobalMinimum_PSOHybrid_WithSeeds)
{
  BoothFunc aFunc;

  math_Vector aLower(size_t{2}, -10.0);
  math_Vector aUpper(size_t{2}, 10.0);

  MathOpt::GlobalConfig aGlobConfig(MathOpt::GlobalStrategy::PSOHybrid, 100);
  aGlobConfig.NbPopulation = 40;

  MathOpt::PSOConfig aPSOConfig;
  aPSOConfig.NbParticles   = 40;
  aPSOConfig.MaxIterations = 100;
  aPSOConfig.Tolerance     = 1.0e-6;

  aPSOConfig.InitMode = MathOpt::PSOInitMode::SeededPlusRandom;

  // Seed near the minimum (1, 3)
  math_Vector aSeedPos(2);
  aSeedPos.ChangeAt(0) = 1.5;
  aSeedPos.ChangeAt(1) = 2.5;
  NCollection_DynamicArray<MathOpt::PSOSeedParticle> aSeeds;
  aSeeds.Append(MathOpt::PSOSeedParticle(aSeedPos));

  MathOpt::PSOStats aStats;
  auto              aResult =
    MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aGlobConfig, &aPSOConfig, &aSeeds, &aStats);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_NEAR(aResult.Solution->At(0), 1.0, THE_TOLERANCE);
  EXPECT_NEAR(aResult.Solution->At(1), 3.0, THE_TOLERANCE);
  EXPECT_NEAR(*aResult.Value, 0.0, THE_TOLERANCE);
  EXPECT_GT(aStats.NbFunctionEvals, 0);
}

TEST(MathOpt_GlobOptTest, PSOHybridKeepsBetterPSOCandidateAfterClamping)
{
  OutsideMinimumFunc aFunc;
  math_Vector        aLower(size_t{1}, 0.0);
  math_Vector        anUpper(size_t{1}, 1.0);

  math_Vector                                        aSeedPosition(size_t{1}, 0.25);
  NCollection_DynamicArray<MathOpt::PSOSeedParticle> aSeeds;
  aSeeds.Append(MathOpt::PSOSeedParticle(aSeedPosition, 0.0));

  MathOpt::PSOConfig aPSOConfig;
  aPSOConfig.NbParticles        = 1;
  aPSOConfig.MaxIterations      = 1;
  aPSOConfig.InitMode           = MathOpt::PSOInitMode::SeededOnly;
  aPSOConfig.VelocityClamp      = 0.0;
  aPSOConfig.PhiPersonal        = 0.0;
  aPSOConfig.PhiGlobal          = 0.0;
  aPSOConfig.PolishBudgetPerDim = 0;

  MathOpt::GlobalConfig         aConfig(MathOpt::GlobalStrategy::PSOHybrid, 2);
  const MathUtils::VectorResult aResult =
    MathOpt::GlobalMinimum(aFunc, aLower, anUpper, aConfig, &aPSOConfig, &aSeeds);

  ASSERT_TRUE(aResult.IsDone());
  EXPECT_DOUBLE_EQ(aResult.Solution->At(0), 0.25);
  EXPECT_DOUBLE_EQ(*aResult.Value, 0.0);
  EXPECT_EQ(aFunc.NbOutsideCalls, 0);
}

TEST(MathOpt_GlobOptTest, PSOHybridAcceptsSingleIterationBudget)
{
  SphereFunc            aFunc;
  math_Vector           aLower(size_t{1}, -1.0);
  math_Vector           anUpper(size_t{1}, 1.0);
  MathOpt::GlobalConfig aConfig(MathOpt::GlobalStrategy::PSOHybrid, 1);
  aConfig.NbPopulation       = 4;
  aConfig.PolishBudgetPerDim = 0;

  const MathUtils::VectorResult aResult = MathOpt::GlobalMinimum(aFunc, aLower, anUpper, aConfig);

  EXPECT_NE(aResult.Status, MathUtils::Status::InvalidInput);
}

TEST(MathOpt_GlobOptTest, GlobalMinimum_NonPSOStrategies_UnchangedWithPSOOptions)
{
  SphereFunc aFunc;

  math_Vector aLower(size_t{2}, -5.0);
  math_Vector aUpper(size_t{2}, 5.0);

  // MultiStart - PSO params should be ignored
  MathOpt::GlobalConfig aConfigMS(MathOpt::GlobalStrategy::MultiStart, 100);
  aConfigMS.NbStarts = 20;

  MathOpt::PSOConfig aPSOConfig;
  aPSOConfig.NbParticles = 100;

  auto aResultMS = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfigMS, &aPSOConfig);
  ASSERT_TRUE(aResultMS.IsDone());
  EXPECT_NEAR(*aResultMS.Value, 0.0, THE_TOLERANCE);

  // DE - PSO params should be ignored
  MathOpt::GlobalConfig aConfigDE(MathOpt::GlobalStrategy::DifferentialEvolution, 100);
  aConfigDE.NbPopulation = 40;

  auto aResultDE = MathOpt::GlobalMinimum(aFunc, aLower, aUpper, aConfigDE, &aPSOConfig);
  ASSERT_TRUE(aResultDE.IsDone());
  EXPECT_NEAR(*aResultDE.Value, 0.0, THE_TOLERANCE);
}