//+------------------------------------------------------------------+
//|                                                     num_diff.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
//+------------------------------------------------------------------+
//| directional scheme options                                       |
//+------------------------------------------------------------------+
enum ENUM_SCHEME_DIRECTION
  {
   SCHEME_1=0,// 1-sided finite difference scheme (Forward or Backward)
   SCHEME_2   // 2-sided finite difference scheme (Central)
  };

//+------------------------------------------------------------------+
//| num points of evaluation                                         |
//+------------------------------------------------------------------+
enum ENUM_DIFF_POINTS
  {
   GRAD_POINT_2=0,     // 2-point numerical differentiation (e.g., forward difference)
   GRAD_POINT_3,       // 3-point numerical differentiation (e.g., central/stencil difference)
   GRAD_POINT_CS,      // Complex step derivative estimation method
   GRAD_POINT_CALLABLE // Exact analytical gradient is available via a callback function
  };

//+------------------------------------------------------------------+
//| num points of evaluation                                         |
//+------------------------------------------------------------------+
enum ENUM_HESS_DIFF_POINTS
  {
   HESS_POINT_2=0,           // 2-point estimation for Hessian calculation
   HESS_POINT_3,             // 3-point estimation for Hessian calculation
   HESS_POINT_CS,            // Complex step differentiation for Hessian estimation
   HESS_POINT_HESS_STRATEGY, // Quasi-Newton Hessian update strategy (e.g., BFGS, L-BFGS)
   HESS_POINT_CALLABLE       // Exact analytical Hessian is available via a callback function
  };

//+------------------------------------------------------------------+
//| struct objective function return                                 |
//+------------------------------------------------------------------+
struct ObjReturn
  {
   double            f;  // Scalar value of the objective function: f(x)
   vector            g;  // Gradient vector: df/dx
   vector            mf; // Multi-valued objective function return (for vector-valued functions)
   matrix            mg; // Jacobian matrix / Multi-valued gradient response

   // Default constructor initializing structures to empty/zero values
   ObjReturn(void)
     {
      f = double(0);
      g = vector::Zeros(0);
      mf = vector::Zeros(0);
      mg = matrix::Zeros(0,0);
     }
   
   // Copy constructor
   ObjReturn(ObjReturn& other)
     {
      f = other.f;
      g = other.g;
      mf = other.mf;
      mg = other.mg;
     }
   
   // Assignment operator overload
   void             operator=(ObjReturn& other)
     {
      f = other.f;
      g = other.g;
      mf = other.mf;
      mg = other.mg;
     }
  };

//+------------------------------------------------------------------+
//| IObjective provides the base interface for any function          |
//| that will be provided to a minimizer routine                     |
//+------------------------------------------------------------------+
interface IObjective
  {
   // Evaluates and returns the vector-valued response of the objective function at point x
   vector objective_function(vector& x);
   
   // Evaluates both the scalar objective value and its gradient at point x
   ObjReturn fun_and_grad(vector& x);
  };

//+------------------------------------------------------------------+
//| differentiation options                                           |
//+------------------------------------------------------------------+
struct GradDiffOptions
  {
   ENUM_DIFF_POINTS method;   // Numerical differentiation technique selection
   vector           rel_step; // Relative step size vector for finite differences
   vector           abs_step; // Absolute step size vector for finite differences
   matrix           bounds;   // Boundary constraints matrix [N x 2] containing [lower, upper] limits

   // Default constructor initializing variables to invalid or empty states
   GradDiffOptions(void)
     {
      method = WRONG_VALUE;
      rel_step = abs_step = vector::Zeros(0);
      bounds = matrix::Zeros(0,0);
     }
   
   // Parameterized constructor
   GradDiffOptions(ENUM_DIFF_POINTS m, vector& relstep, vector& absstep, matrix& bnds)
     {
      method = m;
      rel_step = relstep;
      abs_step = absstep;
      bounds = bnds;
     }
   
   // Copy constructor
   GradDiffOptions(GradDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      bounds = other.bounds;
     }
   
   // Assignment operator overload
   void             operator=(GradDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      bounds = other.bounds;
     }
  };

