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 _MathUtils_FunctorVector_HeaderFile
#define _MathUtils_FunctorVector_HeaderFile

#include <math_Vector.hxx>
#include <math_Matrix.hxx>
#include <MathUtils_Core.hxx>

#include <cmath>
#include <utility>

//! @file MathUtils_FunctorVector.hxx
//! @brief Non-virtual functor classes for N-dimensional (vector) functions.
//!
//! Provides ready-to-use functor classes that work with the template-based
//! math API (MathOpt::Powell, MathOpt::BFGS, MathSys::Newton) without
//! virtual dispatch overhead.

//! Core utilities for modern math solvers.
namespace MathUtils
{

//! Lambda wrapper for N-D objective functions (value only).
//! Wraps a lambda/callable into a functor with Value() method.
//!
//! Usage:
//! @code
//! auto aFunc = VectorLambda([](const math_Vector& x, double& y) {
//!   y = x.At(0)*x.At(0) + x.At(1)*x.At(1);  // Sphere function
//!   return true;
//! });
//! auto aResult = MathOpt::Powell(aFunc, aStartPoint);
//! @endcode
//!
//! @tparam Lambda callable type with signature bool(const math_Vector&, double&)
template <typename Lambda>
class VectorLambda
{
public:
  //! Constructor from lambda/callable.
  //! @param theLambda callable with signature bool(const math_Vector&, double&)
  explicit VectorLambda(Lambda theLambda)
      : myLambda(std::move(theLambda))
  {
  }

  //! Evaluates the function at theX.
  //! @param[in] theX input vector
  //! @param[out] theY function value f(theX)
  //! @return true if evaluation succeeded
  bool Value(const math_Vector& theX, double& theY) const { return myLambda(theX, theY); }

private:
  Lambda myLambda;
};

//! Lambda wrapper for N-D objective functions with gradient.
//! Wraps a lambda/callable into a functor with Value() and Gradient() methods.
//!
//! Usage:
//! @code
//! auto aFunc = VectorLambdaWithGradient(
//!   [](const math_Vector& x, double& y) {
//!     y = x.At(0)*x.At(0) + x.At(1)*x.At(1);
//!     return true;
//!   },
//!   [](const math_Vector& x, math_Vector& g) {
//!     g.ChangeAt(0) = 2.0 * x.At(0);
//!     g.ChangeAt(1) = 2.0 * x.At(1);
//!     return true;
//!   });
//! auto aResult = MathOpt::BFGS(aFunc, aStartPoint);
//! @endcode
//!
//! @tparam ValueLambda callable for value: bool(const math_Vector&, double&)
//! @tparam GradLambda callable for gradient: bool(const math_Vector&, math_Vector&)
template <typename ValueLambda, typename GradLambda>
class VectorLambdaWithGradient
{
public:
  //! Constructor from value and gradient lambdas.
  //! @param theValueLambda callable for value evaluation
  //! @param theGradLambda callable for gradient evaluation
  VectorLambdaWithGradient(ValueLambda theValueLambda, GradLambda theGradLambda)
      : myValueLambda(std::move(theValueLambda)),
        myGradLambda(std::move(theGradLambda))
  {
  }

  //! Evaluates the function value at theX.
  //! @param[in] theX input vector
  //! @param[out] theY function value f(theX)
  //! @return true if evaluation succeeded
  bool Value(const math_Vector& theX, double& theY) const { return myValueLambda(theX, theY); }

  //! Evaluates the gradient at theX.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    return myGradLambda(theX, theG);
  }

  //! Evaluates both value and gradient at theX.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return myValueLambda(theX, theY) && myGradLambda(theX, theG);
  }

private:
  ValueLambda myValueLambda;
  GradLambda  myGradLambda;
};

//! Helper function to create VectorLambda with type deduction.
//! @tparam Lambda callable type (auto-deduced)
//! @param theLambda callable with signature bool(const math_Vector&, double&)
//! @return VectorLambda wrapper
template <typename Lambda>
VectorLambda<Lambda> MakeVector(Lambda theLambda)
{
  return VectorLambda<Lambda>(std::move(theLambda));
}

