EIC code displayed by LXR master/​include/​opencascade/​MathUtils_Deriv.hxx	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-09-28 09:20:53

 // 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_Deriv_HeaderFile
 #define _MathUtils_Deriv_HeaderFile
 
 #include <math_Vector.hxx>
 #include <math_Matrix.hxx>
 #include <MathUtils_Core.hxx>
 
 #include <cmath>
 
 //! Core utilities for modern math solvers.
 namespace MathUtils
 {
 
 //! Central difference derivative approximation for scalar functions.
 //! f'(x) ~= (f(x+h) - f(x-h)) / (2h)
 //! Accuracy: O(h^2)
 //!
 //! @tparam Function type with Value(double theX, double& theF) method
 //! @param theFunc function to differentiate
 //! @param theX point at which to evaluate derivative
 //! @param theDeriv computed derivative value
 //! @param theStep step size (default 1e-8)
 //! @return true if successful
 template <typename Function>
 bool CentralDifference(Function& theFunc, double theX, double& theDeriv, double theStep = 1.0e-8)
 {
   double aFPlus  = 0.0;
   double aFMinus = 0.0;
 
   if (!theFunc.Value(theX + theStep, aFPlus))
   {
     return false;
   }
   if (!theFunc.Value(theX - theStep, aFMinus))
   {
     return false;
   }
 
   theDeriv = (aFPlus - aFMinus) / (2.0 * theStep);
   return true;
 }
 
 //! Forward difference derivative (one-sided).
 //! f'(x) ~= (f(x+h) - f(x)) / h
 //! Accuracy: O(h)
 //! Useful when central difference is not possible (e.g., at boundaries).
 //!
 //! @tparam Function type with Value(double theX, double& theF) method
 //! @param theFunc function to differentiate
 //! @param theX point at which to evaluate derivative
 //! @param theFx function value at theX (if already computed)
 //! @param theDeriv computed derivative value
 //! @param theStep step size (default 1e-8)
 //! @return true if successful
 template <typename Function>
 bool ForwardDifference(Function& theFunc,
                        double    theX,
                        double    theFx,
                        double&   theDeriv,
                        double    theStep = 1.0e-8)
 {
   double aFPlus = 0.0;
 
   if (!theFunc.Value(theX + theStep, aFPlus))
   {
     return false;
   }
 
   theDeriv = (aFPlus - theFx) / theStep;
   return true;
 }
 
 //! Numerical gradient using central differences for N-D functions.
 //! df/dx[i] ~= (f(x + h[i]*e[i]) - f(x - h[i]*e[i])) / (2*h[i])
 //!
 //! @tparam Function type with Value(const math_Vector&, double&) method
 //! @param theFunc function to differentiate
 //! @param theX point at which to evaluate gradient (temporarily modified)
 //! @param theGrad output gradient vector (same dimension as theX)
 //! @param theStep step size (default 1e-8)
 //! @return true if successful
 template <typename Function>
 bool NumericalGradient(Function&    theFunc,
                        math_Vector& theX,
                        math_Vector& theGrad,
                        double       theStep = 1.0e-8)
 {
   const int aLower = theX.Lower();
   const int aUpper = theX.Upper();
 
   for (int i = aLower; i <= aUpper; ++i)
   {
     const double aXi     = theX(i);
     double       aFPlus  = 0.0;
     double       aFMinus = 0.0;
 
     // Forward perturbation
     theX(i) = aXi + theStep;
     if (!theFunc.Value(theX, aFPlus))
     {
       theX(i) = aXi;
       return false;
     }
 
     // Backward perturbation
     theX(i) = aXi - theStep;
     if (!theFunc.Value(theX, aFMinus))
     {
       theX(i) = aXi;
       return false;
     }
 
     // Restore original value
     theX(i) = aXi;
 
     // Central difference
     theGrad(i) = (aFPlus - aFMinus) / (2.0 * theStep);
   }
 
   return true;
 }
 
 //! Numerical gradient with adaptive step size.
 //! Uses step proportional to |x[i]| for better conditioning.
 //!
 //! @tparam Function type with Value(const math_Vector&, double&) method
 //! @param theFunc function to differentiate
 //! @param theX point at which to evaluate gradient
 //! @param theGrad output gradient vector
 //! @param theRelStep relative step size (default 1e-8)
 //! @return true if successful
 template <typename Function>
 bool NumericalGradientAdaptive(Function&    theFunc,
                                math_Vector& theX,
                                math_Vector& theGrad,
                                double       theRelStep = 1.0e-8)
 {
   const int aLower = theX.Lower();
   const int aUpper = theX.Upper();
 
   for (int i = aLower; i <= aUpper; ++i)
   {
     const double aXi = theX(i);
     // Adaptive step: larger for larger |x|, with minimum floor
     const double aStep = theRelStep * std::max(1.0, std::abs(aXi));
 
     double aFPlus  = 0.0;
     double aFMinus = 0.0;
 
     theX(i) = aXi + aStep;
     if (!theFunc.Value(theX, aFPlus))
     {
       theX(i) = aXi;
       return false;
     }
 
     theX(i) = aXi - aStep;
     if (!theFunc.Value(theX, aFMinus))
     {
       theX(i) = aXi;
       return false;
     }
 
     theX(i)    = aXi;
     theGrad(i) = (aFPlus - aFMinus) / (2.0 * aStep);
   }
 
   return true;
 }
 
 //! Numerical Jacobian matrix for vector-valued functions.
 //! J[i,j] = dF[i]/dx[j] ~= (F[i](x + h[j]*e[j]) - F[i](x - h[j]*e[j])) / (2*h[j])
 //!
 //! @tparam Function type with Value(const math_Vector& x, math_Vector& F) method
 //! @param theFunc vector-valued function F: R^n -> R^m
 //! @param theX point at which to evaluate Jacobian (n-dimensional)
 //! @param theJac output Jacobian matrix (m x n)
 //! @param theStep step size (default 1e-8)
 //! @return true if successful
 template <typename Function>
 bool NumericalJacobian(Function&    theFunc,
                        math_Vector& theX,
                        math_Matrix& theJac,
                        double       theStep = 1.0e-8)
 {
   const int aNbRows = theJac.RowNumber();
   const int aNbCols = theJac.ColNumber();
 
   math_Vector aFPlus(1, aNbRows);
   math_Vector aFMinus(1, aNbRows);
 
   for (int j = 1; j <= aNbCols; ++j)
   {
     const int    aIdx = theX.Lower() + j - 1;
     const double aXj  = theX(aIdx);
 
     // Forward perturbation
     theX(aIdx) = aXj + theStep;
     if (!theFunc.Value(theX, aFPlus))
     {
       theX(aIdx) = aXj;
       return false;
     }
 
     // Backward perturbation
     theX(aIdx) = aXj - theStep;
     if (!theFunc.Value(theX, aFMinus))
     {
       theX(aIdx) = aXj;
       return false;
     }
 
     // Restore
     theX(aIdx) = aXj;
 
     // Fill column of Jacobian
     for (int i = 1; i <= aNbRows; ++i)
     {
       theJac(i, j) = (aFPlus(i) - aFMinus(i)) / (2.0 * theStep);
     }
   }
 
   return true;
 }
 
 //! Numerical Hessian matrix using finite differences.
 //! H[i,j] = d^2f/dx[i]dx[j]
 //! Uses central differences on gradient.
 //!
 //! @tparam Function type with Value(const math_Vector&, double&) method
 //! @param theFunc scalar function f: R^n -> R
 //! @param theX point at which to evaluate Hessian (n-dimensional)
 //! @param theHess output Hessian matrix (n x n, symmetric)
 //! @param theStep step size (default 1e-5, larger than gradient step)
 //! @return true if successful
 template <typename Function>
 bool NumericalHessian(Function&    theFunc,
                       math_Vector& theX,
                       math_Matrix& theHess,
                       double       theStep = 1.0e-5)
 {
   const int aLower = theX.Lower();
   const int aUpper = theX.Upper();
 
   double aFx = 0.0;
   if (!theFunc.Value(theX, aFx))
   {
     return false;
   }
 
   // Diagonal elements: d^2f/dx[i]^2 ~= (f(x+h[i]*e[i]) - 2f(x) + f(x-h[i]*e[i])) / h^2
   for (int i = aLower; i <= aUpper; ++i)
   {
     const double aXi     = theX(i);
     double       aFPlus  = 0.0;
     double       aFMinus = 0.0;
 
     theX(i) = aXi + theStep;
     if (!theFunc.Value(theX, aFPlus))
     {
       theX(i) = aXi;
       return false;
     }
 
     theX(i) = aXi - theStep;
     if (!theFunc.Value(theX, aFMinus))
     {
       theX(i) = aXi;
       return false;
     }
 
     theX(i) = aXi;
 
     const int aMatIdx         = i - aLower + 1;
     theHess(aMatIdx, aMatIdx) = (aFPlus - 2.0 * aFx + aFMinus) / (theStep * theStep);
   }
 
   // Off-diagonal elements: d^2f/dx[i]dx[j]
   // ~= (f(x+h[i]*e[i]+h[j]*e[j]) - f(x+h[i]*e[i]-h[j]*e[j])
   //    - f(x-h[i]*e[i]+h[j]*e[j]) + f(x-h[i]*e[i]-h[j]*e[j])) / (4h^2)
   for (int i = aLower; i <= aUpper; ++i)
   {
     for (int j = i + 1; j <= aUpper; ++j)
     {
       const double aXi  = theX(i);
       const double aXj  = theX(j);
       double       aFpp = 0.0, aFpm = 0.0, aFmp = 0.0, aFmm = 0.0;
 
       // f(x + h[i]*e[i] + h[j]*e[j])
       theX(i) = aXi + theStep;
       theX(j) = aXj + theStep;
       if (!theFunc.Value(theX, aFpp))
       {
         theX(i) = aXi;
         theX(j) = aXj;
         return false;
       }
 
       // f(x + h[i]*e[i] - h[j]*e[j])
       theX(j) = aXj - theStep;
       if (!theFunc.Value(theX, aFpm))
       {
         theX(i) = aXi;
         theX(j) = aXj;
         return false;
       }
 
       // f(x - h[i]*e[i] - h[j]*e[j])
       theX(i) = aXi - theStep;
       if (!theFunc.Value(theX, aFmm))
       {
         theX(i) = aXi;
         theX(j) = aXj;
         return false;
       }
 
       // f(x - h[i]*e[i] + h[j]*e[j])
       theX(j) = aXj + theStep;
       if (!theFunc.Value(theX, aFmp))
       {
         theX(i) = aXi;
         theX(j) = aXj;
         return false;
       }
 
       // Restore
       theX(i) = aXi;
       theX(j) = aXj;
 
       const int    aMatI = i - aLower + 1;
       const int    aMatJ = j - aLower + 1;
       const double aHij  = (aFpp - aFpm - aFmp + aFmm) / (4.0 * theStep * theStep);
 
       // Symmetric
       theHess(aMatI, aMatJ) = aHij;
       theHess(aMatJ, aMatI) = aHij;
     }
   }
 
   return true;
 }
 
 //! Second derivative using central difference.
 //! f''(x) ~= (f(x+h) - 2f(x) + f(x-h)) / h^2
 //!
 //! @tparam Function type with Value(double theX, double& theF) method
 //! @param theFunc function to differentiate
 //! @param theX point at which to evaluate second derivative
 //! @param theFx function value at theX (if already computed)
 //! @param theD2f computed second derivative value
 //! @param theStep step size (default 1e-5)
 //! @return true if successful
 template <typename Function>
 bool SecondDerivative(Function& theFunc,
                       double    theX,
                       double    theFx,
                       double&   theD2f,
                       double    theStep = 1.0e-5)
 {
   double aFPlus  = 0.0;
   double aFMinus = 0.0;
 
   if (!theFunc.Value(theX + theStep, aFPlus))
   {
     return false;
   }
   if (!theFunc.Value(theX - theStep, aFMinus))
   {
     return false;
   }
 
   theD2f = (aFPlus - 2.0 * theFx + aFMinus) / (theStep * theStep);
   return true;
 }
 
 } // namespace MathUtils
 
 #endif // _MathUtils_Deriv_HeaderFile

This page was automatically generated by the 2.3.7 LXR engine. The LXR team