//+------------------------------------------------------------------+
//| differentiation options                                           |
//+------------------------------------------------------------------+
struct HessDiffOptions
  {
   bool                  as_linear_operator; // Flag to treat Hessian operations as matrix-free linear operator mappings
   ENUM_HESS_DIFF_POINTS method;             // Hessian evaluation method selection
   vector                rel_step;           // Relative step sizes for Hessian finite differences
   vector                abs_step;           // Absolute step sizes for Hessian finite differences

   // Default constructor
   HessDiffOptions(void)
     {
      method = WRONG_VALUE;
      as_linear_operator = false;
      rel_step = abs_step = vector::Zeros(0);
     }
   
   // Parameterized constructor
   HessDiffOptions(ENUM_HESS_DIFF_POINTS m, vector& relstep, vector& absstep, bool aslinearoperator)
     {
      method = m;
      rel_step = relstep;
      abs_step = absstep;
      as_linear_operator = aslinearoperator;
     }
   
   // Copy constructor
   HessDiffOptions(HessDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      as_linear_operator = other.as_linear_operator;
     }
   
   // Assignment operator overload
   void             operator=(HessDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      as_linear_operator = other.as_linear_operator;
     }
  };

//+---------------------------------------------------------------------------+
//| function objective representing the objective function and its derivatives.|
//+---------------------------------------------------------------------------+
class CFunctor:public IObjective
  {
protected:
   vector            m_xp;              // Previous execution parameters state vector
   vector            m_x;               // Current state variable parameters vector
   ulong             m_n;               // Dimensionality size of the parameter vector
   matrix            m_H;               // System Hessian matrix
   int               m_nfev,m_ngev,m_nhev; // Evaluation counters for function, gradient, and Hessian evaluations
   bool              m_fupdated,m_gupdated,m_hupdated; // State dirty flags tracking whether properties require calculation
   double            m_lowest_f,m_f;    // Globally tracked minimum objective value discovered, and current value
   vector            m_lowest_x,m_g;    // Variable state coordinates yielding m_lowest_f, and current gradient vector
   bool              m_clip;            // Enforcement flag to bound variables before function calculations
   GradDiffOptions   m_grad_options;    // Operational configuration settings for gradients
   HessDiffOptions   m_hess_options;    // Operational configuration settings for Hessians
   
   // Internal method to compute and update the objective function value if state changed
   void              update_fun(void)
     {
      if(!m_fupdated)
        {
         double fx = wrapped_fun(m_x);
         if(fx<m_lowest_f)
           {
            m_lowest_f = fx;
            m_lowest_x = m_x;
           }
         m_f  = fx;
         m_fupdated = true;
        }
     }
   
   // Internal method to calculate and refresh gradient vectors if outdated
   void              update_grad(void)
     {
      if(!m_gupdated)
        {
         if(m_grad_options.method!=GRAD_POINT_CALLABLE)
            update_fun();
         vector ff(1);
         ff[0] = m_f;
         m_g = wrapped_grad(m_x,ff);
         m_gupdated = true;
        }
     }
   
   // Internal method to compute and refresh system Hessian matrices if outdated
   void              update_hess(void)
     {
      if(!m_hupdated)
        {
         if(m_hess_options.method != HESS_POINT_CALLABLE)
           {
            update_grad();
            m_H = wrapped_hess(m_x,m_g);
           }
         else
           {
            vector a = vector::Zeros(0);
            m_H = wrapped_hess(m_x,a);
           }
         m_hupdated = true;
        }
     }
   
   // Internal helper method to reset states when location coordinates shift
   void              update_x(vector& x)
     {
      m_x = x;
      m_fupdated = m_hupdated = m_gupdated = false; // Flag matrices as requiring computation updates
     }

public:
                     CFunctor(void)
     {
      m_fupdated = m_hupdated = m_gupdated = false;
      m_lowest_f = DBL_MAX;
      m_lowest_x = m_g = vector::Zeros(0);
      m_H = matrix::Zeros(0,0);
      m_nfev = m_ngev = m_nhev = 0;
      m_grad_options.method = GRAD_POINT_2;
      m_hess_options.method = HESS_POINT_CALLABLE;
      m_hess_options.as_linear_operator = true;
      m_clip = false;
     }
                    ~CFunctor(void)
     {
     }
   
   // Set configuration structures and boundary parameter limits
   void              setParams(GradDiffOptions& grad_opts, bool clip_to_bounds=true)
     {
      m_clip = clip_to_bounds;
      m_grad_options = grad_opts;
     }
   
   // Toggles out of bounds variable clipping strategy
   void                   setClipOption(bool yes)
     {
      m_clip = yes;
     }
   
   // Strategy implementation method configurer for gradients
   void                   setGradOption(ENUM_DIFF_POINTS grad)
     {
      m_grad_options.method = grad;
     }
   
   // Set constant numerical step sizing for delta changes
   void                    setAbsoluteStep(vector& epsilon)
     {
      m_grad_options.abs_step = epsilon;
      m_hess_options.abs_step = epsilon;
     }
   
   // Set physical box-constraints configuration matrix boundaries
   void                    setBounds(matrix& finite_bounds)
     {
      m_grad_options.bounds = finite_bounds;
     }
   
   // Set scalar relative ratio steps for step assessments
   void                    setRelativeStep(vector& finite_diff_rel_step)
     {
      m_grad_options.rel_step = finite_diff_rel_step;
      m_hess_options.rel_step = finite_diff_rel_step;
     }
   
   // Strategy implementation method configurer for Hessians
   void                   setHessOption(ENUM_HESS_DIFF_POINTS hess)
     {
      m_hess_options.method = hess;
     }
   
   // Validates configuration dimensions, tests math values, ensures solver constraints safety 
   bool                    initialize(vector& x)
     {
      m_x = x;
      m_xp = m_x;
      m_n = x.Size();
      
      // Ensure logic integrity: dynamic Hessian approximation requires at least a known analytical gradient
      if(m_grad_options.method != GRAD_POINT_CALLABLE && m_hess_options.method != HESS_POINT_CALLABLE)
        {
         Print(__FUNCTION__, "Whenever the gradient is estimated via "
               "finite-differences, it is required that"
               " the Hessian "
               "be estimated using one of the "
               "quasi-Newton strategies.");
         return false;
        }

      // Check objective function initial return sanity
      double check = orig_fun(m_x);
      if(MathClassify(check)!=FP_NORMAL)
        {
         Print(__FUNCTION__," check the implementation of the objective function, currently evaluates to an invalid number ");
         return false;
        }

      // Check analytical gradient dimensions if used
      if(m_grad_options.method == GRAD_POINT_CALLABLE)
        {
         vector a  = grad_fun(m_x);
         if(!a.Size())
           {
            Print(__FUNCTION__, " check the implementation of the overriden gradient function, currently evaluates to an empty vector ");
            return false;
           }
        }

      // Correct and populate boundary structural constraints fallback defaults if size mismatch occurs
      if(m_grad_options.bounds.Rows()!=ulong(m_n))
        {
         m_grad_options.bounds.Resize(m_n,2);
         m_grad_options.bounds.Fill(double("inf"));
         for(ulong i  = 0; i<m_n; ++i)
            m_grad_options.bounds[i,0]*=-1.0;
         Print(__FUNCTION__," ","WARNING:Boundary constraints not properly specified.\nOverwritten to (-inf,inf) ");
        }

      update_fun();
      update_grad();

      if(m_hess_options.method == HESS_POINT_CALLABLE)
        {
         vector a = vector::Zeros(0);
         m_H = wrapped_hess(x,a);
         m_hupdated = true;
        }

      return true;
     }

   // Encapsulates objective evaluation loop, handles safety box limits tracking, counts execution passes
   double                  wrapped_fun(vector& x)
     {
      m_nfev += 1;
      if(m_clip)
        {
         matrix bnds = m_grad_options.bounds;
         for(ulong i = 0; i<x.Size(); ++i)
           {
            if(x[i]<bnds[i,0])
               x[i] = bnds[i,0];
            else
               if(x[i]>bnds[i,1])
                  x[i] = bnds[i,1];
           }
        }
      vector copy = x;
      return orig_fun(copy);
     }
   
   // Interface compliance proxy method accessing calculation engines
   vector             objective_function(vector& x)
     {
      vector r(1);
      r[0] = wrapped_fun(x);
      return r;
     }
   
   // Evaluates analytical solutions or forwards routines into numerical step approximations for gradients
   vector             wrapped_grad(vector& x,vector& f0)
     {
      m_ngev += 1;
      if(m_clip)
        {
         matrix bnds = m_grad_options.bounds;
         for(ulong i = 0; i<x.Size(); ++i)
           {
            if(x[i]<bnds[i,0])
               x[i] = bnds[i,0];
            else
               if(x[i]>bnds[i,1])
                  x[i] = bnds[i,1];
           }
        }
      vector copy = x;
      if(m_grad_options.method == GRAD_POINT_CALLABLE)
         return grad_fun(copy);
      IObjective* objective = GetPointer(this);
      matrix ad = approx_derivative(objective,copy,f0,m_grad_options.method,m_grad_options.rel_step,m_grad_options.abs_step,m_grad_options.bounds);
      return ad.Row(0);
     }
   
   // Evaluates analytical solutions or forwards routines into numerical step approximations for Hessians
   matrix             wrapped_hess(vector& x, vector& f0)
     {
      m_nhev += 1;
      if(m_clip)
        {
         matrix bnds = m_grad_options.bounds;
         for(ulong i = 0; i<x.Size(); ++i)
           {
            if(x[i]<bnds[i,0])
               x[i] = bnds[i,0];
            else
               if(x[i]>bnds[i,1])
                  x[i] = bnds[i,1];
           }
        }
      vector copy = x;
      if(m_hess_options.method == HESS_POINT_CALLABLE)
         return hess_fun(copy);
      IObjective* objective = GetPointer(this);
      return approx_derivative(objective,x,f0,m_grad_options.method,m_grad_options.rel_step,m_grad_options.abs_step,m_grad_options.bounds);
     }
   
   // Data parsing helpers extracting specific constraints vectors
   vector             lower_bounds(void)  { return m_grad_options.bounds.Col(0); }
   vector             upper_bounds(void)  { return m_grad_options.bounds.Col(1); }
   vector             initial_params(void) { return m_xp; }

   // Virtual placeholder declarations designed to be overriden by inherited custom user targets
   virtual double    orig_fun(vector& x) { return double("nan"); }
   virtual vector    grad_fun(vector& x) { return vector::Zeros(0); }
   virtual matrix    hess_fun(vector& x) { return matrix::Zeros(0,0); }
   
   // Consolidated extraction layer getting simultaneously computed coordinates pairs
   ObjReturn         fun_and_grad(vector& x)
     {
      vector dif  = MathAbs(m_x - x);
      if(dif.Sum() >= DBL_EPSILON || dif.HasNan())
         update_x(x);
      update_fun();
      update_grad();
      ObjReturn out;
      out.f = m_f;
      out.g = m_g;
      return out;
     }
  };

