875 lines
31 KiB
MQL5
875 lines
31 KiB
MQL5
|
//+------------------------------------------------------------------+
|
||
|
|
//| 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);
|
||
|
|
}
|
||
|
|
//+------------------------------------------------------------------+
|