836 righe
25 KiB
MQL5
836 righe
25 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
|
||
|
|
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);
|
||
|
|
}
|
||
|
|
|
||
|
|
//+------------------------------------------------------------------+
|