//+------------------------------------------------------------------+
//| base class representing function and its derivative               |
//+------------------------------------------------------------------+
class CConstraints:public IObjective
  {
protected:
   int               m_n,m_m;          // m_n: parameter dimension space size, m_m: total quantity of constraint equations
   bool              m_clip,m_initialized; // Status switches managing operation modes and validation configurations
   GradDiffOptions   m_options;        // Settings configuration defining differentiation procedures
   
   // Evaluates constraint functions wrapping execution boundaries securely
   vector            wrapped_func(vector& x)
     {
      if(!m_initialized)
         return vector::Zeros(0);
      if(m_clip)
        {
         matrix bnds = m_options.bounds;
         for(ulong i = 0; i<x.Size(); ++i)
           {
            if(x[i]<bnds[i,0])
               x[i] = bnds[i,0];
            else
               if(x[i]>bnds[i,1])
                  x[i] = bnds[i,1];
           }
        }
      vector copy = x;
      return orig_fun(copy);
     }
   
   // Computes the constraints Jacobian matrix (numerical or analytical derivations)
   matrix            wrapped_grad(vector& x,vector& f0)
     {
      if(!m_initialized)
         return matrix::Zeros(0,0);
      if(m_clip)
        {
         matrix bnds = m_options.bounds;
         for(ulong i = 0; i<x.Size(); ++i)
           {
            if(x[i]<bnds[i,0])
               x[i] = bnds[i,0];
            else
               if(x[i]>bnds[i,1])
                  x[i] = bnds[i,1];
           }
        }
      vector copy = x;
      if(m_options.method == GRAD_POINT_CALLABLE)
         return orig_grad(copy);
      IObjective* fun = (IObjective*)GetPointer(this);
      matrix ad = approx_derivative(fun,copy,f0,m_options.method,m_options.rel_step,m_options.abs_step,m_options.bounds);
      return ad;
     }
     
public:
                     CConstraints(void)
     {
      m_n = m_m = 0;
      m_clip = true;
      m_initialized = false;
      m_options.method = GRAD_POINT_2;
     }
                    ~CConstraints(void)
     {
     }
   
   // Validates systems shapes and initializes structural memory arrays
   bool                    initialize(vector& x)
     {
      if((m_m == 0 && m_n) || (int)x.Size()!=m_n)
        {
         Print(__FUNCTION__," ","Invalid property dimensions.", " n ", m_n, " vs ", x.Size());
         return false;
        }

      vector check = orig_fun(x);
      if(!check.Size() || check.HasNan() || MathClassify(fabs(check.Max())) == FP_INFINITE)
        {
         Print(__FUNCTION__," ","Check the implementation of the constraints function.");
         return false;
        }

      if(m_options.method == GRAD_POINT_CALLABLE)
        {
         matrix a  = orig_grad(x);
         if(a.Rows()!=m_m || a.HasNan() || MathClassify(fabs(a.Max())) == FP_INFINITE)
           {
            Print(__FUNCTION__," ","Check the implementation of the overriden constraints gradient function.");
            return false;
           }
        }
      if(m_options.bounds.Rows()!=ulong(m_n))
        {
         m_options.bounds.Resize(m_n,2);
         m_options.bounds.Fill(double("inf"));
         for(int i  = 0; i<m_n; ++i)
            m_options.bounds[i,0]*=-1.0;
         Print(__FUNCTION__, " ","WARNING:Boundary constraints not properly specified.\nOverwritten to (-inf,inf) ");
        }
      m_initialized = true;
      return m_initialized;
     }

   // Setter methods for custom properties management overrides
   void              setParams(int m, int n, GradDiffOptions& opts, bool clip_to_bounds=true)
     {
      m_initialized = false;
      m_n = fabs(n); m_m = fabs(m); m_clip = clip_to_bounds; m_options = opts;
     }
   void              setGradientOptions(GradDiffOptions& opts) { m_initialized = false; m_options = opts; }
   void              setBounds(matrix& bounds)                { m_initialized = false; m_options.bounds = bounds; }
   void              setNumConstraints(int m)                 { m_initialized = false; m_m = m; }
   void              setDimension(int n)                      { m_initialized = false; m_n = n; }

   // Proxy redirect interfaces fulfilling generic abstract execution layouts
   vector             objective_function(vector& x)           { return wrapped_func(x); }
   
   // Dual mapping executor packing multivariant array packages sequentially
   ObjReturn         fun_and_grad(vector& x)
     {
      ObjReturn out;
      out.mf = objective_function(x);
      out.mg = wrapped_grad(x,out.mf);
      return out;
     }
   
   // Base virtual targets intended for framework adaptation usage expansions
   virtual vector    orig_fun(vector& x)  { return vector::Zeros(0); }
   virtual matrix    orig_grad(vector& x) { return matrix::Zeros(0,0); }
   
   int               numConstraints(void) { return m_m; }
   int               dim(void)            { return m_n; }
   vector            lower_bounds(void)   { return m_options.bounds.Col(0); }
   vector            upper_bounds(void)   { return m_options.bounds.Col(1); }
  };