//! Helper function to create VectorLambdaWithGradient with type deduction.
//! @tparam ValueLambda callable for value
//! @tparam GradLambda callable for gradient
//! @param theValueLambda callable for value evaluation
//! @param theGradLambda callable for gradient evaluation
//! @return VectorLambdaWithGradient wrapper
template <typename ValueLambda, typename GradLambda>
VectorLambdaWithGradient<ValueLambda, GradLambda> MakeVectorWithGradient(ValueLambda theValueLambda,
                                                                         GradLambda  theGradLambda)
{
  return VectorLambdaWithGradient<ValueLambda, GradLambda>(std::move(theValueLambda),
                                                           std::move(theGradLambda));
}

//! Quadratic form functor: f(x) = x^T A x + b^T x + c.
//! Commonly used for testing optimization algorithms.
//!
//! Usage:
//! @code
//! math_Matrix A(2, 2);
//! A.ChangeAt(0, 0) = 2.0; A.ChangeAt(0, 1) = 0.0;
//! A.ChangeAt(1, 0) = 0.0; A.ChangeAt(1, 1) = 2.0;
//! math_Vector b(2);
//! b.ChangeAt(0) = -4.0; b.ChangeAt(1) = -4.0;
//! QuadraticForm aFunc(A, b, 8.0);  // f(x) = 2*x1^2 + 2*x2^2 - 4*x1 - 4*x2 + 8
//! // Minimum at (1, 1) with value 4
//! @endcode
class QuadraticForm
{
public:
  //! Constructor from matrix, vector, and constant.
  //! @param theA quadratic coefficient matrix (must be square)
  //! @param theB linear coefficient vector
  //! @param theC constant term
  QuadraticForm(const math_Matrix& theA, const math_Vector& theB, double theC)
      : myA(theA),
        myB(theB),
        myC(theC)
  {
  }

  //! Evaluates the quadratic form f(x) = x^T A x + b^T x + c.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @return true if evaluation succeeded
  bool Value(const math_Vector& theX, double& theY) const
  {
    theY = myC;
    // x^T A x
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      for (size_t j = 0; j < theX.Size(); ++j)
      {
        theY += theX.At(i) * myA.At(i, j) * theX.At(j);
      }
    }
    // b^T x
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theY += myB.At(i) * theX.At(i);
    }
    return true;
  }

  //! Evaluates the gradient: g = 2*A*x + b.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    // g = (A + A^T) * x + b = 2*A*x + b (for symmetric A)
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theG.ChangeAt(i) = myB.At(i);
      for (size_t j = 0; j < theX.Size(); ++j)
      {
        theG.ChangeAt(i) += (myA.At(i, j) + myA.At(j, i)) * theX.At(j);
      }
    }
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }

private:
  math_Matrix myA;
  math_Vector myB;
  double      myC;
};

//! Rosenbrock function functor (for testing optimization).
//! f(x,y) = (a - x)^2 + b*(y - x^2)^2
//! Default: a = 1, b = 100
//! Global minimum at (a, a^2) = (1, 1) with f = 0.
//!
//! Usage:
//! @code
//! Rosenbrock aRosen;  // Default a=1, b=100
//! math_Vector aStart(2);
//! aStart.ChangeAt(0) = -1.0; aStart.ChangeAt(1) = 1.0;
//! auto aResult = MathOpt::BFGS(aRosen, aStart);
//! // Should converge to (1, 1)
//! @endcode
class Rosenbrock
{
public:
  //! Constructor with parameters.
  //! @param theA parameter a (default 1.0)
  //! @param theB parameter b (default 100.0)
  Rosenbrock(double theA = 1.0, double theB = 100.0)
      : myA(theA),
        myB(theB)
  {
  }

  //! Evaluates the Rosenbrock function.
  //! @param[in] theX input vector (must have length 2)
  //! @param[out] theY function value
  //! @return true if evaluation succeeded
  bool Value(const math_Vector& theX, double& theY) const
  {
    const double x  = theX.At(0);
    const double y  = theX.At(1);
    const double t1 = myA - x;
    const double t2 = y - x * x;
    theY            = t1 * t1 + myB * t2 * t2;
    return true;
  }

