//+------------------------------------------------------------------+
//|                                                     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
   SCHEME_2//2 sided
  };
//+------------------------------------------------------------------+
//| num points of evaluation                                         |
//+------------------------------------------------------------------+
enum ENUM_DIFF_POINTS
  {
   GRAD_POINT_2=0,//2-point
   GRAD_POINT_3,//3-point
   GRAD_POINT_CS,//complex
   GRAD_POINT_CALLABLE//callable
  };
//+------------------------------------------------------------------+
//| num points of evaluation                                         |
//+------------------------------------------------------------------+
enum ENUM_HESS_DIFF_POINTS
  {
   HESS_POINT_2=0,//2-point
   HESS_POINT_3,//3-point
   HESS_POINT_CS,//complex
   HESS_POINT_HESS_STRATEGY,//hessian update strategy
   HESS_POINT_CALLABLE//callable
  };
//+------------------------------------------------------------------+
//|struct objective function return                                  |
//+------------------------------------------------------------------+
struct ObjReturn
  {
   double            f;
   vector            g;
   vector            mf;
   matrix            mg;
                     ObjReturn(void)
     {
      f = double(0);
      g = vector::Zeros(0);
      mf = vector::Zeros(0);
      mg = matrix::Zeros(0,0);
     }
                     ObjReturn(ObjReturn& other)
     {
      f = other.f;
      g = other.g;
      mf = other.mf;
      mg = other.mg;
     }
   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
  {
//---the objective function
   vector objective_function(vector& x);
   ObjReturn fun_and_grad(vector& x);
  };
//+------------------------------------------------------------------+
//|differentiation options                                           |
//+------------------------------------------------------------------+
struct GradDiffOptions
  {
   ENUM_DIFF_POINTS  method;
   vector            rel_step;
   vector            abs_step;
   matrix            bounds;

                     GradDiffOptions(void)
     {
      method = WRONG_VALUE;
      rel_step = abs_step = vector::Zeros(0);
      bounds = matrix::Zeros(0,0);
     }
                     GradDiffOptions(ENUM_DIFF_POINTS m, vector& relstep, vector& absstep, matrix& bnds)
     {
      method = m;
      rel_step = relstep;
      abs_step = absstep;
      bounds = bnds;
     }
                     GradDiffOptions(GradDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      bounds = other.bounds;
     }
   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;
   ENUM_HESS_DIFF_POINTS method;
   vector            rel_step;
   vector            abs_step;

                     HessDiffOptions(void)
     {
      method = WRONG_VALUE;
      as_linear_operator = false;
      rel_step = abs_step = vector::Zeros(0);
     }
                     HessDiffOptions(ENUM_HESS_DIFF_POINTS m, vector& relstep, vector& absstep, bool aslinearoperator)
     {
      method = m;
      rel_step = relstep;
      abs_step = absstep;
      as_linear_operator = aslinearoperator;
     }
                     HessDiffOptions(HessDiffOptions& other)
     {
      method = other.method;
      rel_step = other.rel_step;
      abs_step = other.abs_step;
      as_linear_operator = other.as_linear_operator;
     }
   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;
   vector            m_x;
   ulong             m_n;
   matrix            m_H;
   int               m_nfev,m_ngev,m_nhev;
   bool              m_fupdated,m_gupdated,m_hupdated;
   double            m_lowest_f,m_f;
   vector            m_lowest_x,m_g;
   bool              m_clip;
   GradDiffOptions   m_grad_options;
   HessDiffOptions   m_hess_options;
   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;
        }
     }
   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;
        }
     }
   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;
        }
     }
   void              update_x(vector& x)
     {
      m_x = x;
      m_fupdated = m_hupdated = m_gupdated = false;
     }

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)
     {
     }
   void              setParams(GradDiffOptions& grad_opts, bool clip_to_bounds=true)
     {
      m_clip = clip_to_bounds;
      m_grad_options = grad_opts;
     }
   void                   setClipOption(bool yes)
     {
      m_clip = yes;
     }
   void                   setGradOption(ENUM_DIFF_POINTS grad)
     {
      m_grad_options.method = grad;
     }
   void                    setAbsoluteStep(vector& epsilon)
     {
      m_grad_options.abs_step = epsilon;
      m_hess_options.abs_step = epsilon;
     }
   void                   setBounds(matrix& finite_bounds)
     {
      m_grad_options.bounds = finite_bounds;
     }
   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;
     }
   void                   setHessOption(ENUM_HESS_DIFF_POINTS hess)
     {
      m_hess_options.method = hess;
     }
   bool                   initialize(vector& x)
     {
      m_x = x;
      m_xp = m_x;
      m_n = x.Size();
      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;
        }

      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;
        }

      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;
           }
        }

      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;
     }

   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);
     }
   vector            objective_function(vector& x)
     {
      vector r(1);
      r[0] = wrapped_fun(x);
      return r;
     }
   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);
     }
   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);
     }
   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 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);
     }
   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;
   bool              m_clip,m_initialized;
   GradDiffOptions   m_options;
   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);
     }
   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)
     {
     }
   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;
     }

   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;
     }

   vector            objective_function(vector& x)
     {

      return wrapped_func(x);
     }

   ObjReturn         fun_and_grad(vector& x)
     {
      ObjReturn out;
      out.mf = objective_function(x);
      out.mg = wrapped_grad(x,out.mf);
      return out;
     }
   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:
         return pow(2.220446049250313e-16,0.5);
      case GRAD_POINT_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;
   for(ulong i = 0; i<signx0.Size(); ++i)
     {
      if(x0[i] >= 0.)
         signx0[i]=1.*2-1;
      else
         signx0[i] = 0.0*2-1;
     }
   double rstep = eps_for_method(method);

   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)
            abs_step[i] = rstep*signx0[i]*MathMax(1.,fabs(x0[i]));
     }
   return abs_step;
  }
//+------------------------------------------------------------------+
//|Adjust final 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;
     }

   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:
        {
         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.;
            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:
        {
         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.;
              }
            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);

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

   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;
        }
      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);
           }
         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);
              }

      j_transposed.Row(df/dx,i);
      x1[i] = x2[i] = x0[i];
     }
   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/*sparsity,as linear_operator*/)
  {

   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 shaptes between bounds and x0 ");
      return matrix::Zeros(0,0);
     }

   if(!f0.Size())
      f0 = fun.objective_function(x0);
   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);
        }

   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]));
        }
     }

   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;
     }
   return dense_difference(fun,x0,f0,h,use_one_sided,method);
  }

//+------------------------------------------------------------------+