//+----------------------------------------------------------------------------------+
//| Calculates relative EPS step to use for a given data type and numdiff step method.|
//+----------------------------------------------------------------------------------+
double eps_for_method(ENUM_DIFF_POINTS method)
  {
   switch(method)
     {
      case GRAD_POINT_2:
      case GRAD_POINT_CS:
         // Optimal step for 2-point central/forward schema is roughly sqrt(machine_epsilon)
         return pow(2.220446049250313e-16,0.5);
      case GRAD_POINT_3:
         // Optimal step for 3-point formulations is roughly machine_epsilon^(1/3)
         return pow(2.220446049250313e-16,(1./3.));
     };
   return DBL_EPSILON;
  }

//+---------------------------------------------------------------------------------+
//| Computes an absolute step from a relative step for finite difference calculation.|
//+---------------------------------------------------------------------------------+
vector compute_absolute_step(vector& rel_step,vector& x0, vector& f0, ENUM_DIFF_POINTS method)
  {
   vector signx0 = x0;
   vector abs_step = signx0;
   
   // Extracts directional polarity signs across variables space arrays
   for(ulong i = 0; i<signx0.Size(); ++i)
     {
      if(x0[i] >= 0.)
         signx0[i]=1.*2-1; // Yields 1.0
      else
         signx0[i] = 0.0*2-1; // Yields -1.0
     }
   double rstep = eps_for_method(method);

   // Automatically generates fallback absolute scaling matrices if relative instructions are missing
   if(rel_step.Size()==0)
      for(ulong i = 0; i<abs_step.Size(); ++i)
         abs_step[i] = rstep*signx0[i]*MathMax(1.,fabs(x0[i]));
   else
     {
      abs_step = rstep*signx0*MathAbs(x0);
      vector dx = ((x0+abs_step) - x0);
      for(ulong i = 0; i<abs_step.Size(); ++i)
         if(dx[i] == 0.0) // Handles edge cases near exactly zero points
            abs_step[i] = rstep*signx0[i]*MathMax(1.,fabs(x0[i]));
     }
   return abs_step;
  }