  //! Evaluates the gradient of the Rosenbrock function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    const double x   = theX.At(0);
    const double y   = theX.At(1);
    const double t2  = y - x * x;
    theG.ChangeAt(0) = -2.0 * (myA - x) - 4.0 * myB * x * t2;
    theG.ChangeAt(1) = 2.0 * myB * t2;
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true if evaluation succeeded
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }

private:
  double myA;
  double myB;
};

//! Sphere function functor (for testing optimization).
//! f(x) = sum(x[i]^2) for all i.
//! Global minimum at origin with f = 0.
//!
//! Usage:
//! @code
//! Sphere aSphere;
//! math_Vector aStart(3);
//! aStart.Init(1.0);  // Start at (1, 1, 1)
//! auto aResult = MathOpt::Powell(aSphere, aStart);
//! // Should converge to (0, 0, 0)
//! @endcode
class Sphere
{
public:
  //! Evaluates the sphere function.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    theY = 0.0;
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theY += theX.At(i) * theX.At(i);
    }
    return true;
  }

  //! Evaluates the gradient of the sphere function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theG.ChangeAt(i) = 2.0 * theX.At(i);
    }
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }
};

//! Booth function functor (for testing optimization).
//! f(x,y) = (x + 2y - 7)^2 + (2x + y - 5)^2
//! Global minimum at (1, 3) with f = 0.
//!
//! Usage:
//! @code
//! Booth aBooth;
//! math_Vector aStart(2);
//! aStart.ChangeAt(0) = 0.0; aStart.ChangeAt(1) = 0.0;
//! auto aResult = MathOpt::BFGS(aBooth, aStart);
//! // Should converge to (1, 3)
//! @endcode
class Booth
{
public:
  //! Evaluates the Booth function.
  //! @param[in] theX input vector (must have length 2)
  //! @param[out] theY function value
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    const double x  = theX.At(0);
    const double y  = theX.At(1);
    const double t1 = x + 2.0 * y - 7.0;
    const double t2 = 2.0 * x + y - 5.0;
    theY            = t1 * t1 + t2 * t2;
    return true;
  }

  //! Evaluates the gradient of the Booth function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    const double x   = theX.At(0);
    const double y   = theX.At(1);
    const double t1  = x + 2.0 * y - 7.0;
    const double t2  = 2.0 * x + y - 5.0;
    theG.ChangeAt(0) = 2.0 * t1 + 4.0 * t2;
    theG.ChangeAt(1) = 4.0 * t1 + 2.0 * t2;
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }
};

//! Beale function functor (for testing optimization).
//! f(x,y) = (1.5 - x + xy)^2 + (2.25 - x + xy^2)^2 + (2.625 - x + xy^3)^2
//! Global minimum at (3, 0.5) with f = 0.
class Beale
{
public:
  //! Evaluates the Beale function.
  //! @param[in] theX input vector (must have length 2)
  //! @param[out] theY function value
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    const double x  = theX.At(0);
    const double y  = theX.At(1);
    const double t1 = 1.5 - x + x * y;
    const double t2 = 2.25 - x + x * y * y;
    const double t3 = 2.625 - x + x * y * y * y;
    theY            = t1 * t1 + t2 * t2 + t3 * t3;
    return true;
  }

  //! Evaluates the gradient of the Beale function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    const double x   = theX.At(0);
    const double y   = theX.At(1);
    const double y2  = y * y;
    const double y3  = y2 * y;
    const double t1  = 1.5 - x + x * y;
    const double t2  = 2.25 - x + x * y2;
    const double t3  = 2.625 - x + x * y3;
    theG.ChangeAt(0) = 2.0 * ((y - 1.0) * t1 + (y2 - 1.0) * t2 + (y3 - 1.0) * t3);
    theG.ChangeAt(1) = 2.0 * x * (t1 + 2.0 * y * t2 + 3.0 * y2 * t3);
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }
};

