gavindiaz/include
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
amirghadiri1987
2025-02-07 19:12:13 +03:30
committed by GitHub
parent a9340e5d4e
commit 56b382609f
88 changed files with 274630 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Math/Alglib/alglib.mqh
+38345
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/alglibinternal.mqh
+15934
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/alglibmisc.mqh
+2575
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/ap.mqh
+2549
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/arrayresize.mqh
+104
View File
@@ -0,0 +1,104 @@
//+------------------------------------------------------------------+
//| arrayresize.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
//--- forward declaration
class CRowInt;
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(int &arr[],const int size)
{
int old=ArraySize(arr);
int res=ArrayResize(arr,size);
//--- fill array if necessary
if(res>0 && old<size)
ArrayFill(arr,old,size-old,0);
//--- return result
return(res);
}
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(short &arr[],const int size)
{
int old=ArraySize(arr);
int res=ArrayResize(arr,size);
//--- fill array if necessary
if(res>0 && old<size)
ArrayFill(arr,old,size-old,0);
//--- return result
return(res);
}
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(char &arr[],const int size)
{
int old=ArraySize(arr);
int res=ArrayResize(arr,size);
//--- fill array if necessary
if(res>0 && old<size)
ArrayFill(arr,old,size-old,0);
//--- return result
return(res);
}
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(bool &arr[],const int size)
{
int old=ArraySize(arr);
int res=ArrayResize(arr,size);
//--- fill array if necessary
if(res>0 && old<size)
ArrayFill(arr,old,size-old,false);
//--- return result
return(res);
}
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(string &arr[],const int size)
{
return(ArrayResize(arr,size));
}
//+------------------------------------------------------------------+
//| ArrayResizeAL for Alglib library with MQL4 and MQL5 features |
//+------------------------------------------------------------------+
int ArrayResizeAL(CRowInt &arr[],const int size)
{
return(ArrayResize(arr,size));
}
//+------------------------------------------------------------------+
Math/Alglib/bitconvert.mqh
+411
View File
@@ -0,0 +1,411 @@
//+------------------------------------------------------------------+
//| bitconvert.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
#include "arrayresize.mqh"
//+------------------------------------------------------------------+
//| Converting numbers into an array of bits and vice versa |
//+------------------------------------------------------------------+
class BitConverter
{
public:
static void GetBytes(const int d,uchar &bytes[]);
static void GetBytes(const double d,uchar &bytes[]);
static int ToInt32(uchar &bytes[]);
static double ToDouble(uchar &bytes[]);
static bool IsLittleEndian(void);
};
//+------------------------------------------------------------------+
//| Converting integer to a byte array |
//+------------------------------------------------------------------+
void BitConverter::GetBytes(const int d,uchar &bytes[])
{
//--- create variables
int x;
int q;
int r;
int i;
int div;
//--- allocation
ArrayResize(bytes,4);
for(i=0; i<4; i++)
{
//--- check
if(d>=0)
bytes[i]=0;
else
bytes[i]=255;
}
//--- initialization
q=-1;
r=-1;
i=3;
div=256*256*256;
//--- check
if(d<0)
x=~d;
else
x=d;
//--- converting number
while(i!=-1)
{
//--- quotient
q=x/div;
//--- remainder of division
r=x%div;
//--- get byte
if(d>=0)
bytes[i]+=(uchar)q;
else
bytes[i]-=(uchar)q;
//--- the next iteration is reduced divisor
x=r;
div=div/256;
i--;
}
}
//+------------------------------------------------------------------+
//| Converting double to a byte array |
//+------------------------------------------------------------------+
void BitConverter::GetBytes(const double d,uchar &bytes[])
{
//--- module
double abs_d=MathAbs(d);
//--- number without its fractional
double floor_d=MathFloor(abs_d);
//--- fractional part
double fractional_d=abs_d-floor_d;
//--- variable will store the degree
int power;
//--- exponent shift
double exp_shift;
//--- create variables
int k;
int j;
uchar u;
double step;
double f;
//--- abs_d as bits
bool abs_d_to_bitArray[];
//--- d as bits in format IEEE 754
bool d_to_bitArray[];
//--- allocation
ArrayResizeAL(d_to_bitArray,64);
ArrayResizeAL(bytes,8);
//--- initialization
power=0;
//--- for integer part
while(1)
{
//--- if the number is less than or equal floor_d, we increase the degree
//--- if 2^power > floor_d, then maximal number < floor_d - it 2^(power-1)
if(floor_d>=MathPow(2,power))
power++;
else
break;
}
//--- get power-1
power--;
//--- if power=-1, then floor_d=0
//--- find a negative power for the fractional part
if(power==-1)
{
power=0;
while(1)
{
//--- the same principle as above
if(fractional_d<MathPow(2,power))
power--;
else
break;
}
}
//--- convert decimal to binary
j=0;
step=power;
while(abs_d!=0.0)
{
f=MathPow(2,step);
//--- check
if(f>abs_d)
{
//--- if abs_d < f - this bit is zero
ArrayResize(abs_d_to_bitArray,j+1);
abs_d_to_bitArray[j]=0;
j++;
step-=1;
}
else
{
//--- if abs_d >= f, then this bit is equal one
//--- reduction abs_d
abs_d-=f;
ArrayResize(abs_d_to_bitArray,j+1);
abs_d_to_bitArray[j]=1;
j++;
step-=1;
}
}
//--- according to IEEE 754,
//--- zero bit determines sign of the number, 0 -> '+', 1 -> '-'.
if(d>=0)
d_to_bitArray[0]=0;
else
d_to_bitArray[0]=1;
//--- offset input
exp_shift=1023+power;
//--- bits from the first and 11 are reserved for the shifted exponential
j=1;
for(int i=10; i>=0; i--)
{
if(MathPow(2,i)>exp_shift)
d_to_bitArray[j]=0;
else
{
d_to_bitArray[j]=1;
//--- reduction
exp_shift-=MathPow(2,i);
}
j++;
}
//--- Get the length of the array of the binary representation of abs_d
k=ArraySize(abs_d_to_bitArray);
j=1;
//--- Bits from 12 to 63 are filled with binary representation of abs_b
//--- the first element abs_d_to_bitArray is always 1
for(int i=12; i<64; i++)
{
if(j<k)
{
d_to_bitArray[i]=abs_d_to_bitArray[i-11];
j++;
}
else
d_to_bitArray[i]=0;
}
//--- reverse
ArrayReverse(d_to_bitArray);
//--- converting bit array to byte array
for(int i=0; i<8; i++)
{
u=0;
//--- get byte
for(int t=0; t<8; t++)
u+=(uchar)(d_to_bitArray[i*8+t]*MathPow(2,t));
//--- save byte
bytes[i]=u;
}
}
//+------------------------------------------------------------------+
//| Converting byte array to a integer |
//+------------------------------------------------------------------+
int BitConverter::ToInt32(uchar &bytes[])
{
//--- get size
int size=ArraySize(bytes);
//--- create variables
int d=0;
int mul=256*256*256;
//--- get number
for(int i=size-1; i>=0; i--)
{
d+=bytes[i]*mul;
mul=mul/256;
}
//--- return result
return(d);
}
//+------------------------------------------------------------------+
//| Converting byte array to a double |
//+------------------------------------------------------------------+
double BitConverter::ToDouble(uchar &bytes[])
{
//--- create variables
int s;
//--- exponent shift
int e=0;
//--- mantissa
double m=0;
//--- array of bits in IEEE 754
bool bits[];
ArrayResize(bits,64);
//--- get array of bits from array of bytes
for(int i=0; i<8; i++)
{
for(int j=7; j>=0; j--)
{
//--- if 2 in power >, bits[i*8+j]=0, else bits[i*8+j]=0
if(MathPow(2,j)>bytes[i])
bits[i*8+j]=0;
else
{
bits[i*8+j]=1;
//--- reduction
bytes[i]-=(uchar)MathPow(2,j);
}
}
}
//--- search bits with 1
bool allzero=true;
for(int i=0; i<64; i++)
if(bits[i]==1)
allzero=false;
//--- if all bits are 0, then number is 0
if(allzero==true)
return(0.0);
//--- reverse array
ArrayReverse(bits);
//--- s-the first bit, determines sign of the number
s=bits[0];
//--- calculation exponent shift
for(int i=10; i>=0; i--)
e+=(int)(bits[11-i]*MathPow(2,i));
//--- get mantissa
for(int i=0; i<52; i++)
m+=bits[12+i]*MathPow(2,-1-i);
//--- return result
return(MathPow(-1,s)*MathPow(2,e-1023)*(1+m));
}
//+------------------------------------------------------------------+
//| Byte ordering (forward, backward) |
//+------------------------------------------------------------------+
bool BitConverter::IsLittleEndian(void)
{
//--- forward
return(true);
}
//+------------------------------------------------------------------+
//| Get string from char array |
//+------------------------------------------------------------------+
string GetSelectionString(char &buf[],int startIndex,int lenght)
{
return(CharArrayToString(buf,startIndex,lenght));
}
//+------------------------------------------------------------------+
//| Get sign of number |
//+------------------------------------------------------------------+
double MathSign(const double x)
{
//--- if x>0
if(x>0)
return(1);
//--- if ?==0
if(x==0)
return(0);
//--- ?<0
return(-1);
}
//+------------------------------------------------------------------+
//| Structure stores a variable of type double |
//+------------------------------------------------------------------+
union UDoubleValue
{
double value;
long bits;
UDoubleValue(double dbl): value(dbl) { }
UDoubleValue(long bit_value): bits(bit_value) { }
};
//+------------------------------------------------------------------+
//| Work with infinity and NaN |
//+------------------------------------------------------------------+
class CInfOrNaN
{
public:
//--- checks
static bool IsPositiveInfinity(const double x);
static bool IsNegativeInfinity(const double x);
static bool IsInfinity(const double x);
static bool IsNaN(const double x);
//--- generation values
static double PositiveInfinity(void);
static double NegativeInfinity(void);
static double NaN(void);
};
//+------------------------------------------------------------------+
//| Check for +inf |
//+------------------------------------------------------------------+
bool CInfOrNaN::IsPositiveInfinity(const double x)
{
UDoubleValue val=x;
//--- check
return(val.bits==0x7FF0000000000000);
}
//+------------------------------------------------------------------+
//| Check for -inf |
//+------------------------------------------------------------------+
bool CInfOrNaN::IsNegativeInfinity(const double x)
{
UDoubleValue val=x;
//--- check
return(val.bits==0xFFF0000000000000);
}
//+------------------------------------------------------------------+
//| Check for +-inf |
//+------------------------------------------------------------------+
bool CInfOrNaN::IsInfinity(const double x)
{
return(MathClassify(x)==FP_INFINITE);
}
//+------------------------------------------------------------------+
//| Check for NaN |
//+------------------------------------------------------------------+
bool CInfOrNaN::IsNaN(const double x)
{
return(MathClassify(x)==FP_NAN);
}
//+------------------------------------------------------------------+
//| Return +inf |
//+------------------------------------------------------------------+
double CInfOrNaN::PositiveInfinity(void)
{
UDoubleValue val(0x7FF0000000000000);
return(val.value);
}
//+------------------------------------------------------------------+
//| Return -inf |
//+------------------------------------------------------------------+
double CInfOrNaN::NegativeInfinity(void)
{
UDoubleValue val(0xFFF0000000000000);
return(val.value);
}
//+------------------------------------------------------------------+
//| Return NaN |
//+------------------------------------------------------------------+
double CInfOrNaN::NaN(void)
{
UDoubleValue val(0x7FFFFFFFFFFFFFFF);
return(val.value);
}
//+------------------------------------------------------------------+
Math/Alglib/dataanalysis.mqh
+26387
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/delegatefunctions.mqh
+444
View File
@@ -0,0 +1,444 @@
//+------------------------------------------------------------------+
//| delegatefunctions.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
#include <Object.mqh>
#include "matrix.mqh"
//+------------------------------------------------------------------+
//| Calculates f(arg), stores result to func |
//+------------------------------------------------------------------+
class CNDimensional_Func
{
public:
CNDimensional_Func(void);
~CNDimensional_Func(void);
virtual void Func(double &x[],double &func,CObject &obj);
virtual void Func(CRowDouble &x,double &func,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_Func::CNDimensional_Func(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_Func::~CNDimensional_Func(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Func::Func(double &x[],double &func,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
void CNDimensional_Func::Func(CRowDouble &x,double &func,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| calculates func = f(arg), grad[i] = df(arg)/d(arg[i]) |
//+------------------------------------------------------------------+
class CNDimensional_Grad
{
public:
//--- constructor, destructor
CNDimensional_Grad(void);
~CNDimensional_Grad(void);
//--- virtual method
virtual void Grad(double &x[],double &func,double &grad[],CObject &obj);
virtual void Grad(CRowDouble &x,double &func,CRowDouble &grad,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_Grad::CNDimensional_Grad(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_Grad::~CNDimensional_Grad(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Grad::Grad(double &x[],double &func,double &grad[],
CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Grad::Grad(CRowDouble &x,double &func,CRowDouble &grad,
CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates func = f(arg), grad[i] = df(arg)/d(arg[i]), |
//| hess[i,j] = d2f(arg)/(d(arg[i])*d(arg[j])) |
//+------------------------------------------------------------------+
class CNDimensional_Hess
{
public:
//--- constructor, destructor
CNDimensional_Hess(void);
~CNDimensional_Hess(void);
//--- virtual method
virtual void Hess(double &x[],double &func,double &grad[],CMatrixDouble &hess,CObject &obj);
virtual void Hess(CRowDouble &x,double &func,CRowDouble &grad,CMatrixDouble &hess,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_Hess::CNDimensional_Hess(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_Hess::~CNDimensional_Hess(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Hess::Hess(double &x[],double &func,double &grad[],
CMatrixDouble &hess,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Hess::Hess(CRowDouble &x,double &func,CRowDouble &grad,
CMatrixDouble &hess,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates vector function f(arg), stores result to fi |
//+------------------------------------------------------------------+
class CNDimensional_FVec
{
public:
//--- constructor, destructor
CNDimensional_FVec(void);
~CNDimensional_FVec(void);
//--- virtual method
virtual void FVec(double &x[],double &fi[],CObject &obj);
virtual void FVec(CRowDouble &x,CRowDouble &fi,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_FVec::CNDimensional_FVec(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_FVec::~CNDimensional_FVec(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_FVec::FVec(double &x[],double &fi[],CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_FVec::FVec(CRowDouble &x,CRowDouble &fi,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates f[i] = fi(arg), jac[i,j] = df[i](arg)/d(arg[j]) |
//+------------------------------------------------------------------+
class CNDimensional_Jac
{
public:
//--- constructor, destructor
CNDimensional_Jac(void);
~CNDimensional_Jac(void);
//--- virtual method
virtual void Jac(double &x[],double &fi[],CMatrixDouble &jac,CObject &obj);
virtual void Jac(CRowDouble &x,CRowDouble &fi,CMatrixDouble &jac,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_Jac::CNDimensional_Jac(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_Jac::~CNDimensional_Jac(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Jac::Jac(double &x[],double &fi[],CMatrixDouble &jac,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Jac::Jac(CRowDouble &x,CRowDouble &fi,CMatrixDouble &jac,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates f(p,q), stores result to func |
//+------------------------------------------------------------------+
class CNDimensional_PFunc
{
public:
//--- constructor, destructor
CNDimensional_PFunc(void);
~CNDimensional_PFunc(void);
//--- virtual method
virtual void PFunc(double &c[],double &x[],double &func,CObject &obj);
virtual void PFunc(CRowDouble &c,CRowDouble &x,double &func,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_PFunc::CNDimensional_PFunc(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_PFunc::~CNDimensional_PFunc(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PFunc::PFunc(double &c[],double &x[],double &func,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PFunc::PFunc(CRowDouble &c,CRowDouble &x,double &func,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates func = f(p,q), grad[i] = df(p,q)/d(p[i]) |
//+------------------------------------------------------------------+
class CNDimensional_PGrad
{
public:
//--- constructor, destructor
CNDimensional_PGrad(void);
~CNDimensional_PGrad(void);
//--- virtual method
virtual void PGrad(double &c[],double &x[],double &func,double &grad[],CObject &obj);
virtual void PGrad(CRowDouble &c,CRowDouble &x,double &func,CRowDouble &grad,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_PGrad::CNDimensional_PGrad(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_PGrad::~CNDimensional_PGrad(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PGrad::PGrad(double &c[],double &x[],double &func,
double &grad[],CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PGrad::PGrad(CRowDouble &c,CRowDouble &x,double &func,
CRowDouble &grad,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Calculates func = f(p,q), grad[i] = df(p,q)/d(p[i]), |
//| hess[i,j] = d2f(p,q)/(d(p[i])*d(p[j])) |
//+------------------------------------------------------------------+
class CNDimensional_PHess
{
public:
//--- constructor, destructor
CNDimensional_PHess(void);
~CNDimensional_PHess(void);
//--- virtual method
virtual void PHess(double &c[],double &x[],double &func,double &grad[],CMatrixDouble &hess,CObject &obj);
virtual void PHess(CRowDouble &c,CRowDouble &x,double &func,CRowDouble &grad,CMatrixDouble &hess,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_PHess::CNDimensional_PHess(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_PHess::~CNDimensional_PHess(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PHess::PHess(double &c[],double &x[],double &func,
double &grad[],CMatrixDouble &hess,
CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_PHess::PHess(CRowDouble &c,CRowDouble &x,double &func,
CRowDouble &grad,CMatrixDouble &hess,
CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Callbacks for ODE solvers: calculates dy/dx for given y[] and x |
//+------------------------------------------------------------------+
class CNDimensional_ODE_RP
{
public:
//--- constructor, destructor
CNDimensional_ODE_RP(void);
~CNDimensional_ODE_RP(void);
//--- virtual method
virtual void ODE_RP(double &y[],double x,double &dy[],CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_ODE_RP::CNDimensional_ODE_RP(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_ODE_RP::~CNDimensional_ODE_RP(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_ODE_RP::ODE_RP(double &y[],double x,double &dy[],
CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Callbacks for integrators: calculates f(x) for given x |
//| (additional parameters xminusa and bminusx contain x-a and b-x) |
//+------------------------------------------------------------------+
class CIntegrator1_Func
{
public:
//--- constructor, destructor
CIntegrator1_Func(void);
~CIntegrator1_Func(void);
//--- virtual method
virtual void Int_Func(double x,double xminusa,double bminusx,double &y,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CIntegrator1_Func::CIntegrator1_Func(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CIntegrator1_Func::~CIntegrator1_Func(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CIntegrator1_Func::Int_Func(double x,double xminusa,double bminusx,
double &y,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| Callbacks for progress reports: reports current position of |
//| optimization algo |
//+------------------------------------------------------------------+
class CNDimensional_Rep
{
public:
//--- constructor, destructor
CNDimensional_Rep(void);
~CNDimensional_Rep(void);
//--- virtual method
virtual void Rep(double &arg[],double func,CObject &obj);
virtual void Rep(CRowDouble &arg,double func,CObject &obj);
};
//+------------------------------------------------------------------+
//| Constructor without parameters |
//+------------------------------------------------------------------+
CNDimensional_Rep::CNDimensional_Rep(void)
{
}
//+------------------------------------------------------------------+
//| Destructor |
//+------------------------------------------------------------------+
CNDimensional_Rep::~CNDimensional_Rep(void)
{
}
//+------------------------------------------------------------------+
//| Empty function body |
//+------------------------------------------------------------------+
void CNDimensional_Rep::Rep(double &arg[],double func,CObject &obj)
{
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
void CNDimensional_Rep::Rep(CRowDouble &arg,double func,CObject &obj)
{
}
//+------------------------------------------------------------------+
Math/Alglib/diffequations.mqh
+849
View File
@@ -0,0 +1,849 @@
//+------------------------------------------------------------------+
//| diffequations.mqh |
//| Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//| Copyright 2012-2023, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5 |
//| |
//| The features of the library include: |
//| - Linear algebra (direct algorithms, EVD, SVD) |
//| - Solving systems of linear and non-linear equations |
//| - Interpolation |
//| - Optimization |
//| - FFT (Fast Fourier Transform) |
//| - Numerical integration |
//| - Linear and nonlinear least-squares fitting |
//| - Ordinary differential equations |
//| - Computation of special functions |
//| - Descriptive statistics and hypothesis testing |
//| - Data analysis - classification, regression |
//| - Implementing linear algebra algorithms, interpolation, etc. |
//| in high-precision arithmetic (using MPFR) |
//| |
//| This file is free software; you can redistribute it and/or |
//| modify it under the terms of the GNU General Public License as |
//| published by the Free Software Foundation (www.fsf.org); either |
//| version 2 of the License, or (at your option) any later version. |
//| |
//| This program is distributed in the hope that it will be useful, |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
//| GNU General Public License for more details. |
//+------------------------------------------------------------------+
#include "matrix.mqh"
#include "alglibinternal.mqh"
//+------------------------------------------------------------------+
//| Auxiliary class for CODESolver |
//+------------------------------------------------------------------+
class CODESolverState
{
public:
int m_n;
int m_m;
double m_xscale;
double m_h;
double m_eps;
bool m_fraceps;
int m_repterminationtype;
int m_repnfev;
int m_solvertype;
bool m_needdy;
double m_x;
RCommState m_rstate;
//--- arrays
double m_yc[];
double m_escale[];
double m_xg[];
double m_y[];
double m_dy[];
double m_yn[];
double m_yns[];
double m_rka[];
double m_rkc[];
double m_rkcs[];
//--- matrices
CMatrixDouble m_ytbl;
CMatrixDouble m_rkb;
CMatrixDouble m_rkk;
public:
CODESolverState(void);
~CODESolverState(void) {}
void Copy(CODESolverState &obj);
};
//+------------------------------------------------------------------+
//| Constructor |
//+------------------------------------------------------------------+
CODESolverState::CODESolverState(void)
{
m_n=0;
m_m=0;
m_xscale=0;
m_h=0;
m_eps=0;
m_fraceps=false;
m_repterminationtype=0;
m_repnfev=0;
m_solvertype=0;
m_needdy=false;
m_x=0;
}
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CODESolverState::Copy(CODESolverState &obj)
{
//--- copy variables
m_n=obj.m_n;
m_m=obj.m_m;
m_xscale=obj.m_xscale;
m_h=obj.m_h;
m_eps=obj.m_eps;
m_fraceps=obj.m_fraceps;
m_repterminationtype=obj.m_repterminationtype;
m_repnfev=obj.m_repnfev;
m_solvertype=obj.m_solvertype;
m_needdy=obj.m_needdy;
m_x=obj.m_x;
m_rstate.Copy(obj.m_rstate);
//--- copy arrays
ArrayCopy(m_yc,obj.m_yc);
ArrayCopy(m_escale,obj.m_escale);
ArrayCopy(m_xg,obj.m_xg);
ArrayCopy(m_y,obj.m_y);
ArrayCopy(m_dy,obj.m_dy);
ArrayCopy(m_yn,obj.m_yn);
ArrayCopy(m_yns,obj.m_yns);
ArrayCopy(m_rka,obj.m_rka);
ArrayCopy(m_rkc,obj.m_rkc);
ArrayCopy(m_rkcs,obj.m_rkcs);
//--- copy matrices
m_ytbl=obj.m_ytbl;
m_rkb=obj.m_rkb;
m_rkk=obj.m_rkk;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CODESolverState |
//+------------------------------------------------------------------+
class CODESolverStateShell
{
private:
CODESolverState m_innerobj;
public:
//--- constructors, destructor
CODESolverStateShell(void) {}
CODESolverStateShell(CODESolverState &obj) { m_innerobj.Copy(obj); }
~CODESolverStateShell(void) {}
//--- methods
bool GetNeedDY(void);
void SetNeedDY(const bool b);
double GetX(void);
void SetX(const double d);
CODESolverState *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable needdy |
//+------------------------------------------------------------------+
bool CODESolverStateShell::GetNeedDY(void)
{
return(m_innerobj.m_needdy);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable needdy |
//+------------------------------------------------------------------+
void CODESolverStateShell::SetNeedDY(const bool b)
{
m_innerobj.m_needdy=b;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable x |
//+------------------------------------------------------------------+
double CODESolverStateShell::GetX(void)
{
return(m_innerobj.m_x);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable x |
//+------------------------------------------------------------------+
void CODESolverStateShell::SetX(const double d)
{
m_innerobj.m_x=d;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CODESolverState *CODESolverStateShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Auxiliary class for CODESolver |
//+------------------------------------------------------------------+
class CODESolverReport
{
public:
//--- class variables
int m_nfev;
int m_terminationtype;
//--- constructor, destructor
CODESolverReport(void) { ZeroMemory(this); }
~CODESolverReport(void) {}
//--- copy
void Copy(CODESolverReport &obj);
};
//+------------------------------------------------------------------+
//| Copy |
//+------------------------------------------------------------------+
void CODESolverReport::Copy(CODESolverReport &obj)
{
//--- copy variables
m_nfev=obj.m_nfev;
m_terminationtype=obj.m_terminationtype;
}
//+------------------------------------------------------------------+
//| This class is a shell for class CODESolverReport |
//+------------------------------------------------------------------+
class CODESolverReportShell
{
private:
CODESolverReport m_innerobj;
public:
//--- constructor, destructor
CODESolverReportShell(void) {}
CODESolverReportShell(CODESolverReport &obj) { m_innerobj.Copy(obj); }
~CODESolverReportShell(void) {}
//--- methods
int GetNFev(void);
void SetNFev(const int i);
int GetTerminationType(void);
void SetTerminationType(const int i);
CODESolverReport *GetInnerObj(void);
};
//+------------------------------------------------------------------+
//| Returns the value of the variable nfev |
//+------------------------------------------------------------------+
int CODESolverReportShell::GetNFev(void)
{
return(m_innerobj.m_nfev);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable nfev |
//+------------------------------------------------------------------+
void CODESolverReportShell::SetNFev(const int i)
{
m_innerobj.m_nfev=i;
}
//+------------------------------------------------------------------+
//| Returns the value of the variable terminationtype |
//+------------------------------------------------------------------+
int CODESolverReportShell::GetTerminationType(void)
{
return(m_innerobj.m_terminationtype);
}
//+------------------------------------------------------------------+
//| Changing the value of the variable terminationtype |
//+------------------------------------------------------------------+
void CODESolverReportShell::SetTerminationType(const int i)
{
m_innerobj.m_terminationtype=i;
}
//+------------------------------------------------------------------+
//| Return object of class |
//+------------------------------------------------------------------+
CODESolverReport *CODESolverReportShell::GetInnerObj(void)
{
return(GetPointer(m_innerobj));
}
//+------------------------------------------------------------------+
//| Solution of ordinary differential equations |
//+------------------------------------------------------------------+
class CODESolver
{
public:
//--- class constants
static const double m_odesolvermaxgrow;
static const double m_odesolvermaxshrink;
//--- public methods
static void ODESolverRKCK(double &y[],const int n,double &x[],const int m,const double eps,const double h,CODESolverState &state);
static void ODESolverResults(CODESolverState &state,int &m,double &xtbl[],CMatrixDouble &ytbl,CODESolverReport &rep);
static bool ODESolverIteration(CODESolverState &state);
private:
//--- private method
static void ODESolverInit(int solvertype,double &y[],const int n,double &x[],const int m,const double eps,double h,CODESolverState &state);
//--- auxiliary functions for ODESolverIteration
static void Func_lbl_rcomm(CODESolverState &state,int n,int m,int i,int j,int k,int klimit,bool gridpoint,double xc,double v,double h,double h2,double err,double maxgrowpow);
static bool Func_lbl_6(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
static bool Func_lbl_8(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
static bool Func_lbl_10(CODESolverState &state,int &n,int &m,int &i,int &j,int &k,int &klimit,bool &gridpoint,double &xc,double &v,double &h,double &h2,double &err,double &maxgrowpow);
};
//+------------------------------------------------------------------+
//| Initialize constants |
//+------------------------------------------------------------------+
const double CODESolver::m_odesolvermaxgrow=3.0;
const double CODESolver::m_odesolvermaxshrink=10.0;
//+------------------------------------------------------------------+
//| Cash-Karp adaptive ODE solver. |
//| This subroutine solves ODE Y'=f(Y,x) with initial conditions |
//| Y(xs)=Ys (here Y may be single variable or vector of N variables)|
//| INPUT PARAMETERS: |
//| Y - initial conditions, array[0..N-1]. |
//| contains values of Y[] at X[0] |
//| N - system size |
//| X - points at which Y should be tabulated, |
//| array[0..M-1] integrations starts at X[0], ends |
//| at X[M-1], intermediate values at X[i] are |
//| returned too. |
//| SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!|
//| M - number of intermediate points + first point + |
//| last point: |
//| * M>2 means that you need both Y(X[M-1]) and M-2 |
//| values at intermediate points |
//| * M=2 means that you want just to integrate from |
//| X[0] to X[1] and don't interested in |
//| intermediate values. |
//| * M=1 means that you don't want to integrate :) |
//| it is degenerate case, but it will be handled |
//| correctly. |
//| * M<1 means error |
//| Eps - tolerance (absolute/relative error on each step |
//| will be less than Eps). When passing: |
//| * Eps>0, it means desired ABSOLUTE error |
//| * Eps<0, it means desired RELATIVE error. |
//| Relative errors are calculated with respect to |
//| maximum values of Y seen so far. Be careful to |
//| use this criterion when starting from Y[] that |
//| are close to zero. |
//| H - initial step lenth, it will be adjusted |
//| automatically after the first step. If H=0, step |
//| will be selected automatically (usualy it will |
//| be equal to 0.001 of min(x[i]-x[j])). |
//| OUTPUT PARAMETERS |
//| State - structure which stores algorithm state between |
//| subsequent calls of OdeSolverIteration. Used |
//| for reverse communication. This structure should |
//| be passed to the OdeSolverIteration subroutine. |
//| SEE ALSO |
//| AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults|
//+------------------------------------------------------------------+
void CODESolver::ODESolverRKCK(double &y[],const int n,double &x[],
const int m,const double eps,
const double h,CODESolverState &state)
{
//--- check
if(!CAp::Assert(n>=1,"ODESolverRKCK: N<1!"))
return;
//--- check
if(!CAp::Assert(m>=1,"ODESolverRKCK: M<1!"))
return;
//--- check
if(!CAp::Assert(CAp::Len(y)>=n,"ODESolverRKCK: Length(Y)<N!"))
return;
//--- check
if(!CAp::Assert(CAp::Len(x)>=m,"ODESolverRKCK: Length(X)<M!"))
return;
//--- check
if(!CAp::Assert(CApServ::IsFiniteVector(y,n),"ODESolverRKCK: Y contains infinite or NaN values!"))
return;
//--- check
if(!CAp::Assert(CApServ::IsFiniteVector(x,m),"ODESolverRKCK: Y contains infinite or NaN values!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(eps),"ODESolverRKCK: Eps is not finite!"))
return;
//--- check
if(!CAp::Assert(eps!=0.0,"ODESolverRKCK: Eps is zero!"))
return;
//--- check
if(!CAp::Assert(CMath::IsFinite(h),"ODESolverRKCK: H is not finite!"))
return;
//--- function call
ODESolverInit(0,y,n,x,m,eps,h,state);
}
//+------------------------------------------------------------------+
//| ODE solver results |
//| Called after OdeSolverIteration returned False. |
//| INPUT PARAMETERS: |
//| State - algorithm state (used by OdeSolverIteration). |
//| OUTPUT PARAMETERS: |
//| M - number of tabulated values, M>=1 |
//| XTbl - array[0..M-1], values of X |
//| YTbl - array[0..M-1,0..N-1], values of Y in X[i] |
//| Rep - solver report: |
//| * Rep.TerminationType completetion code: |
//| * -2 X is not ordered by |
//| ascending/descending or there are |
//| non-distinct X[], i.e. X[i]=X[i+1] |
//| * -1 incorrect parameters were specified |
//| * 1 task has been solved |
//| * Rep.NFEV contains number of function |
//| calculations |
//+------------------------------------------------------------------+
void CODESolver::ODESolverResults(CODESolverState &state,int &m,
double &xtbl[],CMatrixDouble &ytbl,
CODESolverReport &rep)
{
//--- create variables
double v=0;
int i=0;
int i_=0;
//--- initialization
m=0;
rep.m_terminationtype=state.m_repterminationtype;
//--- check
if(rep.m_terminationtype>0)
{
//--- change values
m=state.m_m;
rep.m_nfev=state.m_repnfev;
//--- allocation
ArrayResize(xtbl,state.m_m);
v=state.m_xscale;
//--- calculation
for(i_=0; i_<=state.m_m-1; i_++)
xtbl[i_]=v*state.m_xg[i_];
//--- allocation
ytbl.Resize(state.m_m,state.m_n);
for(i=0; i<=state.m_m-1; i++)
{
for(i_=0; i_<=state.m_n-1; i_++)
ytbl.Set(i,i_,state.m_ytbl[i][i_]);
}
}
else
rep.m_nfev=0;
}
//+------------------------------------------------------------------+
//| Internal initialization subroutine |
//+------------------------------------------------------------------+
void CODESolver::ODESolverInit(int solvertype,double &y[],const int n,
double &x[],const int m,const double eps,
double h,CODESolverState &state)
{
//--- create variables
int i=0;
double v=0;
int i_=0;
//--- Prepare RComm
state.m_rstate.ia.Resize(6);
ArrayResizeAL(state.m_rstate.ba,1);
state.m_rstate.ra.Resize(6);
state.m_rstate.stage=-1;
state.m_needdy=false;
//--- check parameters.
if((n<=0 || m<1) || eps==0.0)
{
state.m_repterminationtype=-1;
return;
}
//--- check
if(h<0.0)
h=-h;
//--- quick exit if necessary.
//--- after this block we assume that M>1
if(m==1)
{
//--- change values
state.m_repnfev=0;
state.m_repterminationtype=1;
state.m_ytbl.Resize(1,n);
for(i_=0; i_<=n-1; i_++)
state.m_ytbl.Set(0,i_,y[i_]);
//--- allocation
ArrayResize(state.m_xg,m);
for(i_=0; i_<=m-1; i_++)
state.m_xg[i_]=x[i_];
//--- exit the function
return;
}
//--- check again: correct order of X[]
if(x[1]==x[0])
{
state.m_repterminationtype=-2;
return;
}
for(i=1; i<=m-1; i++)
{
//--- check
if((x[1]>x[0] && x[i]<=x[i-1]) || (x[1]<x[0] && x[i]>=x[i-1]))
{
state.m_repterminationtype=-2;
return;
}
}
//--- auto-select H if necessary
if(h==0.0)
{
v=MathAbs(x[1]-x[0]);
for(i=2; i<=m-1; i++)
v=MathMin(v,MathAbs(x[i]-x[i-1]));
h=0.001*v;
}
//--- store parameters
state.m_n=n;
state.m_m=m;
state.m_h=h;
state.m_eps=MathAbs(eps);
state.m_fraceps=eps<0.0;
//--- allocation
ArrayResize(state.m_xg,m);
for(i_=0; i_<=m-1; i_++)
state.m_xg[i_]=x[i_];
//--- check
if(x[1]>x[0])
state.m_xscale=1;
else
{
state.m_xscale=-1;
for(i_=0; i_<=m-1; i_++)
state.m_xg[i_]=-1*state.m_xg[i_];
}
//--- allocation
ArrayResize(state.m_yc,n);
for(i_=0; i_<=n-1; i_++)
state.m_yc[i_]=y[i_];
//--- change values
state.m_solvertype=solvertype;
state.m_repterminationtype=0;
//--- Allocate arrays
ArrayResize(state.m_y,n);
ArrayResize(state.m_dy,n);
}
//+------------------------------------------------------------------+
//| Iterative method |
//+------------------------------------------------------------------+
bool CODESolver::ODESolverIteration(CODESolverState &state)
{
//--- create variables
int n=0;
int m=0;
int i=0;
int j=0;
int k=0;
double xc=0;
double v=0;
double h=0;
double h2=0;
bool gridpoint;
double err=0;
double maxgrowpow=0;
int klimit=0;
int i_=0;
//--- This code initializes locals by:
//--- * random values determined during code
//--- generation - on first subroutine call
//--- * values from previous call - on subsequent calls
if(state.m_rstate.stage>=0)
{
//--- initialization
n=state.m_rstate.ia[0];
m=state.m_rstate.ia[1];
i=state.m_rstate.ia[2];
j=state.m_rstate.ia[3];
k=state.m_rstate.ia[4];
klimit=state.m_rstate.ia[5];
gridpoint=state.m_rstate.ba[0];
xc=state.m_rstate.ra[0];
v=state.m_rstate.ra[1];
h=state.m_rstate.ra[2];
h2=state.m_rstate.ra[3];
err=state.m_rstate.ra[4];
maxgrowpow=state.m_rstate.ra[5];
}
else
{
//--- initialization
n=-983;
m=-989;
i=-834;
j=900;
k=-287;
klimit=364;
gridpoint=false;
xc=-338;
v=-686;
h=912;
h2=585;
err=497;
maxgrowpow=-271;
}
//--- check
if(state.m_rstate.stage==0)
{
//--- change values
state.m_needdy=false;
state.m_repnfev=state.m_repnfev+1;
v=h*state.m_xscale;
for(i_=0; i_<=n-1; i_++)
state.m_rkk.Set(k,i_,v*state.m_dy[i_]);
//--- update YN/YNS
v=state.m_rkc[k];
for(i_=0; i_<=n-1; i_++)
state.m_yn[i_]=state.m_yn[i_]+v*state.m_rkk[k][i_];
v=state.m_rkcs[k];
for(i_=0; i_<=n-1; i_++)
state.m_yns[i_]=state.m_yns[i_]+v*state.m_rkk[k][i_];
k=k+1;
return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//--- Routine body
//--- prepare
if(state.m_repterminationtype!=0)
return(false);
//--- change values
n=state.m_n;
m=state.m_m;
h=state.m_h;
maxgrowpow=MathPow(m_odesolvermaxgrow,5);
state.m_repnfev=0;
//--- some preliminary checks for internal errors
//--- after this we assume that H>0 and M>1
if(!CAp::Assert(state.m_h>0.0,"ODESolver: internal error"))
return(false);
//--- check
if(!CAp::Assert(m>1,"ODESolverIteration: internal error"))
return(false);
//--- choose solver
if(state.m_solvertype!=0)
return(false);
//--- Cask-Karp solver
//--- Prepare coefficients table.
//--- Check it for errors
ArrayResize(state.m_rka,6);
//--- calculation
state.m_rka[0]=0;
state.m_rka[1]=1.0/5.0;
state.m_rka[2]=3.0/10.0;
state.m_rka[3]=3.0/5.0;
state.m_rka[4]=1;
state.m_rka[5]=7.0/8.0;
state.m_rkb.Resize(6,5);
state.m_rkb.Set(1,0,1.0/5.0);
state.m_rkb.Set(2,0,3.0/40.0);
state.m_rkb.Set(2,1,9.0/40.0);
state.m_rkb.Set(3,0,3.0/10.0);
state.m_rkb.Set(3,1,-(9.0/10.0));
state.m_rkb.Set(3,2,6.0/5.0);
state.m_rkb.Set(4,0,-(11.0/54.0));
state.m_rkb.Set(4,1,5.0/2.0);
state.m_rkb.Set(4,2,-(70.0/27.0));
state.m_rkb.Set(4,3,35.0/27.0);
state.m_rkb.Set(5,0,1631.0/55296.0);
state.m_rkb.Set(5,1,175.0/512.0);
state.m_rkb.Set(5,2,575.0/13824.0);
state.m_rkb.Set(5,3,44275.0/110592.0);
state.m_rkb.Set(5,4,253.0/4096.0);
//--- allocation
ArrayResize(state.m_rkc,6);
//--- calculation
state.m_rkc[0]=37.0/378.0;
state.m_rkc[1]=0;
state.m_rkc[2]=250.0/621.0;
state.m_rkc[3]=125.0/594.0;
state.m_rkc[4]=0;
state.m_rkc[5]=512.0/1771.0;
//--- allocation
ArrayResize(state.m_rkcs,6);
//--- calculation
state.m_rkcs[0]=2825.0/27648.0;
state.m_rkcs[1]=0;
state.m_rkcs[2]=18575.0/48384.0;
state.m_rkcs[3]=13525.0/55296.0;
state.m_rkcs[4]=277.0/14336.0;
state.m_rkcs[5]=1.0/4.0;
state.m_rkk.Resize(6,n);
//--- Main cycle consists of two iterations:
//--- * outer where we travel from X[i-1] to X[i]
//--- * inner where we travel inside [X[i-1],X[i]]
state.m_ytbl.Resize(m,n);
ArrayResize(state.m_escale,n);
ArrayResize(state.m_yn,n);
ArrayResize(state.m_yns,n);
//--- change value
xc=state.m_xg[0];
for(i_=0; i_<=n-1; i_++)
state.m_ytbl.Set(0,i_,state.m_yc[i_]);
for(j=0; j<=n-1; j++)
state.m_escale[j]=0;
i=1;
//--- check
if(i>m-1)
{
state.m_repterminationtype=1;
//--- return result
return(false);
}
//--- begin inner iteration
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get |
//| rid of the operator unconditional jump goto |
//+------------------------------------------------------------------+
void CODESolver::Func_lbl_rcomm(CODESolverState &state,int n,int m,
int i,int j,int k,int klimit,
bool gridpoint,double xc,double v,
double h,double h2,double err,
double maxgrowpow)
{
//--- save
state.m_rstate.ia.Set(0,n);
state.m_rstate.ia.Set(1,m);
state.m_rstate.ia.Set(2,i);
state.m_rstate.ia.Set(3,j);
state.m_rstate.ia.Set(4,k);
state.m_rstate.ia.Set(5,klimit);
state.m_rstate.ba[0]= gridpoint;
state.m_rstate.ra.Set(0,xc);
state.m_rstate.ra.Set(1,v);
state.m_rstate.ra.Set(2,h);
state.m_rstate.ra.Set(3,h2);
state.m_rstate.ra.Set(4,err);
state.m_rstate.ra.Set(5,maxgrowpow);
}
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get |
//| rid of the operator unconditional jump goto |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_6(CODESolverState &state,int &n,int &m,
int &i,int &j,int &k,int &klimit,
bool &gridpoint,double &xc,double &v,
double &h,double &h2,double &err,
double &maxgrowpow)
{
//--- truncate step if needed (beyond right boundary).
//--- determine should we store X or not
if(xc+h>=state.m_xg[i])
{
h=state.m_xg[i]-xc;
gridpoint=true;
}
else
gridpoint=false;
//--- Update error scale maximums
//--- These maximums are initialized by zeros,
//--- then updated every iterations.
for(j=0; j<=n-1; j++)
state.m_escale[j]=MathMax(state.m_escale[j],MathAbs(state.m_yc[j]));
//--- make one step:
//--- 1. calculate all info needed to do step
//--- 2. update errors scale maximums using values/derivatives
//--- obtained during (1)
//--- Take into account that we use scaling of X to reduce task
//--- to the form where x[0] < x[1] < ... < x[n-1]. So X is
//--- replaced by x=xscale*t,and dy/dx=f(y,x) is replaced
//--- by dy/dt=xscale*f(y,xscale*t).
for(int i_=0; i_<=n-1; i_++)
state.m_yn[i_]=state.m_yc[i_];
for(int i_=0; i_<=n-1; i_++)
state.m_yns[i_]=state.m_yc[i_];
k=0;
//--- function call, return result
return(Func_lbl_8(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get |
//| rid of the operator unconditional jump goto |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_8(CODESolverState &state,int &n,int &m,
int &i,int &j,int &k,int &klimit,
bool &gridpoint,double &xc,double &v,
double &h,double &h2,double &err,
double &maxgrowpow)
{
//--- check
if(k>5)
return(Func_lbl_10(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
//--- prepare data for the next update of YN/YNS
state.m_x=state.m_xscale*(xc+state.m_rka[k]*h);
for(int i_=0; i_<=n-1; i_++)
state.m_y[i_]=state.m_yc[i_];
//--- calculation
for(j=0; j<=k-1; j++)
{
v=state.m_rkb[k][j];
for(int i_=0; i_<=n-1; i_++)
state.m_y[i_]=state.m_y[i_]+v*state.m_rkk[j][i_];
}
state.m_needdy=true;
state.m_rstate.stage=0;
//--- Saving state
Func_lbl_rcomm(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow);
//--- return result
return(true);
}
//+------------------------------------------------------------------+
//| Auxiliary function for ODESolverIteration. Is a product to get |
//| rid of the operator unconditional jump goto |
//+------------------------------------------------------------------+
bool CODESolver::Func_lbl_10(CODESolverState &state,int &n,int &m,
int &i,int &j,int &k,int &klimit,
bool &gridpoint,double &xc,double &v,
double &h,double &h2,double &err,
double &maxgrowpow)
{
//--- estimate error
err=0;
for(j=0; j<=n-1; j++)
{
//--- check
if(!state.m_fraceps)
{
//--- absolute error is estimated
err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j]));
}
else
{
//--- Relative error is estimated
v=state.m_escale[j];
//--- check
if(v==0.0)
v=1;
err=MathMax(err,MathAbs(state.m_yn[j]-state.m_yns[j])/v);
}
}
//--- calculate new step,restart if necessary
if(maxgrowpow*err<=state.m_eps)
h2=m_odesolvermaxgrow*h;
else
h2=h*MathPow(state.m_eps/err,0.2);
//--- check
if(h2<h/m_odesolvermaxshrink)
h2=h/m_odesolvermaxshrink;
//--- check
if(err>state.m_eps)
{
h=h2;
//--- begin inner iteration
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//--- advance position
xc=xc+h;
for(int i_=0; i_<=n-1; i_++)
state.m_yc[i_]=state.m_yn[i_];
//--- update H
h=h2;
//--- break on grid point
if(gridpoint)
{
//--- save result
for(int i_=0; i_<=n-1; i_++)
state.m_ytbl.Set(i,i_,state.m_yc[i_]);
i=i+1;
//--- check
if(i>m-1)
{
state.m_repterminationtype=1;
return(false);
}
//--- begin inner iteration
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//--- begin inner iteration
return(Func_lbl_6(state,n,m,i,j,k,klimit,gridpoint,xc,v,h,h2,err,maxgrowpow));
}
//+------------------------------------------------------------------+
Math/Alglib/fasttransforms.mqh
+2657
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/integration.mqh
+2980
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/interpolation.mqh
+34261
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/linalg.mqh
+37629
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/matrix.mqh
+1025
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/optimization.mqh
+51766
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/solvers.mqh
+7154
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/specialfunctions.mqh
+6534
View File
File diff suppressed because it is too large Load Diff
Math/Alglib/statistics.mqh
+13650
View File
File diff suppressed because it is too large Load Diff