//+------------------------------------------------------------------+
//| Adjust finite difference scheme to the presence of bounds.       |
//+------------------------------------------------------------------+
vector adjust_scheme_to_bounds(vector& x0, vector& h, int num_steps,ENUM_SCHEME_DIRECTION scheme, vector& lb, vector& ub, vector &one_sided)
  {
   switch(scheme)
     {
      case SCHEME_1:
         one_sided = vector::Ones(h.Size());
         break;
      case SCHEME_2:
         one_sided = vector::Ones(h.Size());
         h = MathAbs(h);
         break;
     }

   // Skip evaluations entirely if boundary limits are set to infinite domains
   bool all_true = true;
   for(ulong i = 0; i<x0.Size(); ++i)
      if(lb[i] != -double("inf") || ub[i] != double("inf"))
        {
         all_true = false;
         break;
        }

   if(all_true)
      return h;

   vector h_total  = h * double(num_steps);
   vector h_adjusted = h;
   vector lower_dist = x0 - lb;
   vector upper_dist = ub - x0;

   int forward,backward,fitting,violated,central, adjusted_central;
   forward = backward = violated = fitting = central  = false;
   double x = 0.;
   double min_dist = 0.;
   
   switch(scheme)
     {
      case SCHEME_1: // Adjustments for 1-sided scheme (swaps sign if upper limit hit)
        {
         for(ulong i = 0; i<h.Size(); ++i)
           {
            x = x0[i] + h_total[i];
            violated = int(x<lb[i]|x>ub[i]);
            fitting  = int(fabs(h_total[i])<=MathMax(lower_dist[i],upper_dist[i]));
            if(bool(violated & fitting))
               h_adjusted[i]*=-1.; // Flip the direction to stay within bounds
            forward = int((upper_dist[i] >= lower_dist[i]) & ~fitting);
            if(forward)
               h_adjusted[i] = upper_dist[i]/double(num_steps);
            backward = int((upper_dist[i]<lower_dist[i]) & ~fitting);
            if(backward)
               h_adjusted[i] = -lower_dist[i]/double(num_steps);
           }
        }
      break;
      case SCHEME_2: // Adjustments for 2-sided scheme (fallback to one-sided if bound is close)
        {
         for(ulong i = 0; i<h.Size(); ++i)
           {
            central = int(((lower_dist[i]>=h_total[i]) & (upper_dist[i] >= h_total[i])));
            forward = int(((upper_dist[i]>=lower_dist[i]) & ~central));
            if(forward)
              {
               h_adjusted[i] = MathMin(h[i],0.5*upper_dist[i]/double(num_steps));
               one_sided[i] = 1.; // Enforce transformation shift to 1-sided estimation layout
              }
            backward = int(((upper_dist[i]<lower_dist[i]) & ~central));
            if(backward)
              {
               h_adjusted[i] = -1.* MathMin(h[i],0.5*lower_dist[i]/double(num_steps));
               one_sided[i] = 1.0;
              }
            min_dist = MathMin(upper_dist[i],lower_dist[i])/double(num_steps);
            adjusted_central = int((~central & (fabs(h_adjusted[i])<=min_dist)));
            if(adjusted_central)
              {
               h_adjusted[i] = min_dist;
               one_sided[i] = 0.;
              }
           }
        }
      break;
     }
   return h_adjusted;
  }