//! Himmelblau function functor (for testing optimization).
//! f(x,y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2
//! Has four local minima, all with f = 0:
//! (3.0, 2.0), (-2.805118, 3.131312), (-3.779310, -3.283186), (3.584428, -1.848126)
class Himmelblau
{
public:
  //! Evaluates the Himmelblau function.
  //! @param[in] theX input vector (must have length 2)
  //! @param[out] theY function value
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    const double x  = theX.At(0);
    const double y  = theX.At(1);
    const double t1 = x * x + y - 11.0;
    const double t2 = x + y * y - 7.0;
    theY            = t1 * t1 + t2 * t2;
    return true;
  }

  //! Evaluates the gradient of the Himmelblau function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    const double x   = theX.At(0);
    const double y   = theX.At(1);
    const double t1  = x * x + y - 11.0;
    const double t2  = x + y * y - 7.0;
    theG.ChangeAt(0) = 4.0 * x * t1 + 2.0 * t2;
    theG.ChangeAt(1) = 2.0 * t1 + 4.0 * y * t2;
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }
};

//! Rastrigin function functor (for testing global optimization).
//! f(x) = A*n + sum(x[i]^2 - A*cos(2*pi*x[i])) for all i
//! Default: A = 10
//! Global minimum at origin with f = 0.
//! Highly multimodal - challenging for local optimizers.
class Rastrigin
{
public:
  //! Constructor with parameter.
  //! @param theA parameter A (default 10.0)
  explicit Rastrigin(double theA = 10.0)
      : myA(theA)
  {
  }

  //! Evaluates the Rastrigin function.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    const size_t n = theX.Size();
    theY           = myA * static_cast<double>(n);
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theY += theX.At(i) * theX.At(i) - myA * std::cos(THE_2PI * theX.At(i));
    }
    return true;
  }

  //! Evaluates the gradient of the Rastrigin function.
  //! @param[in] theX input vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    for (size_t i = 0; i < theX.Size(); ++i)
    {
      theG.ChangeAt(i) = 2.0 * theX.At(i) + myA * THE_2PI * std::sin(THE_2PI * theX.At(i));
    }
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }

private:
  double myA;
};

//! Ackley function functor (for testing global optimization).
//! f(x) = -a*exp(-b*sqrt(sum(x[i]^2)/n)) - exp(sum(cos(c*x[i]))/n) + a + e
//! Default: a = 20, b = 0.2, c = 2*pi
//! Global minimum at origin with f = 0.
class Ackley
{
public:
  //! Constructor with parameters.
  //! @param theA parameter a (default 20.0)
  //! @param theB parameter b (default 0.2)
  //! @param theC parameter c (default 2*pi)
  Ackley(double theA = 20.0, double theB = 0.2, double theC = THE_2PI)
      : myA(theA),
        myB(theB),
        myC(theC)
  {
  }

  //! Evaluates the Ackley function.
  //! @param[in] theX input vector
  //! @param[out] theY function value
  //! @return false for empty or non-finite input, parameters, or result
  bool Value(const math_Vector& theX, double& theY) const
  {
    const size_t aSize = theX.Size();
    if (aSize == 0 || !std::isfinite(myA) || !std::isfinite(myB) || !std::isfinite(myC))
    {
      return false;
    }

    double aSumSq  = 0.0;
    double aSumCos = 0.0;
    for (size_t i = 0; i < aSize; ++i)
    {
      const double aValue = theX.At(i);
      if (!std::isfinite(aValue) || !std::isfinite(myC * aValue))
      {
        return false;
      }
      aSumSq += aValue * aValue;
      aSumCos += std::cos(myC * aValue);
      if (!std::isfinite(aSumSq) || !std::isfinite(aSumCos))
      {
        return false;
      }
    }
    const double aDivisor = static_cast<double>(aSize);
    theY = -myA * std::exp(-myB * std::sqrt(aSumSq / aDivisor)) - std::exp(aSumCos / aDivisor) + myA
           + THE_EULER_NUMBER;
    return std::isfinite(theY);
  }

private:
  double myA;
  double myB;
  double myC;
};

//! Linear system residual functor: f(x) = ||Ax - b||^2.
//! Useful for solving overdetermined linear systems via optimization.
//!
//! Usage:
//! @code
//! math_Matrix A(3, 2);  // 3x2 overdetermined system
//! math_Vector b(3);
//! // ... fill A and b ...
//! LinearResidual aRes(A, b);
//! math_Vector aStart(2);
//! aStart.Init(0.0);
//! auto aResult = MathOpt::BFGS(aRes, aStart);
//! @endcode
class LinearResidual
{
public:
  //! Constructor from matrix and right-hand side.
  //! @param theA coefficient matrix (m x n)
  //! @param theB right-hand side vector (m)
  LinearResidual(const math_Matrix& theA, const math_Vector& theB)
      : myA(theA),
        myB(theB)
  {
  }

