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.

#ifndef _MathOpt_Brent_HeaderFile
#define _MathOpt_Brent_HeaderFile

#include <MathUtils_Types.hxx>
#include <MathUtils_Config.hxx>
#include <MathUtils_Core.hxx>
#include <MathUtils_Bracket.hxx>
#include <MathOpt_Utils.hxx>

#include <cmath>

//! Optimization algorithms for scalar and vector functions.
namespace MathOpt
{
using namespace MathUtils;

//! Brent's method for 1D minimization.
//! Combines golden section search with parabolic interpolation.
//! Guaranteed to converge for unimodal functions within the given interval.
//!
//! Algorithm:
//! 1. Maintain bracket [a, b] with interior point x where f(x) < f(a), f(x) < f(b)
//! 2. Try parabolic interpolation using three points
//! 3. If parabolic step is rejected, use golden section step
//! 4. Update bracket and repeat until convergence
//!
//! @tparam Function type with Value(double theX, double& theF) method
//! @param theFunc function to minimize
//! @param theLower lower bound of search interval
//! @param theUpper upper bound of search interval
//! @param theConfig solver configuration
//! @return result containing minimum location and value
template <typename Function>
ScalarResult Brent(Function&     theFunc,
                   double        theLower,
                   double        theUpper,
                   const Config& theConfig = Config())
{
  ScalarResult aResult;

  if (!Utils::IsValidConfig(theConfig) || !std::isfinite(theLower) || !std::isfinite(theUpper)
      || theLower >= theUpper)
  {
    aResult.Status = Status::InvalidInput;
    return aResult;
  }

  double aA = theLower;
  double aB = theUpper;

  // Initial point using golden section
  double aX = aA + MathUtils::THE_GOLDEN_SECTION * (aB - aA);
  double aW = aX;
  double aV = aX;

  double       aFx          = 0.0;
  const Status aValueStatus = Utils::ValueStatus(theFunc, aX, aFx);
  if (aValueStatus != Status::OK)
  {
    aResult.Status = aValueStatus;
    return aResult;
  }
  double aFw = aFx;
  double aFv = aFx;

  double aD = 0.0; // Current step
  double aE = 0.0; // Previous step

  for (uint32_t anIter = 0; anIter < theConfig.MaxIterations; ++anIter)
  {
    const double aXm   = 0.5 * (aA + aB);
    const double aTol1 = theConfig.XTolerance * std::abs(aX) + MathUtils::THE_ZERO_TOL / 10.0;
    const double aTol2 = 2.0 * aTol1;

    aResult.NbIterations = anIter + 1;

    // Check convergence
    if (std::abs(aX - aXm) <= (aTol2 - 0.5 * (aB - aA)))
    {
      aResult.Status = Status::OK;
      aResult.Root   = aX;
      aResult.Value  = aFx;
      return aResult;
    }

    double aU            = 0.0;
    bool   aUseParabolic = false;

    // Try parabolic interpolation if step is large enough
    if (std::abs(aE) > aTol1)
    {
      // Parabolic fit through x, w, v
      const double aR = (aX - aW) * (aFx - aFv);
      double       aQ = (aX - aV) * (aFx - aFw);
      double       aP = (aX - aV) * aQ - (aX - aW) * aR;
      aQ              = 2.0 * (aQ - aR);

      if (aQ > 0.0)
      {
        aP = -aP;
      }
      else
      {
        aQ = -aQ;
      }

      const double aETmp = aE;
      aE                 = aD;

      // Check if parabolic step is acceptable
      if (std::abs(aP) < std::abs(0.5 * aQ * aETmp) && aP > aQ * (aA - aX) && aP < aQ * (aB - aX))
      {
        aD = aP / aQ;
        aU = aX + aD;

        // Don't evaluate too close to bounds
        if ((aU - aA) < aTol2 || (aB - aU) < aTol2)
        {
          aD = MathUtils::SignTransfer(aTol1, aXm - aX);
        }
        aUseParabolic = true;
      }
    }

    if (!aUseParabolic)
    {
      // Golden section step
      aE = (aX < aXm) ? (aB - aX) : (aA - aX);
      aD = MathUtils::THE_GOLDEN_SECTION * aE;
    }

    // Ensure step is at least aTol1
    if (std::abs(aD) >= aTol1)
    {
      aU = aX + aD;
    }
    else
    {
      aU = aX + MathUtils::SignTransfer(aTol1, aD);
    }

    double       aFu               = 0.0;
    const Status aProbeValueStatus = Utils::ValueStatus(theFunc, aU, aFu);
    if (aProbeValueStatus != Status::OK)
    {
      aResult.Status = aProbeValueStatus;
      aResult.Root   = aX;
      aResult.Value  = aFx;
      return aResult;
    }

    // Update bracket and best points
    if (aFu <= aFx)
    {
      if (aU < aX)
      {
        aB = aX;
      }
      else
      {
        aA = aX;
      }

      aV  = aW;
      aW  = aX;
      aX  = aU;
      aFv = aFw;
      aFw = aFx;
      aFx = aFu;
    }
    else
    {
      if (aU < aX)
      {
        aA = aU;
      }
      else
      {
        aB = aU;
      }

      if (aFu <= aFw || aW == aX)
      {
        aV  = aW;
        aW  = aU;
        aFv = aFw;
        aFw = aFu;
      }
      else if (aFu <= aFv || aV == aX || aV == aW)
      {
        aV  = aU;
        aFv = aFu;
      }
    }
  }

  // Maximum iterations reached
  aResult.Status = Status::MaxIterations;
  aResult.Root   = aX;
  aResult.Value  = aFx;
  return aResult;
}

//! Golden section search for 1D minimization.
//! Simpler than Brent but with guaranteed linear convergence.
//! Does not attempt parabolic interpolation.
//!
//! @tparam Function type with Value(double theX, double& theF) method
//! @param theFunc function to minimize
//! @param theLower lower bound of search interval
//! @param theUpper upper bound of search interval
//! @param theConfig solver configuration
//! @return result containing minimum location and value
template <typename Function>
ScalarResult Golden(Function&     theFunc,
                    double        theLower,
                    double        theUpper,
                    const Config& theConfig = Config())
{
  ScalarResult aResult;

  if (!Utils::IsValidConfig(theConfig) || !std::isfinite(theLower) || !std::isfinite(theUpper)
      || theLower >= theUpper)
  {
    aResult.Status = Status::InvalidInput;
    return aResult;
  }

  double aA = theLower;
  double aB = theUpper;

  // Initialize interior points
  double aX1 = aA + MathUtils::THE_GOLDEN_SECTION * (aB - aA);
  double aX2 = aA + (1.0 - MathUtils::THE_GOLDEN_SECTION) * (aB - aA);

  double aF1 = 0.0;
  double aF2 = 0.0;

  const Status aFirstValueStatus = Utils::ValueStatus(theFunc, aX1, aF1);
  if (aFirstValueStatus != Status::OK)
  {
    aResult.Status = aFirstValueStatus;
    return aResult;
  }
  const Status aSecondValueStatus = Utils::ValueStatus(theFunc, aX2, aF2);
  if (aSecondValueStatus != Status::OK)
  {
    aResult.Status = aSecondValueStatus;
    return aResult;
  }

  for (uint32_t anIter = 0; anIter < theConfig.MaxIterations; ++anIter)
  {
    aResult.NbIterations = anIter + 1;

    // Check convergence
    if ((aB - aA) < theConfig.XTolerance * (std::abs(aX1) + std::abs(aX2)))
    {
      aResult.Status = Status::OK;
      if (aF1 < aF2)
      {
        aResult.Root  = aX1;
        aResult.Value = aF1;
      }
      else
      {
        aResult.Root  = aX2;
        aResult.Value = aF2;
      }
      return aResult;
    }

    if (aF1 < aF2)
    {
      // Minimum is in [a, x2]
      aB                        = aX2;
      aX2                       = aX1;
      aF2                       = aF1;
      aX1                       = aA + MathUtils::THE_GOLDEN_SECTION * (aB - aA);
      const Status aValueStatus = Utils::ValueStatus(theFunc, aX1, aF1);
      if (aValueStatus != Status::OK)
      {
        aResult.Status = aValueStatus;
        aResult.Root   = aX2;
        aResult.Value  = aF2;
        return aResult;
      }
    }
    else
    {
      // Minimum is in [x1, b]
      aA                        = aX1;
      aX1                       = aX2;
      aF1                       = aF2;
      aX2                       = aA + (1.0 - MathUtils::THE_GOLDEN_SECTION) * (aB - aA);
      const Status aValueStatus = Utils::ValueStatus(theFunc, aX2, aF2);
      if (aValueStatus != Status::OK)
      {
        aResult.Status = aValueStatus;
        aResult.Root   = aX1;
        aResult.Value  = aF1;
        return aResult;
      }
    }
  }

  // Maximum iterations reached
  aResult.Status = Status::MaxIterations;
  if (aF1 < aF2)
  {
    aResult.Root  = aX1;
    aResult.Value = aF1;
  }
  else
  {
    aResult.Root  = aX2;
    aResult.Value = aF2;
  }
  return aResult;
}

//! Brent's method with automatic bracket search.
//! First attempts to bracket a minimum, then applies Brent's method.
//!
//! @tparam Function type with Value(double theX, double& theF) method
//! @param theFunc function to minimize
//! @param theGuess initial guess
//! @param theStep initial step size for bracket search
//! @param theConfig solver configuration
//! @return result containing minimum location and value
template <typename Function>
ScalarResult BrentWithBracket(Function&     theFunc,
                              double        theGuess,
                              double        theStep   = 1.0,
                              const Config& theConfig = Config())
{
  ScalarResult aResult;

  if (!Utils::IsValidConfig(theConfig) || !std::isfinite(theGuess) || !std::isfinite(theStep)
      || theStep == 0.0 || !std::isfinite(theGuess + theStep))
  {
    aResult.Status = Status::InvalidInput;
    return aResult;
  }

  // Try to bracket minimum
  MathUtils::MinBracketResult aBracket =
    MathUtils::BracketMinimum(theFunc, theGuess, theGuess + theStep);

  if (!aBracket.IsValid)
  {
    aResult.Status = aBracket.Status;
    return aResult;
  }
  if (!std::isfinite(aBracket.A) || !std::isfinite(aBracket.C) || aBracket.A == aBracket.C)
  {
    aResult.Status = Status::NumericalError;
    return aResult;
  }

  // Apply Brent's method on the bracket
  return Brent(theFunc, aBracket.A, aBracket.C, theConfig);
}

} // namespace MathOpt

#endif // _MathOpt_Brent_HeaderFile