//+------------------------------------------------------------------+
//| dense difference                                                 |
//+------------------------------------------------------------------+
matrix dense_difference(IObjective* fun, vector& x0, vector& f0, vector& h, vector& use_one_sided, ENUM_DIFF_POINTS method)
  {
   ulong m = f0.Size();
   ulong n = x0.Size();

   matrix j_transposed = matrix::Zeros(n,m); // Target transpose container matrix storage

   vector x1 = x0;
   vector x2 = x0;
   vector df = vector::Zeros(f0.Size());

   // Iteratively loops through parameters dimension vectors evaluating stencil formulas
   for(ulong i = 0; i<h.Size(); ++i)
     {
      double dx = 1.e-12;
      if(method == GRAD_POINT_2)
        {
         x1[i] += h[i];
         dx = x1[i] - x0[i];
         df = fun.objective_function(x1) - f0; // Forward formula df = f(x+h) - f(x)
        }
      else
         if(method == GRAD_POINT_3 && use_one_sided[i]!=0.0)
           {
            x1[i] += h[i];
            x2[i] += 2. * h[i];
            dx = x2[i] - x0[i];
            df = -3.0 * f0 + 4 * fun.objective_function(x1) - fun.objective_function(x2); // One-sided 3-point formula
           }
         else
            if(method == GRAD_POINT_3 && use_one_sided[i]==0.0)
              {
               x1[i] -= h[i];
               x2[i] += h[i];
               dx = x2[i] - x1[i];
               df = fun.objective_function(x2) - fun.objective_function(x1); // Central 3-point formula
              }

      j_transposed.Row(df/dx,i); // Set finite approximation row slice output mappings
      x1[i] = x2[i] = x0[i];     // Resets working indexes trackers state
     }
   return j_transposed.Transpose();
  }