  //! Evaluates the residual ||Ax - b||^2.
  //! @param[in] theX solution vector (n)
  //! @param[out] theY squared residual norm
  //! @return true (always succeeds)
  bool Value(const math_Vector& theX, double& theY) const
  {
    theY = 0.0;
    for (size_t i = 0; i < myA.RowSize(); ++i)
    {
      double aResidual = -myB.At(i);
      for (size_t j = 0; j < myA.ColSize(); ++j)
      {
        aResidual += myA.At(i, j) * theX.At(j);
      }
      theY += aResidual * aResidual;
    }
    return true;
  }

  //! Evaluates the gradient: g = 2 * A^T * (Ax - b).
  //! @param[in] theX solution vector
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Gradient(const math_Vector& theX, math_Vector& theG) const
  {
    // First compute residual r = Ax - b
    math_Vector aResidual(myA.RowSize());
    for (size_t i = 0; i < myA.RowSize(); ++i)
    {
      aResidual.ChangeAt(i) = -myB.At(i);
      for (size_t j = 0; j < myA.ColSize(); ++j)
      {
        aResidual.ChangeAt(i) += myA.At(i, j) * theX.At(j);
      }
    }
    // Then compute g = 2 * A^T * r
    for (size_t j = 0; j < myA.ColSize(); ++j)
    {
      theG.ChangeAt(j) = 0.0;
      for (size_t i = 0; i < myA.RowSize(); ++i)
      {
        theG.ChangeAt(j) += 2.0 * myA.At(i, j) * aResidual.At(i);
      }
    }
    return true;
  }

  //! Evaluates both value and gradient.
  //! @param[in] theX solution vector
  //! @param[out] theY squared residual norm
  //! @param[out] theG gradient vector
  //! @return true (always succeeds)
  bool Values(const math_Vector& theX, double& theY, math_Vector& theG) const
  {
    return Value(theX, theY) && Gradient(theX, theG);
  }

private:
  math_Matrix myA;
  math_Vector myB;
};

//! Nonlinear system functor: F(x) = [f1(x), f2(x), ..., fn(x)].
//! Lambda wrapper for systems of nonlinear equations.
//!
//! Usage:
//! @code
//! auto aSys = SystemLambda([](const math_Vector& x, math_Vector& f) {
//!   f.ChangeAt(0) = x.At(0)*x.At(0) + x.At(1)*x.At(1) - 1.0;  // Circle
//!   f.ChangeAt(1) = x.At(0) - x.At(1);                         // Line y=x
//!   return true;
//! });
//! auto aResult = MathSys::Newton(aSys, aStartPoint);
//! @endcode
//!
//! @tparam Lambda callable with signature bool(const math_Vector&, math_Vector&)
template <typename Lambda>
class SystemLambda
{
public:
  //! Constructor from lambda.
  //! @param theLambda callable for system evaluation
  //! @param theNbEquations number of equations in the system
  SystemLambda(Lambda theLambda, size_t theNbEquations)
      : myLambda(std::move(theLambda)),
        myNbEquations(theNbEquations)
  {
  }

  //! Returns the number of equations.
  //! @return number of equations
  size_t NbEquations() const { return myNbEquations; }

  //! Evaluates the system F(x).
  //! @param[in] theX input vector
  //! @param[out] theF system values
  //! @return true if evaluation succeeded
  bool Value(const math_Vector& theX, math_Vector& theF) const { return myLambda(theX, theF); }

private:
  Lambda myLambda;
  size_t myNbEquations;
};

//! Helper function to create SystemLambda with type deduction.
//! @tparam Lambda callable type
//! @param theLambda callable for system evaluation
//! @param theNbEquations number of equations
//! @return SystemLambda wrapper
template <typename Lambda>
SystemLambda<Lambda> MakeSystem(Lambda theLambda, size_t theNbEquations)
{
  return SystemLambda<Lambda>(std::move(theLambda), theNbEquations);
}

} // namespace MathUtils

#endif // _MathUtils_FunctorVector_HeaderFile