//+---------------------------------------------------------------------------------------+
//| Compute finite difference approximation of the derivatives of a vector-valued function. |
//+---------------------------------------------------------------------------------------+
matrix approx_derivative(IObjective* fun,vector& x0,vector &f0,ENUM_DIFF_POINTS method, vector &rel_step,vector &abs_step, matrix& bounds)
  {
   // Safety pointer execution health check verification
   if(CheckPointer(fun)==POINTER_INVALID)
     {
      Print(__FUNCTION__, " fun variable is an invalid pointer ");
      return matrix::Zeros(0,0);
     }

   vector lb,ub;
   lb = bounds.Col(0);
   ub = bounds.Col(1);

   if(lb.Size()!=x0.Size() ||ub.Size()!=x0.Size())
     {
      Print(__FUNCTION__, " inconsistent shapes between bounds and x0 ");
      return matrix::Zeros(0,0);
     }

   // Ensure structural function values exist
   if(!f0.Size())
      f0 = fun.objective_function(x0);
      
   // Bounds containment sanity validation checks
   for(ulong i = 0; i<x0.Size(); ++i)
      if(x0[i]<lb[i] || x0[i]>ub[i])
        {
         Print(__FUNCTION__, " x0 violates bound constraints ");
         return matrix::Zeros(0,0);
        }

   // Step evaluation logic execution passes
   vector h;
   if(!abs_step.Size())
      h = compute_absolute_step(rel_step,x0,f0,method);
   else
     {
      h = abs_step;
      vector signx0 = vector::Zeros(x0.Size());
      for(ulong i = 0; i<x0.Size(); ++i)
        {
         if(x0[i]>=0.0)
            signx0[i] = 1.0*2.-1.;
         else
            signx0[i] = 0.0*2.-1.;
         if(((x0[i]+h[i]) - x0[i]) == 0.0)
            h[i] = eps_for_method(method)*signx0[i]*MathMax(1.,fabs(x0[i]));
        }
     }

   // Schemes alignment corrections based on active borders settings
   vector use_one_sided;
   switch(method)
     {
      case GRAD_POINT_2:
         h = adjust_scheme_to_bounds(x0,h,1,SCHEME_1,lb,ub,use_one_sided);
         break;
      case GRAD_POINT_3:
         h = adjust_scheme_to_bounds(x0,h,1,SCHEME_2,lb,ub,use_one_sided);
         break;
      case GRAD_POINT_CS:
         use_one_sided = vector::Zeros(x0.Size());
         break;
     }
     
   // Triggers execution calculation pipelines over defined matrices grids
   return dense_difference(fun,x0,f0,h,use_one_sided,method);
  }
//+------------------------------------------------------------------+