ufranco/Article-20589-Volatility-Modeling-Part-1
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
Article-20589-Volatility-Mo.../Include/np.mqh

4651 lines
290 KiB
MQL5
Raw Permalink Blame History

//+------------------------------------------------------------------+
//| np.mqh |
//| Copyright 2024, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link "https://www.mql5.com"
#include<Graphics\Graphic.mqh>
#include<Math\Stat\Normal.mqh>
#include<Math\Alglib\alglib.mqh>
#define _NP_
#define BEGIN 0
#define STEP 1
#define END LONG_MAX
#define BEGIN_REVERSE -1
#define STEP_REVERSE -1
#define END_REVERSE LONG_MIN
//+------------------------------------------------------------------+
//| library constants |
//+------------------------------------------------------------------+
namespace np
{
namespace internal
{
//+------------------------------------------------------------------+
//|find index of a sorted array such that arr[i] <= key < arr[i + 1] |
//+------------------------------------------------------------------+
int binary_search(const double key, const double &arr[])
{
int imin = 0;
int imax = ArraySize(arr)-1;
int mid = -1;
/* Handle keys outside of the arr range first */
if(key > arr[arr.Size()-1])
{
return ArraySize(arr);
}
else
if(key < arr[0])
{
return -1;
}
while(imin<=imax)
{
mid = (imin+imax)/2;
if(mid<int(arr.Size()-1) && arr[mid]<=key && arr[mid+1]>key)
return mid;
else
{
if(key<arr[mid])
imax = mid - 1;
else
imin = mid + 1;
}
}
return ArraySize(arr) - 2;
}
//+------------------------------------------------------------------+
//| cast string data to vector format doubles and dates |
//+------------------------------------------------------------------+
vector formatCol(string &Arr[],bool handle_dates,int date_column, int column_index)
{
int size = ArraySize(Arr);
vector ret(size);
//---
if(handle_dates && date_column == column_index)
{
for(int i=0; i<size; ++i)
ret[i] = (double)StringToTime(Arr[i]);
}
else
{
for(int i=0; i<size; ++i)
ret[i] = StringToDouble(Arr[i]);
}
//---
return ret;
}
//+------------------------------------------------------------------+
//|get column |
//+------------------------------------------------------------------+
template<typename T>
void col(const T &Matrix[], T &Col[], int column, int cols)
{
int rows = ArraySize(Matrix)/cols;
ArrayResize(Col,rows);
int start = 0;
for(int i=0; i<cols; ++i)
{
start = i;
if(i != column-1)
continue;
else
for(int j=0; j<rows; ++j)
{
//printf("ColMatrix[%d} Matrix{%d]",j,start);
Col[j] = Matrix[start];
start += cols;
}
}
}
//+------------------------------------------------------------------+
//| product for arrays |
//+------------------------------------------------------------------+
template<typename T>
T prodArray(T &in_array[], uint from=0)
{
T result = T(1);
for(uint i = from; i<in_array.Size(); ++i)
result*=in_array[i];
return result;
}
//+------------------------------------------------------------------+
//|calculate size of container |
//+------------------------------------------------------------------+
long calculatesize(long initial,long stop,long step)
{
//---calculate len
long len = fabs(initial-stop)/fabs(step);
//---
if(len*fabs(step) < fabs(initial-stop))
len += (fabs(initial-stop) - len*fabs(step))/fabs(step);
//---
return len;
}
//+------------------------------------------------------------------+
//| process start parameter |
//+------------------------------------------------------------------+
long normalizestart(long start,ulong size)
{
long s = (start<0 && ulong(fabs(start))<size)?long(size)+start:(start<0 && ulong(fabs(start))>=size)?0:start;
return s;
}
//+------------------------------------------------------------------+
//|process stop parameter |
//+------------------------------------------------------------------+
long normalizestop(long stop, ulong size)
{
long s = (stop>0 && ulong(stop)>=size)?long(size):(stop<0 && ulong(fabs(stop))>size)?-1:(stop<0 && ulong(fabs(stop))<=size)?long(size)+stop:stop;
return s;
}
//+------------------------------------------------------------------+
//| validate all calculated params |
//+------------------------------------------------------------------+
bool invalidparams(long start,long stop, long step, ulong size)
{
if(step==0 || start==stop || ulong(start)>=size || start<0 || ulong(fabs(stop))>size || (step<0 && start<stop) || (step>0 && start>stop))
return true;
return false;
}
//+------------------------------------------------------------------+
//| quantiles helper function |
//+------------------------------------------------------------------+
vector quantile(vector &in,vector &abprobs,vector &probs)
{
vector x = in;
if(!sort(x))
{
Print(__FUNCTION__, " sort failure ");
vector::Zeros(probs.Size());
}
ulong n = x.Size();
if(!x.Size())
{
Print(__FUNCTION__, " invalid parameter, empty vector ");
vector::Zeros(probs.Size());
}
vector aleph = (double(n)*probs+abprobs);
vector k = aleph;
if(!k.Clip(1.0, double(n-1)))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(probs.Size());
}
k = floor(k);
vector gamma = aleph-k;
if(!gamma.Clip(0.0, 1.0))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(probs.Size());
}
vector k1 = k-1.0;
ulong index[];
if(!np::vecAsArray(k1,index))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(probs.Size());
}
vector x1 = np::select(k1,index);
ArrayFree(index);
if(!np::vecAsArray(k,index))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(probs.Size());
}
vector x2 = np::select(k,index);
return (1.0-gamma)*x1+gamma*x2;
}
}
//+------------------------------------------------------------------+
//| Return arrays representing indices of a grid |
//+------------------------------------------------------------------+
void indices(ulong rows, ulong cols, matrix &outRC[])
{
ulong N = 2;
if(ArrayResize(outRC,2)<0 || !outRC[0].Resize(rows,cols) || !outRC[1].Resize(rows,cols))
{
Print(__FUNCTION__, " resize error ");
return;
}
vector v = np::arange(rows);
outRC[0] = np::repeat_vector_as_rows_cols(v,cols,false);
v = np::arange(cols);
outRC[1] = np::repeat_vector_as_rows_cols(v,rows,true);
}
//+------------------------------------------------------------------+
//| The function sorts (in) vector and indices[] simultaneously |
//| using QuickSort algorithm. After sorting the initial ordering |
//| of the elements is located in the indices[] array. |
//| |
//| Arguments: |
//| in : vector with values to sort |
//| ascending : direction of sort (true:ascending) |
//| (false:descending) |
//| first : First element index |
//| last : Last element index |
//| |
//| Return value: None |
//+------------------------------------------------------------------+
template <typename T>
bool quickSort(vector<T> &in,bool ascending,long first,long last)
{
long i,j;
T p_double,t_double;
//--- check
if(first<0 || last<0)
{
return false;
}
//--- sort
i=first;
j=last;
while(i<last)
{
//--- ">>1" is quick division by 2
p_double=in[(first+last)>>1];
while(i<j)
{
while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
{
//--- control the output of the in bounds
if(i==in.Size()-1)
break;
i++;
}
while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
{
//--- control the output of the in bounds
if(j==0)
break;
j--;
}
if(i<=j)
{
//-- swap elements i and j
t_double=in[i];
in[i]=in[j];
in[j]=t_double;
//--
i++;
//--- control the output of the in bounds
if(j==0)
break;
j--;
}
}
if(first<j)
quickSort(in,ascending,first,j);
first=i;
j=last;
}
return true;
}
//+------------------------------------------------------------------+
//| The function sorts (in) vector and indices[] simultaneously |
//| using QuickSort algorithm. After sorting the initial ordering |
//| of the elements is located in the indices[] array. |
//| |
//| Arguments: |
//| in : vector with values to sort |
//| ascending : direction of sort (true:ascending) |
//| (false:descending) |
//| indices : array of indexes |
//| first : First element index |
//| last : Last element index |
//| |
//| Return value: None |
//+------------------------------------------------------------------+
template <typename T>
bool quickSortIndices(vector<T>& _in,bool ascending,long &indices[],long first,long last)
{
vector<T> in = _in;
long i,j,t_int;
T p_double,t_double;
//--- check
if(first<0 || last<0)
{
return false;
}
//--- sort
i=first;
j=last;
while(i<last)
{
//--- ">>1" is quick division by 2
p_double=in[(first+last)>>1];
while(i<j)
{
while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
{
//--- control the output of the in bounds
if(i==in.Size()-1)
break;
i++;
}
while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
{
//--- control the output of the in bounds
if(j==0)
break;
j--;
}
if(i<=j)
{
//-- swap elements i and j
t_double=in[i];
in[i]=in[j];
in[j]=t_double;
//-- swap indices i and j
t_int=indices[i];
indices[i]=indices[j];
indices[j]=t_int;
i++;
//--- control the output of the in bounds
if(j==0)
break;
j--;
}
}
if(first<j)
quickSortIndices(in,ascending,indices,first,j);
first=i;
j=last;
}
return true;
}
//+------------------------------------------------------------------+
//| sort a vector |
//+------------------------------------------------------------------+
template<typename T>
bool sort(vector<T>&in, bool ascending=true)
{
if(in.Size()<1)
return false;
return quickSort(in,ascending,0,in.Size()-1);
}
//+------------------------------------------------------------------+
//| sort matrix by columns or rows |
//+------------------------------------------------------------------+
template<typename T>
bool sort(matrix<T>&in, bool by_rows, bool ascending=true)
{
if(by_rows)
{
for(ulong i = 0; i<in.Rows(); ++i)
{
vector row = in.Row(i);
if(!sort(row,ascending))
return false;
in.Row(row,i);
}
}
else
{
for(ulong i = 0; i<in.Cols(); ++i)
{
vector col = in.Col(i);
if(!sort(col,ascending))
return false;
in.Col(col,i);
}
}
return true;
}
//+------------------------------------------------------------------+
//| extract a portion of vector within a range at specific interval |
//| with support for negative indexing python style |
//| returns new vector from start to stop-1 at step intervals |
//| start - index of first element (inclusive) |
//| stop - index past last element (exclusive) |
//| step - step size |
//+------------------------------------------------------------------+
template<typename T>
vector<T> sliceVector(vector<T> &in, const long index_start=BEGIN, const long index_stop = END, const long index_step = STEP)
{
//---local variables
long st,sp,stp,size;
//---check function parameters
st = internal::normalizestart(index_start,in.Size());
sp = internal::normalizestop(index_stop,in.Size());
stp = index_step;
//---error handling
if(internal::invalidparams(st,sp,stp,in.Size()))
{
Print(__FUNCTION__," invalid function parameters ");
return vector<T>::Zeros(0);
}
//---calculate size of new vector
size = internal::calculatesize(st,sp,stp);
//---initialize output vector
vector<T> out(ulong(size));
//---assign elements to output
for(long pos = st, index=0; index<size; ++index,pos+=stp)
out[index] = in[pos];
//---done
return out;
}
//+------------------------------------------------------------------+
//| set subset of contiguous container elements to a certain value |
//+------------------------------------------------------------------+
template<typename T>
bool fillVector(vector<T> &in, T value, const long index_start=BEGIN, const long index_stop = END, const long index_step = STEP)
{
//---local variables
long st,sp,stp,size;
//---check function parameters
st = internal::normalizestart(index_start,in.Size());
sp = internal::normalizestop(index_stop,in.Size());
stp = index_step;
//---error handling
if(internal::invalidparams(st,sp,stp,in.Size()))
{
Print(__FUNCTION__," invalid function parameters ");
return false;
}
//---calculate size of new vector
size = internal::calculatesize(st,sp,stp);
//---assign elements to output
for(long pos = st, index=0; index<size; ++index,pos+=stp)
in[pos] = value;
//---done
return true;
}
//+------------------------------------------------------------------+
//| reverse a vector |
//+------------------------------------------------------------------+
template<typename T>
bool reverseVector(vector<T> &in)
{
if(in.Size()<2)
return in.Size()!=0;
ulong size = in.Size()-1;
ulong half = (in.Size())>>1;
for(ulong i = 0; i<half; ++i)
{
T temp = in[i];
in[i] = in[size - i];
in[size - i] = temp;
}
return true;
}
//+------------------------------------------------------------------+
//| extract a portion a matrix within a range at specific intervals |
//| with support for negative indexing python style |
//| returns new matrix |
//| row_start - row index of first row (inclusive) |
//| row_stop - row index past last row (exclusive) |
//| row_step - row step size |
//| column_start - columnindex of first column (inclusive) |
//| column_stop - column index past last columun (exclusive) |
//| column_step - column step size |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrix(matrix<T> &in, const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP,const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
{
//---local variables
long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
rst =internal::normalizestart(row_start,in.Rows());
rsp =internal::normalizestop(row_stop,in.Rows());
rstp = row_step;
//---error handling
if(internal::invalidparams(rst,rsp,rstp,in.Rows()))
{
Print(__FUNCTION__," invalid function parameters for row specification ");
return matrix<T>::Zeros(0,0);
}
//---calculate size of new matrix rows
rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
cst = internal::normalizestart(column_start,in.Cols());
csp = internal::normalizestop(column_stop,in.Cols());
cstp = column_step;
//---error handling
if(internal::invalidparams(cst,csp,cstp,in.Cols()))
{
Print(__FUNCTION__," invalid function parameters for column specification ");
return matrix<T>::Zeros(0,0);
}
//---calculate size of new matrix columns
csize = internal::calculatesize(cst,csp,cstp);
//---initialize output matrix
matrix<T> out(ulong(rsize),ulong(csize));
//---assign elements to output
for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
out[row_index][column_index] = in[rpos][cpos];
//---done
return out;
}
//+------------------------------------------------------------------+
//| fill a portion a matrix within a range at specific intervals |
//| with support for negative indexing python style |
//| edits container in place |
//| row_start - row index of first row (inclusive) |
//| row_stop - row index past last row (exclusive) |
//| row_step - row step size |
//| column_start - columnindex of first column (inclusive) |
//| column_stop - column index past last columun (exclusive) |
//| column_step - column step size |
//+------------------------------------------------------------------+
template<typename T>
bool fillMatrix(matrix<T> &in, const T value,const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP,const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
{
//---local variables
long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
rst =internal::normalizestart(row_start,in.Rows());
rsp =internal::normalizestop(row_stop,in.Rows());
rstp = row_step;
//---error handling
if(internal::invalidparams(rst,rsp,rstp,in.Rows()))
{
Print(__FUNCTION__," invalid function parameters for row specification ");
return false;
}
//---calculate size of new matrix rows
rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
cst = internal::normalizestart(column_start,in.Cols());
csp = internal::normalizestop(column_stop,in.Cols());
cstp = column_step;
//---error handling
if(internal::invalidparams(cst,csp,cstp,in.Cols()))
{
Print(__FUNCTION__," invalid function parameters for column specification ");
return false;
}
//---calculate size of new matrix columns
csize = internal::calculatesize(cst,csp,cstp);
//---assign elements to output
for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
in[rpos][cpos] = value;
//---done
return true;
}
//+------------------------------------------------------------------+
//| extract a portion of matrix by selecting range of columns and |
//| optional interval |
//| with support for negative indexing python style |
//| returns new matrix |
//| column_start - columnindex of first column (inclusive) |
//| column_stop - column index past last columun (exclusive) |
//| column_step - column step size |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrixCols(matrix<T> &in, const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
{
return sliceMatrix(in,BEGIN,END,STEP,column_start,column_stop,column_step);
}
//+------------------------------------------------------------------+
//|reverse columns of matrix |
//+------------------------------------------------------------------+
template<typename T>
bool reverseMatrixCols(matrix<T> &in)
{
ulong cols = in.Cols()-1;
if(in.Cols()<2)
return in.Cols()!=0;
ulong cols_end = (in.Cols())>>1;
for(ulong i = 0; i<cols_end; ++i)
if(!in.SwapCols(i,cols-i))
return false;
return true;
}
//+------------------------------------------------------------------+
//| extract a portion of matrix by selecting range of rows and |
//| optional interval |
//| with support for negative indexing python style |
//| returns new matrix |
//| row_start - row index of first row (inclusive) |
//| row_stop - row index past last row (exclusive) |
//| row_step - row step size |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrixRows(matrix<T> &in, const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP)
{
return sliceMatrix(in,row_start,row_stop,row_step);
}
//+------------------------------------------------------------------+
//|reverse rows of matrix |
//+------------------------------------------------------------------+
template<typename T>
bool reverseMatrixRows(matrix<T> &in)
{
ulong rows = in.Rows()-1;
if(in.Rows()<2)
return in.Rows()!=0;
ulong rows_end = (in.Rows())>>1;
for(ulong i = 0; i<rows_end; ++i)
if(!in.SwapRows(i,rows-i))
return false;
return true;
}
//+------------------------------------------------------------------+
//| assign values to a portion of matrix selected within range |
//| of columns and |
//| optional interval |
//| with support for negative indexing python style |
//| dest_start - column index of first column (inclusive) |
//| dest_stop - column index past last columun (exclusive) |
//| dest_step - column step size |
//| indexes refer to positions in copyto |
//+------------------------------------------------------------------+
template <typename T>
bool matrixCopyCols(matrix<T> &copyto,matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
{
return matrixCopy(copyto,copyfrom,BEGIN,END,STEP,dest_start,dest_stop,dest_step);
}
//+------------------------------------------------------------------+
//| assign values to a portion of matrix selected within range |
//| of rows and |
//| optional interval |
//| with support for negative indexing python style |
//| dest_start - row index of first row (inclusive) |
//| dest_stop - row index past last row (exclusive) |
//| dest_step - row step size |
//| indexes refer to positions in copyto |
//+------------------------------------------------------------------+
template <typename T>
bool matrixCopyRows(matrix<T> &copyto, matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
{
return matrixCopy(copyto,copyfrom,dest_start,dest_stop,dest_step);
}
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as collection |
//| of indices |
//+------------------------------------------------------------------+
template <typename T, typename P>
vector<T> select(vector<T> &in,const P &indices[])
{
P size = (P)MathMin(in.Size(),indices.Size());
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=P(in.Size()))
{
Print(__FUNCTION__, " invalid function parameter ");
return vector<T>::Zeros(1);
}
vector<T> out(size);
for(P i=0; i<size; ++i)
out[i] = in[indices[i]];
return out;
}
//+------------------------------------------------------------------+
//|return arbitrary columns of a matrix specified as |
//|collection of indices |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixCols(matrix<T> &in,const long &indices[])
{
ulong size = (ulong)MathMin(in.Cols(),indices.Size());
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=long(in.Cols()))
{
Print(__FUNCTION__," invalid function parameters ");
return matrix<T>::Zeros(0,0);
}
matrix<T> out(in.Rows(),size);
for(ulong i=0; i<size; ++i)
out.Col(in.Col(indices[i]),i);
return out;
}
//+------------------------------------------------------------------+
//|return arbitrary rows of a matrix specified as |
//|collection of indices |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixRows(matrix<T> &in,const long &indices[])
{
ulong size = (ulong)MathMin(in.Rows(),indices.Size());
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=long(in.Rows()))
{
Print(__FUNCTION__," invalid function parameters ");
return matrix<T>::Zeros(0,0);
}
matrix<T> out(size,in.Cols());
for(ulong i=0; i<size; ++i)
out.Row(in.Row(indices[i]),i);
return out;
}
//+------------------------------------------------------------------+
//|copy vector copyfrom into specified positions of vector copyto |
//|with support for negative indexing python style |
//+------------------------------------------------------------------+
template <typename T>
bool vectorCopy(vector<T> &copyto,vector<T> &copyfrom,const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
{
//---local variables
long st,sp,stp,size;
//---check function parameters
st = internal::normalizestart(dest_start,copyto.Size());
sp = internal::normalizestop(dest_stop,copyto.Size());
stp = dest_step;
//---error handling
if(internal::invalidparams(st,sp,stp,copyto.Size()))
{
Print(__FUNCTION__," invalid function parameters ");
return false;
}
//---calculate size of new vector
size = internal::calculatesize(st,sp,stp);
//---error handling
if(copyfrom.Size()<ulong(size))
{
Print(__FUNCTION__," source vector is of insufficient size ");
return false;
}
//---assign elements
for(long pos = st,index=0; index<size; ++index,pos+=stp)
copyto[pos] = copyfrom[index];
//---done
return true;
}
//+------------------------------------------------------------------+
//| assign selected elements of vector the value /fill_value/ |
//+------------------------------------------------------------------+
template <typename T>
bool vectorFill(vector<T> &copyto, T fill_value,const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
{
//---local variables
long st,sp,stp,size;
//---check function parameters
st = internal::normalizestart(dest_start,copyto.Size());
sp = internal::normalizestop(dest_stop,copyto.Size());
stp = dest_step;
//---error handling
if(internal::invalidparams(st,sp,stp,copyto.Size()))
{
Print(__FUNCTION__," invalid function parameters ");
return false;
}
//---calculate size of new vector
size = internal::calculatesize(st,sp,stp);
//---assign elements
for(long pos = st,index=0; index<size; ++index,pos+=stp)
copyto[pos] = fill_value;
//---done
return true;
}
//+------------------------------------------------------------------+
//|copy matrix copyfrom into specified positions of matrix copyto |
//|with support for negative indexing python style |
//+------------------------------------------------------------------+
template<typename T>
bool matrixCopy(matrix<T> &copyto,matrix &copyfrom, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP,const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
{
//---local variables
long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
rst =internal::normalizestart(row_dest_start,copyto.Rows());
rsp =internal::normalizestop(row_dest_stop,copyto.Rows());
rstp = row_dest_step;
//---error handling
if(internal::invalidparams(rst,rsp,rstp,copyto.Rows()))
{
Print(__FUNCTION__," invalid function parameters for row specification ");
return false;
}
//---calculate size of new matrix rows
rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
cst = internal::normalizestart(column_dest_start,copyto.Cols());
csp = internal::normalizestop(column_dest_stop,copyto.Cols());
cstp = column_dest_step;
//---error handling
if(internal::invalidparams(cst,csp,cstp,copyto.Cols()))
{
Print(__FUNCTION__," invalid function parameters for column specification ");
return false;
}
//---calculate size of new matrix columns
csize = internal::calculatesize(cst,csp,cstp);
//---error handling
if(copyfrom.Cols() < ulong(csize) || copyfrom.Rows() < ulong(rsize))
{
Print(__FUNCTION__, " source matrix has insufficient number of columns and/or rows ");
return false;
}
//---assign elements
for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
copyto[rpos][cpos]=copyfrom[row_index][column_index];
//---done
return true;
}
//+------------------------------------------------------------------+
//| assign selected elements of vector the value /fill_value/ |
//+------------------------------------------------------------------+
template<typename T>
bool matrixFill(matrix<T> &copyto,T fill_value, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP,const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
{
//---local variables
long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
rst =internal::normalizestart(row_dest_start,copyto.Rows());
rsp =internal::normalizestop(row_dest_stop,copyto.Rows());
rstp = row_dest_step;
//---error handling
if(internal::invalidparams(rst,rsp,rstp,copyto.Rows()))
{
Print(__FUNCTION__," invalid function parameters for row specification ");
return false;
}
//---calculate size of new matrix rows
rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
cst = internal::normalizestart(column_dest_start,copyto.Cols());
csp = internal::normalizestop(column_dest_stop,copyto.Cols());
cstp = column_dest_step;
//---error handling
if(internal::invalidparams(cst,csp,cstp,copyto.Cols()))
{
Print(__FUNCTION__," invalid function parameters for column specification ");
return false;
}
//---calculate size of new matrix columns
csize = internal::calculatesize(cst,csp,cstp);
//---assign elements
for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
copyto[rpos][cpos]=fill_value;
//---done
return true;
}
//+------------------------------------------------------------------+
//| returns vector |
//+------------------------------------------------------------------+
vector arange(const ulong size, double value = 0.0, double step_value = 1.0)
{
vector out(size);
for(ulong i = 0; i<size; ++i, value+=step_value)
out[i] = value;
return out;
}
//+------------------------------------------------------------------+
//| returns vector |
//+------------------------------------------------------------------+
vector range(int stop_value, int start_value = 0, int step_value = 1)
{
double size = ceil(double(stop_value - start_value)/double(step_value));
vector out(ulong(size));
double value = double(start_value);
double step = double(step_value);
for(ulong i = 0; i<ulong(size); ++i, value+=step)
out[i] = value;
return out;
}
//+------------------------------------------------------------------+
//| outputs array |
//+------------------------------------------------------------------+
template<typename T>
bool arange(T &in_out[],const int size, T value = 0, T step_value = 1)
{
if(!size)
return false;
if(ArrayResize(in_out,fabs(size))!=fabs(size))
{
Print(__FUNCTION__, " error ", GetLastError());
return false;
}
for(int i = 0; i<size; ++i, value+=step_value)
in_out[i] = value;
return true;
}
//+------------------------------------------------------------------+
//| equivalent to numpy linspace(), outputs vector |
//+------------------------------------------------------------------+
template<typename T>
vector<T> linspace(T start, T stop, ulong num=50)
{
vector<T> out(num);
T chunk = (stop-start)/T(num-1);
for(ulong i = 0; i<out.Size(); ++i)
out[i] = start+(T(i)*chunk);
return out;
}
//+------------------------------------------------------------------+
//| fill diagonal |
//+------------------------------------------------------------------+
template<typename T>
bool fillDiagonal(matrix<T> &in_out, T fill_value)
{
if(in_out.Rows()!=in_out.Cols())
{
Print(__FUNCTION__ "Invalid Input. Matrix is not square ");
return false;
}
for(ulong i = 0; i<in_out.Rows(); ++i)
in_out[i][i] = fill_value;
return true;
}
//+------------------------------------------------------------------+
//| copy 3d container |
//+------------------------------------------------------------------+
bool copy3D(matrix &copy_from[], matrix &copy_to[])
{
if(ArraySize(copy_to) != ArraySize(copy_from))
ArrayResize(copy_to,(int)copy_from.Size());
for(uint i = 0; i<copy_to.Size(); ++i)
copy_to[i] = copy_from[i];
return true;
}
//+-----------------------------------------------------------------------+
//| Generate a Vandermonde matrix |
//| in : vector |
//| num_columns : ulong, default: size of in input vector. |
//| ascending : bool, (default: false) Order of the powers of the columns.|
//| If true, the powers increase |
//| from left to right, if false (the default) they are reversed |
//+-----------------------------------------------------------------------+
template<typename T>
matrix<T> vander(vector<T> &in, ulong num_columns = 0,bool ascending = false)
{
ulong n = num_columns?num_columns:in.Size();
matrix<T> out(in.Size(),n);
for(ulong i = 0; i<n; ++i)
out.Col(pow(in,ascending?i:n-1-i),i);
return out;
}
//+------------------------------------------------------------------+
//|generic vector to array |
//+------------------------------------------------------------------+
template<typename T,typename P>
bool vecAsArray(const vector<T> &in, P &out[])
{
//---
if(in.Size()<1)
{
Print(__FUNCTION__," Empty vector");
return false;
}
//---
if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
{
Print(__FUNCTION__," resize error ", GetLastError());
return false;
}
//---
for(uint i = 0; i<out.Size(); ++i)
out[i] = P(in[i]);
//---
return true;
//---
}
//+------------------------------------------------------------------+
//|generic matrix to 1-D array |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, S &out[])
{
//---
if(in.Cols()+in.Rows()<1)
{
Print(__FUNCTION__," Empty vector");
return false;
}
//---
if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
{
Print(__FUNCTION__," resize error ", GetLastError());
return false;
}
//---
for(uint i = 0; i<out.Size(); ++i)
out[i]=S(in.Flat(i));
//---
return true;
//---
}
//+------------------------------------------------------------------+
//| difference vector according to parameter degreee |
//+------------------------------------------------------------------+
template<typename T>
vector<T> diff(vector<T> &in, ulong degree = 1)
{
//---
if(in.Size()<1 || degree<1 || in.Size()<degree+1)
{
Print(__FUNCTION__," Invalid parameter ");
return vector<T>::Zeros(0);
}
//---
vector<T> yy,zz;
//---
yy = sliceVector(in,long(degree));
zz = sliceVector(in,BEGIN,long(in.Size()-degree));
//---
return yy-zz;
}
//+------------------------------------------------------------------+
//| difference matrix about the rows or columns |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> diff(matrix<T> &in, ulong degree = 1, bool by_row = true)
{
//---
if(in.Cols()+in.Rows()<1 || degree<1 || (by_row && in.Cols()<degree+1) || (!by_row && in.Rows()<degree+1))
{
Print(__FUNCTION__," Invalid function parameter ");
return matrix::Zeros(0,0);
}
//---
matrix<T> yy,zz;
//---
if(by_row)
{
yy = sliceMatrix(in,BEGIN,END,STEP,long(degree));
zz = sliceMatrix(in,BEGIN,END,STEP,BEGIN,long(in.Cols()-degree));
}
else
{
yy = sliceMatrix(in,long(degree));
zz = sliceMatrix(in,BEGIN,long(in.Rows()-degree));
}
//---
return yy-zz;
}
//+------------------------------------------------------------------+
//| generate vector of num_repetitions of input vector elements |
//+------------------------------------------------------------------+
template<typename T>
vector<T> repeat_vector_elements(vector<T> &vect,ulong num_repetitions)
{
vector<T> out(num_repetitions*vect.Size());
for(ulong i = 0; i<out.Size(); i+=vect.Size())
vectorCopy(out,vect,i,i+vect.Size());
return out;
}
//+------------------------------------------------------------------+
//| generate vector of num_repetitions of input matrix elements |
//+------------------------------------------------------------------+
template<typename T>
vector<T> repeat_matrix_elements(matrix<T> &in, ulong num_repetitions)
{
vector<T> out(num_repetitions*in.Cols()*in.Rows());
for(ulong i = 0; i<out.Size(); i+=(in.Cols()*num_repetitions))
{
vector copy = repeat_vector_elements(in.Row(i/(in.Cols()*num_repetitions)),num_repetitions);
vectorCopy(out,copy,i,i+copy.Size());
}
return out;
}
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as rows or columns |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_vector_as_rows_cols(vector<T>&in, ulong num_repetitions, bool as_rows = true)
{
matrix<T> out(as_rows?num_repetitions:in.Size(),as_rows?in.Size():num_repetitions);
for(ulong i = 0; i<num_repetitions; ++i)
as_rows?out.Row(in,i):out.Col(in,i);
return out;
}
//+------------------------------------------------------------------+
//| generate matrix of num_repetitions of input columns or rows |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_matrix_rows_cols(matrix<T> &in, ulong num_repetitions, bool by_rows)
{
if(num_repetitions<1)
{
Print(__FUNCTION__, " invalid function parameter ");
return matrix<T>::Zeros(0,0);
}
matrix<T> out((by_rows)?num_repetitions*in.Rows():in.Rows(),(!by_rows)?num_repetitions*in.Cols():in.Cols());
if(by_rows)
for(ulong i = 0; i<out.Rows(); ++i)
out.Row(in.Row(i/num_repetitions),i);
else
for(ulong i = 0; i<out.Cols(); ++i)
out.Col(in.Col(i/num_repetitions),i);
return out;
}
//+------------------------------------------------------------------+
//| generate matrix of an arbitrary number of repetitions of each row|
//| or column as specified by elements of input array |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_matrix_rows_cols(matrix<T> &in, ulong &reps[],bool by_rows)
{
ulong max = 0 ;
for(uint i = 0; i<reps.Size(); ++i)
max+=reps[i];
if((by_rows && in.Rows()!=ulong(reps.Size())) ||
(!by_rows && in.Cols()!=ulong(reps.Size()))||
reps[ArrayMinimum(reps)]<1)
{
Print(__FUNCTION__, " invalid function parameters ");
return matrix<T>::Zeros(0,0);
}
matrix<T> out((by_rows)?max:in.Rows(),(!by_rows)?max:in.Cols());
if(by_rows)
for(ulong i = 0, j = 0; i<out.Rows(); i+=reps[j],++j)
{
matrix sliced =sliceRows(in,j,j+1);
matrix rpl = repeat_matrix_rows_cols(sliced,reps[j],true);
copyRows(out,rpl,i,i+reps[j]);
}
else
for(ulong i = 0, j = 0; i<out.Cols(); i+=reps[j],++j)
{
matrix sliced =sliceCols(in,j,j+1);
matrix rpl = repeat_matrix_rows_cols(sliced,reps[j],false);
copyCols(out,rpl,i,i+reps[j]);
}
return out;
}
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as columns |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> vectorAsColumnMatrix(vector<T>&in, ulong num_cols)
{
return repeat_vector_as_rows_cols(in,num_cols,false);
}
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as rows |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> vectorAsRowMatrix(vector<T>&in, ulong num_rows)
{
return repeat_vector_as_rows_cols(in,num_rows,true);
}
//+------------------------------------------------------------------+
//| The function extracts the unique values from the vector. |
//| |
//| Arguments: |
//| in : vector with values |
//| |
//| |
//| Return value: out vector of unique values . |
//+------------------------------------------------------------------+
template<typename T>
vector<T> unique(vector<T> &in)
{
//--- check array size
ulong size=in.Size();
if(size==0)
{
Print(__FUNCTION__," Invalid function parameter ");
return vector<T>::Zeros(1);
}
//--- prepare additional in
bool checked[];
if(ArrayResize(checked,int(size))!=int(size))
{
Print(__FUNCTION__," error ",GetLastError());
return vector<T>::Zeros(0) ;
}
ArrayFill(checked,0,int(size),false);
//--- prepare out in
vector<T>out(size);
//--- find unique elements
ulong unique_count=0;
T value=0;
while(true)
{
bool flag=false;
for(ulong i=unique_count; i<size; ++i)
{
if(!flag && !checked[i])
{
value=in[i];
out[unique_count]=in[i];
++unique_count;
checked[i]=true;
flag=true;
}
else
if(flag && value==in[i])
checked[i]=true;
}
if(!flag)
break;
}
//--- resize target in
out.Resize(unique_count);
//---
return(out);
}
//+------------------------------------------------------------------+
//| The function extracts the unique values from the vector. |
//| |
//| Arguments: |
//| in : vector with values |
//| out : vector of unique values |
//| indices: index positions in original vector |
//| counts : number of occurances of each unique value |
//| Return bool: true on success and false on error |
//+------------------------------------------------------------------+
template<typename T>
bool unique(vector<T> &in, vector<T> &out, long &indices[],long &counts[])
{
//--- check vector size
ulong size=in.Size();
if(size==0)
{
Print(__FUNCTION__," Invalid function parameter ");
return false;
}
//--- prepare internal and output buffers
bool checked[];
if(ArrayResize(checked,int(size))!=int(size) ||
ArrayResize(indices,int(size))!=int(size) ||
ArrayResize(counts,int(size))!=int(size) ||
!out.Resize(size))
{
Print(__FUNCTION__, " error ", GetLastError());
return(false);
}
ArrayFill(checked,0,int(size),false);
ArrayFill(indices,0,int(size),-1);
ArrayFill(counts,0,int(size),1);
//--- prepare out vector
out.Resize(size);
//--- find unique elements
ulong unique_count=0;
T value=0;
while(true)
{
bool flag=false;
for(ulong i=unique_count; i<size; ++i)
{
if(!flag && !checked[i])
{
value=in[i];
out[unique_count]=in[i];
indices[unique_count]=long(i);
++unique_count;
checked[i]=true;
flag=true;
}
else
if(flag && value==in[i])
{
checked[i]=true;
++counts[unique_count-1];
}
}
if(!flag)
break;
}
//---
return (out.Resize(unique_count) && ArrayResize(indices,int(unique_count))==int(unique_count) && ArrayResize(counts,int(unique_count))==int(unique_count));
//---
}
//+------------------------------------------------------------------+
//| shuffle elements of vector |
//+------------------------------------------------------------------+
template<typename T>
vector shuffleVector(vector<T> &in,CHighQualityRandStateShell* rstate=NULL)
{
vector<T>out(in.Size());
bool destroy = false;
//---
if(rstate == NULL)
{
destroy = true;
rstate = new CHighQualityRandStateShell();
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
}
//---
ulong k=0;
ulong size=in.Size();
out=in;
int error = 0;
T tmp;
//---
for(ulong i=0; i<size; ++i)
{
k=(int)(CAlglib::HQRndUniformR(rstate)*size);
//---
if(k>=size)
k=size-1;
//---
tmp = out[i];
out[i]=out[k];
out[k]=tmp;
}
//---
if(destroy)
delete rstate;
//---
return out;
}
//+------------------------------------------------------------------+
//| shuffle elements of vector |
//+------------------------------------------------------------------+
template<typename T>
bool shuffleArray(T &in[],CHighQualityRandStateShell* rstate=NULL)
{
bool destroy = false;
//---
if(rstate == NULL)
{
destroy = true;
rstate = new CHighQualityRandStateShell();
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
}
//---
ulong k=0;
uint size=in.Size();
int error = 0;
T tmp;
//---
for(ulong i=0; i<size; ++i)
{
k=(int)(CAlglib::HQRndUniformR(rstate)*size);
//---
if(k>=size)
k=size-1;
//---
tmp = in[i];
in[i]=in[k];
in[k]=tmp;
}
//---
if(destroy)
delete rstate;
//---
return true;
}
//+------------------------------------------------------------------+
//| shuffle elements of vector |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> shuffleMatrix(matrix<T> &in,bool by_rows=true,CHighQualityRandStateShell* rstate=NULL)
{
matrix<T>out(in.Rows(), in.Cols());
bool destroy = false;
//---
if(rstate == NULL)
{
destroy = true;
rstate = new CHighQualityRandStateShell();
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
}
//---
ulong size=by_rows?in.Rows():in.Cols();
vector<T> vec;
out=in;
//---
for(ulong i=0; i<size; ++i)
{
vec = by_rows?in.Row(i):in.Col(i);
vec = shuffleVector(vec,rstate);
if((by_rows && !out.Row(vec,i)) ||
(!by_rows && !out.Col(vec,i)))
{
Print(__FUNCTION__, " error ", GetLastError());
return matrix::Zeros(size,size);
}
}
//---
if(destroy)
delete rstate;
//---
return out;
}
//+------------------------------------------------------------------+
//|Integrate using the composite trapezoidal rule |
//+------------------------------------------------------------------+
template<typename T>
T trapezoid(vector<T>&y, vector<T>&x,T dx=1)
{
vector<T>d;
if(x.Size()==0)
{
d = vector::Zeros(y.Size());
d.Fill(dx);
}
else
d = diff(x);
vector<T> ret = (d*(sliceVector(y,1) + sliceVector(y,0,-1)))/2.0;
return ret.Sum();
}
//+---------------------------------------------------------------------------------+
//| One-dimensional linear interpolation for monotonically increasing sample points.|
//+---------------------------------------------------------------------------------+
template<typename T>
vector<T>interp(vector<T> &xv,vector<T>& xpoints,vector<T> &ypoints,double left=DBL_MIN,double right=DBL_MAX, double period=EMPTY_VALUE)
{
vector<T>out = vector<T>::Zeros(xv.Size());
if(period != EMPTY_VALUE)
{
if(period == 0.0)
{
Print(__FUNCTION__," period should be non zero");
return out;
}
vector<T> x = xv;
vector<T> xp = xpoints;
vector<T> yp = ypoints;
vector vp(x.Size());
vp.Fill(MathAbs(period));
x = MathMod(xv,vp);
xp = MathMod(xpoints,vp);
long indices[];
if(!arange(indices,int(xp.Size())) || !quickSortIndices(xp,true,indices,0,long(xp.Size()-1)))
{
Print(__FUNCTION__, " ", __LINE__);
return out;
}
xp = select(xp,indices);
yp = select(yp,indices);
if(!xpoints.Resize((xp.Size()*2)+1) || !ypoints.Resize((yp.Size()*2)+1))
{
Print(__FUNCTION__, " ", __LINE__, " error ", GetLastError());
return out;
}
if(!vectorCopy(xpoints,xp,xp.Size(),xp.Size()*2) || !vectorCopy(ypoints,yp,yp.Size(),yp.Size()*2))
{
Print(__FUNCTION__, " ", __LINE__);
return out;
}
xpoints[xpoints.Size()-1] = xp[0]+MathAbs(period);
ypoints[ypoints.Size()-1] = yp[0];
if(!reverseVector(xp))
{
Print(__FUNCTION__, " ", __LINE__);
return out;
}
vector temp = xp-vp;
if(!vectorCopy(xpoints,temp,0,temp.Size()))
{
Print(__FUNCTION__, " ", __LINE__);
return out;
}
if(!reverseVector(yp) && !vectorCopy(ypoints,yp,0,yp.Size()))
{
Print(__FUNCTION__, " ", __LINE__);
return out;
}
xv = x;
}
double lval,rval;
double xvals[];
if(!vecAsArray(xpoints,xvals) || !ArraySort(xvals))
{
Print(__FUNCTION__, " ", __LINE__," error ", GetLastError());
return out;
}
lval=(left==DBL_MIN)?ypoints[0]:left;
rval=(right==DBL_MAX)?ypoints[ypoints.Size()-1]:right;
if(xpoints.Size()==1)
{
for(ulong i =0; i<out.Size(); ++i)
out[i] = xv[i]<xpoints[0]?lval:(xv[i]>xpoints[0])?rval:ypoints[0];
}
else
{
vector slopes=vector::Zeros(0);
if(xpoints.Size()<=xv.Size())
slopes = vector::Zeros(xpoints.Size()-1);
if(slopes.Size())
{
vector ypd = sliceVector(ypoints,1) - sliceVector(ypoints,0,ypoints.Size()-1);
vector xpd = sliceVector(xpoints,1) - sliceVector(xpoints,0,xpoints.Size()-1);
slopes = ypd/xpd;
}
int j = 0;
for(ulong i = 0; i<xv.Size(); ++i)
{
if(MathClassify(xv[i])==FP_NAN)
out[i] = xv[i];
j = internal::binary_search(xv[i],xvals);
if(j == -1)
out[i] = lval;
else
if(j == int(xpoints.Size()))
out[i] = rval;
else
if(j == int(xpoints.Size()-1))
out[i] = ypoints[j];
else
if(xpoints[j] == xv[i])
out[i] = ypoints[j];
else
{
double slope = slopes.Size()?slopes[j]:(ypoints[j+1] - ypoints[j]) / (xpoints[j+1] - xpoints[j]);
out[i] = slope*(xv[i] - xpoints[j])+ypoints[j];
if(MathClassify(out[i])==FP_NAN)
{
out[i] = slope*(xv[i] - xpoints[j+1]) + ypoints[j+1];
if(MathClassify(out[i])==FP_NAN && ypoints[j] == ypoints[j+1])
out[i] = ypoints[j];
}
}
}
}
return out;
}
//+------------------------------------------------------------------+
//| matrix multiplication |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> matmul(const matrix<T>& matrix_a, const matrix<T>& matrix_b)
{
matrix<T> matrix_c = matrix<T>::Zeros(0,0);
if(matrix_a.Cols()!=matrix_b.Rows())
return(matrix_c);
ulong M=matrix_a.Rows();
ulong K=matrix_a.Cols();
ulong N=matrix_b.Cols();
matrix_c=matrix<T>::Zeros(M,N);
for(ulong m=0; m<M; ++m)
for(ulong k=0; k<K; ++k)
for(ulong n=0; n<N; ++n)
matrix_c[m][n]+=matrix_a[m][k]*matrix_b[k][n];
return(matrix_c);
}
//+------------------------------------------------------------------+
//| Matrix multiply with vector |
//+------------------------------------------------------------------+
template<typename T>
vector<T> matmul(const matrix<T>& matrix_a, const vector<T>& vector_b)
{
matrix<T> matrix_b = matrix<T>::Zeros(vector_b.Size(),1);
matrix_b.Col(vector_b,0);
matrix<T>out = matmul(matrix_a,matrix_b);
return out.Col(0);
}
//+------------------------------------------------------------------+
//| Find indices where elements should be inserted to maintain order |
//+------------------------------------------------------------------+
template<typename T>
bool searchsorted(vector<T>&sorted,vector<T>&sample,bool right_side, vector &out)
{
out = vector::Zeros(sample.Size());
vector<T> samplesorted = sample;
if(!np::sort(samplesorted,sorted[0]<=sorted[sorted.Size()-1]))
{
Print(__FUNCTION__, " error ", GetLastError());
return false;
}
for(ulong j = 0; j<sample.Size(); ++j)
{
for(ulong i = 1; i<sorted.Size(); ++i)
{
if((!right_side && sorted[i-1]<sample[j] && sample[j]<=sorted[i]) ||
(right_side && sorted[i-1]<=sample[j] && sample[j]<sorted[i]))
{
out[j] = double(i);
break;
}
else
{
if((!right_side && sample[j]<=sorted.Min()) || (right_side && sample[j]<sorted.Min()))
{
out[j] = 0;
break;
}
else
{
if((!right_side && sample[j]>sorted.Max()) || (right_side && sample[j]>=sorted.Max()))
{
out[j] = double(sorted.Size());
break;
}
}
}
}
}
return true;
}
//+------------------------------------------------------------------+
//|matrix and vector addition with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> add(matrix<T>&mat, vector<T>&vec)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size() == mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(mat.Row(i)+vec,i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(mat.Col(i)+vec,i);
return out;
}
//+------------------------------------------------------------------+
//|matrix and vector subtraction with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> minus(matrix<T>&mat, vector<T>&vec)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size() == mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(mat.Row(i)-vec,i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(mat.Col(i)-vec,i);
return out;
}
//+------------------------------------------------------------------+
//|matrix and vector subtraction with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> minus(vector<T>&vec, matrix<T>&mat)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size()==mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(vec-mat.Row(i),i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(vec-mat.Col(i),i);
return out;
}
//+------------------------------------------------------------------+
//| vector subtraction with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
vector<T> minus(vector<T>& va, vector<T>& vb)
{
if(va.Size()==vb.Size())
return va-vb;
else
{
if(va.Size()!=1 && vb.Size()!=1)
{
Print(__FUNCTION__, " operands could not be broadcast together ");
return vector<T>::Zeros(0);
}
else
{
if(vb.Size()>va.Size())
return va[0] - vb;
else
return va - vb[0];
}
}
}
//+------------------------------------------------------------------+
//| vector addition with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
vector<T> add(vector<T>& va, vector<T>& vb)
{
if(va.Size()==vb.Size())
return va+vb;
else
{
if(va.Size()!=1 && vb.Size()!=1)
{
Print(__FUNCTION__, " operands could not be broadcast together ");
return vector<T>::Zeros(0);
}
else
{
if(vb.Size()>va.Size())
return va[0] + vb;
else
return va + vb[0];
}
}
}
//+------------------------------------------------------------------+
//| vector multiplication with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
vector<T> multiply(vector<T>& va, vector<T>& vb)
{
if(va.Size()==vb.Size())
return va*vb;
else
{
if(va.Size()!=1 && vb.Size()!=1)
{
Print(__FUNCTION__, " operands could not be broadcast together ");
return vector<T>::Zeros(0);
}
else
{
if(vb.Size()>va.Size())
return va[0]*vb;
else
return va*vb[0];
}
}
}
//+------------------------------------------------------------------+
//| vector division with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
vector<T> divide(vector<T>& va, vector<T>& vb)
{
if(va.Size()==vb.Size())
return va/vb;
else
{
if(va.Size()!=1 && vb.Size()!=1)
{
Print(__FUNCTION__, " operands could not be broadcast together ");
return vector<T>::Zeros(0);
}
else
{
if(vb.Size()>va.Size())
return va[0]/vb;
else
return va/vb[0];
}
}
}
//+------------------------------------------------------------------+
//|matrix and vector mulitiplication with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> multiply(matrix<T>&mat, vector<T>&vec)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size() == mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(mat.Row(i)*vec,i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(mat.Col(i)*vec,i);
return out;
}
//+------------------------------------------------------------------+
//|matrix and vector division with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> divide(matrix<T>&mat, vector<T>&vec)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size() == mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(mat.Row(i)/vec,i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(mat.Col(i)/vec,i);
return out;
}
//+------------------------------------------------------------------+
//|matrix and vector division with broadcasting |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> divide(vector<T>&vec, matrix<T>&mat)
{
if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
{
Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
return matrix<T>::Zeros(0,0);
}
matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
if(vec.Size() == mat.Cols())
for(ulong i = 0; i<mat.Rows(); ++i)
out.Row(vec/mat.Row(i),i);
else
for(ulong i = 0; i<mat.Cols(); ++i)
out.Col(vec/mat.Col(i),i);
return out;
}
//+------------------------------------------------------------------+
//| multivariate normal |
//+------------------------------------------------------------------+
matrix multivariate_normal(vector& mean, matrix& covar, ulong size)
{
if(mean.Size()!=covar.Rows() || covar.Rows()!=covar.Cols())
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(0,0);
}
matrix out = matrix::Zeros(size,mean.Size());
CMatrixDouble t = covar;
double singular_values[];
CMatrixDouble U,Vt;
matrix Ux,Vx;
vector Sx;
ulong test = 0;
if(!CAlglib::RMatrixSVD(t,t.Rows(),t.Cols(),2,2,2,singular_values,U,Vt))
{
Print(__FUNCTION__," error ", GetLastError());
return out;
}
else
{
Sx.Assign(singular_values);
Ux = U.ToMatrix();
Vx = Vt.ToMatrix();
matrix vs = np::multiply(Vx.Transpose(),Sx);
vs = vs.MatMul(Vx);
test = vs.Compare(covar,1.e-8);
}
if(test)
{
Print(__FUNCTION__, " covariance is not symmetric positive-semidefinite");
return out;
}
matrix factor = np::multiply(Vx,sqrt(Sx));
CHighQualityRandStateShell* rstate=new CHighQualityRandStateShell();
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
CAlglib::HQRndNormalM(rstate,int(size),int(covar.Cols()),out);
delete rstate;
out = out.MatMul(factor.Transpose());
out = np::add(out,mean);
return out;
}
//+------------------------------------------------------------------+
//| Flattens the array |
//+------------------------------------------------------------------+
template<typename T>
vector<T>ravelMultiIndex(vector<T> &in[], ulong&dims[])
{
vector<T> vdims;
if(!vdims.Assign(dims))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(1);
}
vector<T> vec = np::sliceVector(vdims,1);
if(!np::reverseVector(vec))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(1);
}
vec = vec.CumProd();
if(!np::reverseVector(vec))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(1);
}
vector<T> nvec = vector::Ones(vec.Size()+1);
if(!vectorCopy(nvec,vec,1))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(1);
}
matrix<T>out(in.Size(),in[0].Size());
for(ulong i = 0; i<out.Rows(); ++i)
if(!out.Row(in[i],i))
{
Print(__FUNCTION__, " error ", GetLastError());
return vector::Zeros(1);
}
return nvec.MatMul(out);
}
//+------------------------------------------------------------------+
//| get the bin edges |
//+------------------------------------------------------------------+
template<typename T>
vector<T> binEdges(vector<T>&in, ulong numbins=3)
{
T first_edge = in.Min();
T last_edge = in.Max();
return np::linspace(first_edge,last_edge,numbins);
}
/*+------------------------------------------------------------------+
//| numpy digitize for scalar values |
//+------------------------------------------------------------------+
template<typename T>
long digitize(T sample, vector<T>&edges, bool right_side=false)
{
for(long i = 1; i<long(edges.Size()); ++i)
{
if((right_side && edges[i-1]<sample && sample<=edges[i]) ||
(!right_side && edges[i-1]<=sample && sample<edges[i]))
return i;
else
{
if((right_side && sample<=edges.Min()) || (!right_side && sample<edges.Min()))
return i;
else
{
if((right_side && sample>edges.Max()) || (!right_side && sample>=edges.Max()))
return long(edges.Size());
}
}
}
return -1;
}
//+------------------------------------------------------------------+
//| Compute the multidimensional histogram of some data. |
//+------------------------------------------------------------------+
template<typename T>
bool histogramdd(matrix<T>&in, ulong bins, vector &hist, vector &edges[])
{
hist = vector::Zeros(ulong(pow(bins,in.Cols())));
if(edges.Size())
ArrayFree(edges);
if(!bins)
bins=3;
if(ArrayResize(edges,int(in.Cols()))<0)
{
Print(__FUNCTION__, " error ", GetLastError());
return false;
}
for(ulong i = 0; i<in.Cols(); ++i)
edges[i] = np::linspace(in.Col(i).Min(), in.Col(i).Max(),bins+1);
for(ulong i = 0; i<in.Rows(); ++i)
{
vector<T>row = in.Row(i);
long final_index=0;
for(ulong j =0; j<row.Size(); ++j)
{
long index = np::digitize(row[j],edges[j]) - 1;
if(ulong(index) == (edges[j].Size()-1))
index -=1;
final_index+=(index*long(pow(bins,row.Size()-(j+1))));
}
if(final_index>=long(hist.Size()) || final_index<0)
{
Print(__FUNCTION__, " hist index out of bounds. Index ", final_index, " hist size ", hist.Size());
return false;
}
hist[final_index] +=1.0;
}
return true;
}
//+------------------------------------------------------------------+
//| Compute the multidimensional histogram of some data. |
//+------------------------------------------------------------------+
template<typename T>
bool histogramdd(matrix<T>&in, ulong &bins[], vector &hist, vector &edges[])
{
hist = vector::Zeros(np::internal::prodArray(bins));
if(edges.Size())
ArrayFree(edges);
if(ulong(bins.Size())<in.Cols())
{
Print(__FUNCTION__, " invalid parameter, bins array is of insufficient size ");
return false;
}
if(bins[ArrayMinimum(bins)]==0)
{
Print(__FUNCTION__, " invalid parameter, bins array contains zero value ");
return false;
}
if(ArrayResize(edges,int(in.Cols()))<0)
{
Print(__FUNCTION__, " error ", GetLastError());
return false;
}
for(ulong i = 0; i<in.Cols(); ++i)
edges[i] = np::linspace(in.Col(i).Min(), in.Col(i).Max(),bins[i]+1);
for(ulong i = 0; i<in.Rows(); ++i)
{
vector<T>row = in.Row(i);
long final_index=0;
for(ulong j =0; j<row.Size(); ++j)
{
long index = np::digitize(row[j],edges[j]) - 1;
if(ulong(index) == (edges[j].Size()-1))
index -=1;
final_index+=(index*long(np::internal::prodArray(bins,uint(j+1))));
}
if(final_index>=long(hist.Size()) || final_index<0)
{
Print(__FUNCTION__, " hist index out of bounds. Index ", final_index, " hist size ", hist.Size());
return false;
}
hist[final_index] +=1.0;
}
return true;
}
//+------------------------------------------------------------------+
//| Compute the one dimensional histogram of some data. |
//+------------------------------------------------------------------+
template<typename T>
bool histogram(vector<T>&in, ulong bins, vector &hist, vector &edges)
{
hist = vector::Zeros(in.Size());
if(!bins)
bins=3;
if(in.Size()/2 < bins)
{
Print(__FUNCTION__, " invalid parameters, too few observations relative to the number of bins ");
return false;
}
edges = np::linspace(in.Min(), in.Max(),bins+1);
for(ulong j =0; j<in.Size(); ++j)
{
long index = np::digitize(in[j],edges) - 1;
if(ulong(index) == (edges.Size()-1))
index -=1;
hist[index] +=1.0;
}
return true;
}
//+------------------------------------------------------------------+
//|Computes empirical quantiles for a matrix. |
//+------------------------------------------------------------------+
matrix quantiles(matrix &in,vector &probs, double alphap=0.4,double betap=0.4, ulong axis = 0, double lower=-DBL_MIN, double upper=DBL_MAX)
{
matrix d = in;
if(!d.Clip(lower,upper))
{
Print(__FUNCTION__, " Clip() failure ", GetLastError());
return matrix::Zeros(probs.Size(),in.Cols());
}
vector m = alphap + probs*(1.0-alphap-betap);
matrix out(axis?in.Rows():probs.Size(),axis?probs.Size():in.Cols());
ulong step = (axis)?in.Rows():in.Cols();
for(ulong i = 0; i<step; ++i)
{
vector vec = (axis)?d.Row(i):d.Col(i);
vector outvec = np::internal::quantile(vec,m,probs);
if((axis && !out.Row(outvec,i)) || (!axis && !out.Col(outvec,i)))
{
Print(__FUNCTION__, " vector insertion failure ", GetLastError());
return matrix::Zeros(out.Rows(),out.Cols());
}
}
return out;
}*/
//+------------------------------------------------------------------+
//|BitwiseNot |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseNot(vector<T>&left)
{
int aleft[];
int res[];
if(left.Size()<1)
{
Print(__FUNCTION__, " invalid function input ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size());
}
if(!MathBitwiseNot(aleft,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseAnd |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseAnd(vector<T>&left,vector<T>&right)
{
int aleft[];
int aright[];
int res[];
if(left.Size()!=right.Size() || right.Size()<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size());
}
if(!MathBitwiseAnd(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseOr |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseOr(vector<T>&left,vector<T>&right)
{
int aleft[];
int aright[];
int res[];
if(left.Size()!=right.Size() || right.Size()<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size())
}
if(!MathBitwiseOr(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseXor |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseXor(vector<T>&left,vector<T>&right)
{
int aleft[];
int aright[];
int res[];
if(left.Size()!=right.Size() || right.Size()<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size())
}
if(!MathBitwiseXor(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseShiftL |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseShiftL(vector<T>&left,int shift)
{
int aleft[];
int res[];
if(!left.Size() || shift<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size())
}
if(!MathBitwiseShiftL(aleft,shift,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseShiftR |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseShiftR(vector<T>&left,int shift)
{
int aleft[];
int res[];
if(!left.Size() || shift<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return vector::Zeros(left.Size());
}
if(!vecAsArray(left,aleft))
{
Print(__FUNCTION__, " vector conversion failed ");
return vector::Zeros(left.Size())
}
if(!MathBitwiseShiftR(aleft,shift,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return vector::Zeros(left.Size());
}
vector<T> vec;
if(!vec.Assign(res))
{
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
return vector::Zeros(left.Size());
}
return vec;
}
//+------------------------------------------------------------------+
//|BitwiseNot |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseNot(matrix<T>&left)
{
int aleft[];
int res[];
if(!left.Rows() || !left.Cols())
{
Print(__FUNCTION__, " invalid function input ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseNot(aleft,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//|BitwiseAnd |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseAnd(matrix<T>&left,matrix<T>&right)
{
int aleft[];
int aright[];
int res[];
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseAnd(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//|BitwiseOr |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseOr(matrix<T>&left,matrix<T>&right)
{
int aleft[];
int aright[];
int res[];
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseOr(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//|BitwiseXor |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseXor(matrix<T>&left,matrix<T>&right)
{
int aleft[];
int aright[];
int res[];
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseXor(aleft,aright,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//|BitwiseShiftL |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseShiftL(matrix<T>&left,int shift)
{
int aleft[];
int res[];
if(!left.Rows() || !left.Cols() || shift<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseShiftL(aleft,shift,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//|BitwiseShiftR |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseShiftR(matrix<T>&left,int shift)
{
int aleft[];
int res[];
if(!left.Rows() || !left.Cols() || shift<1)
{
Print(__FUNCTION__, " invalid function inputs ");
return matrix::Zeros(left.Rows(),left.Cols());
}
if(!matAsArray(left,aleft))
{
Print(__FUNCTION__, " matrix conversion failed ");
return matrix::Zeros(left.Rows(),left.Cols())
}
if(!MathBitwiseShiftR(aleft,shift,res))
{
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
matrix<T> mat;
if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
{
Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
return matrix::Zeros(left.Rows(),left.Cols());
}
return mat;
}
//+------------------------------------------------------------------+
//| plot a matrix |
//+------------------------------------------------------------------+
template<typename T>
CGraphic* plotMatrix(matrix<T>&in,string plot_name,string legend="col", string x_axis_label="x_axis",string y_axis_label = "y_axis", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=550, int y2=310, bool chart_show=true, int linewidth=5, ENUM_CURVE_TYPE curve_type=CURVE_LINES)
{
CGraphic* graph = new CGraphic();
ChartRedraw(chart_id);
ChartSetInteger(chart_id, CHART_SHOW, chart_show);
if(!graph.Create(chart_id, plot_name, sub_win, x1, y1, x2, y2))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
int cols = (int)in.Cols(), rows = (int)in.Rows();
string colnames[];
int split = StringSplit(legend,StringGetCharacter(",",0),colnames);
if(split<cols)
Print(__FUNCTION__," WARNING: Invalid legend labels ");
T x[], y[];
ArrayResize(x, rows);
for(int i=0; i<rows; ++i)
x[i] = T(i+1);
CColorGenerator generator;
for(int i=0; i<cols; ++i)
{
if(!np::vecAsArray(in.Col(i), y))
{
Print(__FUNCTION__, " error ", GetLastError());
ObjectDelete(chart_id,plot_name);
delete graph;
return NULL;
}
CCurve* newcurve = graph.CurveAdd(x, y, generator.Next(),curve_type,split<cols?legend+string(i+1):colnames[i]);
newcurve.PointsFill(points_fill);
newcurve.LinesWidth(linewidth);
}
graph.XAxis().Name(x_axis_label);
graph.XAxis().NameSize(13);
graph.YAxis().Name(y_axis_label);
graph.YAxis().NameSize(13);
graph.FontSet("Lucida Console", 13);
graph.CurvePlotAll();
graph.Update();
return graph;
}
//+------------------------------------------------------------------+
//| plot a matrix |
//+------------------------------------------------------------------+
template<typename T>
CGraphic* plotMatrices(matrix<T>&x,matrix<T>&y,string plot_name, bool by_column = false, string x_axis_label="x_axis",string y_axis_label = "y_axis", string legend="col", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=550, int y2=310, bool chart_show=true)
{
if((by_column && x.Rows()!=y.Rows()) || (!by_column && x.Cols()!=y.Cols()))
{
Print(__FUNCTION__," invalid inputs ");
return NULL;
}
CGraphic* graph = new CGraphic();
ChartRedraw(chart_id);
ChartSetInteger(chart_id, CHART_SHOW, chart_show);
if(!graph.Create(chart_id, plot_name, sub_win, x1, y1, x2, y2))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
int cols = (int)y.Cols(), rows = (int)y.Rows();
string colnames[];
int split = StringSplit(legend,StringGetCharacter(",",0),colnames);
if((by_column && split<cols) || (!by_column && split<rows))
Print(__FUNCTION__," WARNING: Invalid legend labels ");
T xx[], yy[];
CColorGenerator generator;
int ncurves = (by_column)?cols:rows;
for(int i=0; i<ncurves; ++i)
{
if((by_column && (!np::vecAsArray(y.Col(i), yy) || !np::vecAsArray(x.Col(i), xx))) ||
(!by_column && (!np::vecAsArray(y.Row(i), yy) || !np::vecAsArray(x.Row(i), xx))))
{
Print(__FUNCTION__, " error ", GetLastError());
ObjectDelete(chart_id,plot_name);
delete graph;
return NULL;
}
CCurve* newcurve = graph.CurveAdd(xx, yy, generator.Next(), CURVE_POINTS_AND_LINES,((by_column && split<cols) || (!by_column && split<rows))?legend+string(i+1):colnames[i]);
newcurve.PointsFill(points_fill);
newcurve.LinesWidth(ncurves-i);
}
graph.XAxis().Name(x_axis_label);
graph.XAxis().NameSize(13);
graph.YAxis().Name(y_axis_label);
graph.YAxis().NameSize(13);
graph.FontSet("Lucida Console", 13);
graph.CurvePlotAll();
graph.Update();
return graph;
}
//+------------------------------------------------------------------+
//| plot a matrix |
//+------------------------------------------------------------------+
template<typename T>
CGraphic* plotXY(vector<T>&x,vector<T>&y,string plot_name, string x_axis_label="x_axis",string y_axis_label = "y_axis", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=750, int y2=500, bool chart_show=true,int line_width = 3, ENUM_CURVE_TYPE curve_type=CURVE_POINTS_AND_LINES)
{
T xa[];
T ya[];
if(x.Size()!=y.Size() || !np::vecAsArray(x,xa) || !np::vecAsArray(y,ya))
{
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
CGraphic* graph = new CGraphic();
ChartRedraw(chart_id);
ChartSetInteger(chart_id, CHART_SHOW, chart_show);
if(!graph.Create(chart_id, plot_name, sub_win, x1, y1, x2, y2))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
CColorGenerator generator;
CCurve* newcurve = graph.CurveAdd(xa, ya, generator.Next(),curve_type,plot_name);
newcurve.PointsFill(points_fill);
newcurve.LinesWidth(line_width);
graph.XAxis().Name(x_axis_label);
graph.XAxis().NameSize(13);
graph.YAxis().Name(y_axis_label);
graph.YAxis().NameSize(13);
graph.FontSet("Lucida Console", 13);
graph.CurvePlotAll();
graph.Update();
return graph;
}
//+------------------------------------------------------------------+
//| plot vectors |
//+------------------------------------------------------------------+
template<typename T>
void plotxy(vector<T>&x,vector<T>&y,string plot_name, string x_axis_label="x_axis",string y_axis_label = "y_axis", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=750, int y2=500, bool chart_show=true,int line_width = 3, ENUM_CURVE_TYPE curve_type=CURVE_POINTS_AND_LINES,int duration = 20)
{
CGraphic *plot = plotXY(x,y,plot_name,x_axis_label,y_axis_label,points_fill,chart_id,sub_win,x1,y1,x2,y2,chart_show,line_width,curve_type);
if(plot!=NULL)
{
Sleep(duration*1000);
plot.Destroy();
delete plot;
ChartRedraw(chart_id);
return;
}
}
//+------------------------------------------------------------------+
//| plot a scatterplot |
//+------------------------------------------------------------------+
template<typename T>
CGraphic* scatterplot(vector<T>&x,vector<T>&y,string plot_name, string x_axis_label="x_axis",string y_axis_label = "y_axis", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=550, int y2=310, bool chart_show=true)
{
T xa[];
T ya[];
if(x.Size()!=y.Size() || !np::vecAsArray(x,xa) || !np::vecAsArray(y,ya))
{
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
CGraphic* graph = new CGraphic();
ChartRedraw(chart_id);
ChartSetInteger(chart_id, CHART_SHOW, chart_show);
if(!graph.Create(chart_id, plot_name, sub_win, x1, y1, x2, y2))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
CColorGenerator generator;
CCurve* newcurve = graph.CurveAdd(xa, ya, generator.Next(), CURVE_POINTS,plot_name);
newcurve.PointsFill(points_fill);
graph.XAxis().Name(x_axis_label);
graph.XAxis().NameSize(15);
graph.XAxis().ValuesSize(15);
graph.YAxis().Name(y_axis_label);
graph.YAxis().NameSize(15);
graph.YAxis().ValuesSize(15);
graph.FontSet("Lucida Console", 18);
graph.CurvePlotAll();
graph.Update();
return graph;
}
//+------------------------------------------------------------------+
//| plot a scatterplot |
//+------------------------------------------------------------------+
template<typename T>
void scatter_plot(vector<T>&x,vector<T>&y,string plot_name, string x_axis_label="x_axis",string y_axis_label = "y_axis", bool points_fill = true,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=550, int y2=310, bool chart_show=true,int duration=20)
{
CGraphic *plot = scatterplot(x,y,plot_name,x_axis_label,y_axis_label,points_fill,chart_id,sub_win,x1,y1,x2,y2,chart_show);
if(plot!=NULL)
{
Sleep(duration*1000);
plot.Destroy();
delete plot;
ChartRedraw(chart_id);
return;
}
}
//+------------------------------------------------------------------+
//| plot a scatterplot |
//+------------------------------------------------------------------+
template<typename T>
CGraphic* histogram_plot(vector<T>&x,ulong num_bins,string plot_name, string x_axis_label="x_axis",string y_axis_label = "y_axis", int histogram_width = 6,long chart_id=0, int sub_win=0,int x1=30, int y1=40, int x2=550, int y2=310, bool chart_show=true)
{
T xa[];
T ya[];
if(x.Size()<num_bins*10)
{
Print(__FUNCTION__," invalid parameter: num_bins. Invalid relative to size of dataset");
return NULL;
}
double range=x.Max()-x.Min();
double width=range/double(num_bins);
if(width==0)
{
Print(__FUNCTION__," width is zero ");
return NULL;
}
ArrayResize(xa,(int)num_bins);
ArrayResize(ya,(int)num_bins);
//--- define the interval centers
for(ulong i=0; i<num_bins; ++i)
{
xa[i]=x.Min()+(i+0.5)*width;
ya[i]=0.0;
}
//--- fill the frequencies of falling within the interval
for(ulong i=0; i<x.Size(); ++i)
{
ulong ind=ulong((x[i]-x.Min())/width);
if(ind>=num_bins)
ind=num_bins-1;
++ya[ind];
}
CGraphic* graph = new CGraphic();
ChartRedraw(chart_id);
ChartSetInteger(chart_id, CHART_SHOW, chart_show);
if(!graph.Create(chart_id, plot_name, sub_win, x1, y1, x2, y2))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return NULL;
}
CColorGenerator generator;
CCurve* newcurve = graph.CurveAdd(xa, ya, generator.Next(), CURVE_HISTOGRAM,plot_name);
newcurve.HistogramWidth(histogram_width);
graph.XAxis().Name(x_axis_label);
graph.XAxis().NameSize(15);
graph.XAxis().ValuesSize(15);
graph.YAxis().Name(y_axis_label);
graph.YAxis().NameSize(15);
graph.YAxis().ValuesSize(15);
graph.FontSet("Lucida Console", 18);
graph.CurvePlotAll();
graph.Update();
return graph;
}
//+------------------------------------------------------------------+
//| write matrix to csv file |
//+------------------------------------------------------------------+
template <typename T>
bool matrix2csv(string csv_name, matrix<T> &matrix_, string header_string="", bool common=false, int digits=8)
{
FileDelete(csv_name);
int handle = FileOpen(csv_name,FILE_WRITE|FILE_CSV|FILE_ANSI|(common?FILE_COMMON:FILE_ANSI),",",CP_UTF8);
if(header_string == "")
for(ulong i=0; i<matrix_.Cols(); ++i)
header_string += IntegerToString(long(i)) + (i==matrix_.Cols()-1?"":",");
if(handle == INVALID_HANDLE)
{
printf("Invalid %s handle Error %d ",csv_name,GetLastError());
return (false);
}
string concstring;
vector<T> row = {};
datetime time_start = GetTickCount(), current_time;
string header[];
ushort u_sep;
u_sep = StringGetCharacter(",",0);
StringSplit(header_string,u_sep, header);
vector<T> colsinrows = matrix_.Row(0);
if(ArraySize(header) != (int)colsinrows.Size())
{
printf("headers=%d and columns=%d from the matrix vary is size ",ArraySize(header),colsinrows.Size());
return false;
}
//---
string header_str = "";
for(int i=0; i<ArraySize(header); ++i)
header_str += header[i] + (i+1 == colsinrows.Size() ? "" : ",");
FileWrite(handle,header_str);
FileSeek(handle,0,SEEK_SET);
for(ulong i=0; i<matrix_.Rows() && !IsStopped(); ++i)
{
ZeroMemory(concstring);
row = matrix_.Row(i);
for(ulong j=0, cols =1; j<row.Size() && !IsStopped(); ++j, ++cols)
{
current_time = GetTickCount();
concstring += (string)NormalizeDouble(row[j],digits) + (cols == matrix_.Cols() ? "" : ",");
}
FileSeek(handle,0,SEEK_END);
FileWrite(handle,concstring);
}
FileClose(handle);
return (true);
}
//+------------------------------------------------------------------+
//| read csv file into matrix |
//+------------------------------------------------------------------+
matrix readcsv(string file_name,bool common=false,string delimiter=",",bool handle_date_column=false,int date_column_index=-1,bool skip_header=true)
{
string Arr[],csv_header[];
int all_size = 0;
int cols_total=0;
int handle = FileOpen(file_name,FILE_SHARE_READ|FILE_CSV|FILE_ANSI|(common?FILE_COMMON:FILE_ANSI),delimiter);
if(handle == INVALID_HANDLE)
{
printf("Invalid %s handle Error %d ",file_name,GetLastError());
Print(GetLastError()==0?" TIP | File Might be in use Somewhere else or in another Directory":"");
}
else
{
int column = 0;
int rows = skip_header?0:1;
while(!FileIsEnding(handle) && !IsStopped())
{
string data = FileReadString(handle);
//---
if((rows == 0 && skip_header) || (rows == 1 && !skip_header))
{
ArrayResize(csv_header,column+1);
csv_header[column] = skip_header?data:string(column);
}
++column;
if(rows>0) //Avoid the first column which contains the column's header
{
++all_size;
ArrayResize(Arr,all_size);
Arr[all_size-1] = data;
}
//---
if(FileIsLineEnding(handle))
{
cols_total=column;
++rows;
column = 0;
}
}
FileClose(handle);
}
int rows =all_size/cols_total;
Comment("");
matrix mat(rows, cols_total);
string Col[];
vector col_vector;
for(int i=0; i<cols_total; ++i)
{
internal::col(Arr, Col, i+1, cols_total);
mat.Col(internal::formatCol(Col,handle_date_column,date_column_index,i), i);
}
return(mat);
}
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the |
//| elements of the vector. |
//| |
//| Arguments: |
//| in : vector |
//| count : Number of elements to choose |
//| result : Output vector |
//| |
//| Return value: true if successful, otherwise false. |
//+------------------------------------------------------------------+
template<typename T>
bool sampleVector(vector<T> &in,const ulong count,vector<T> &result)
{
if(!in.Size())
return false;
//--- prepare target array and calculate values
if(result.Size()<count)
if(!result.Resize(count))
return(false);
for(ulong i=0; i<count; ++i)
{
int ind=(int)(in.Size()*MathRand()/32768);
result[i]=in[ind];
}
//---
return true;
}
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the |
//| elements of the vector with replacement. |
//| |
//| Arguments: |
//| in : vector with values |
//| count : Number of elements to choose |
//| replace : If true, Sampling with Replacement will be used |
//| result : Output vector |
//| |
//| Return value: true if successful, otherwise false. |
//+------------------------------------------------------------------+
template<typename T>
bool sampleVector(vector<T> &in,const ulong count,const bool replace,vector<T> &result)
{
//--- Sampling with Replacement
if(replace)
return sampleVector(in,count,result);
if(!in.Size())
return(false);
//--- prepare target array
if(result.Size()<count)
if(!result.Resize(count))
return(false);
//--- unique values needed, prepare indices
ulong indices[];
if(!arange(indices,int(in.Size())))
return(false);
for(ulong i=0; i<count; ++i)
{
//--- select random index and swap items [j] and [i]
ulong j=i+ulong((ulong)(in.Size()-i)*MathRand()/32768);
if(j!=i)
{
ulong t=indices[i];
indices[i]=indices[j];
indices[j]=t;
}
}
//--- select data according to indices
for(ulong i=0; i<count; ++i)
result[i]=in[indices[i]];
//---
return(true);
}
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the |
//| elements of the vector. |
//| |
//| Arguments: |
//| in : vector |
//| count : Number of elements to choose |
//| result : Output vector |
//| |
//| Return value: true if successful, otherwise false. |
//+------------------------------------------------------------------+
template<typename T>
bool sampleArray(T &array[],const ulong count,T &result[])
{
if(!array.Size())
return false;
//--- prepare target array and calculate values
if(result.Size()<uint(count))
if(ArrayResize(result,int(count))!=int(count))
return(false);
for(ulong i=0; i<count; ++i)
{
int ind=(int)(array.Size()*MathRand()/32768);
result[i]=array[ind];
}
//---
return true;
}
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the |
//| elements of the vector with replacement. |
//| |
//| Arguments: |
//| array : array with values |
//| count : Number of elements to choose |
//| replace : If true, Sampling with Replacement will be used |
//| result : Output vector |
//| |
//| Return value: true if successful, otherwise false. |
//+------------------------------------------------------------------+
template<typename T>
bool sampleArray(T &array[],const ulong count,const bool replace,T &result[])
{
//--- Sampling with Replacement
if(replace)
return sampleArray(array,count,result);
if(!array.Size())
return(false);
//--- prepare target array
if(result.Size()<uint(count))
if(ArrayResize(result,int(count))!=int(count))
return(false);
//--- unique values needed, prepare indices
ulong indices[];
if(!arange(indices,int(array.Size())))
return(false);
for(ulong i=0; i<count; ++i)
{
//--- select random index and swap items [j] and [i]
ulong j=i+ulong((ulong)(array.Size()-i)*MathRand()/32768);
if(j!=i)
{
ulong t=indices[i];
indices[i]=indices[j];
indices[j]=t;
}
}
//--- select data according to indices
for(ulong i=0; i<count; ++i)
result[i]=array[indices[i]];
//---
return(true);
}
//+------------------------------------------------------------------+
//| EmpiricalProbabilityDensityEmpirical |
//+------------------------------------------------------------------+
template<typename T>
bool epdf(const vector<T> &vec,const ulong count,vector<T> &x,vector<T> &pdf)
{
if(count<=1)
return(false);
ulong size=vec.Size();
if(size==0)
return(false);
//--- check NaN values
for(ulong i=0; i<size; ++i)
{
if(!MathIsValidNumber(vec[i]))
return(false);
}
//--- prepare output arrays
if(x.Size()<count)
if(!x.Resize(count))
return(false);
if(pdf.Size()<count)
if(pdf.Resize(count))
return(false);
//--- search for min,max and range
double minv=vec.Min();
double maxv=vec.Max();
double range=maxv-minv;
if(range==0)
return(false);
//--- calculate probability density of the empirical distribution
for(ulong i=0; i<count; ++i)
{
x[i]=minv+i*range/(count-1);
pdf[i]=0;
}
for(ulong i=0; i<size; ++i)
{
double v=(vec[i]-minv)/range;
ulong ind=ulong((v*(count-1)));
++pdf[ind];
}
//--- normalize values
double dx=range/count;
double sum=0;
for(ulong i=0; i<count; ++i)
sum+=pdf[i]*dx;
if(sum==0)
return(false);
double coef=1.0/sum;
for(ulong i=0; i<count; ++i)
pdf[i]*=coef;
//---
return(true);
}
//+------------------------------------------------------------------+
//| CumulativeDistributionEmpirical |
//+------------------------------------------------------------------+
template<typename T>
bool ecdf(const vector <T> &vec,const ulong count,vector<T> &x,vector<T> &cdf)
{
if(count<1)
return(false);
ulong size=vec.Size();
if(size==0)
return(false);
//--- check NaN values
if(vec.HasNan())
{
Print(__FUNCTION__," input data contains invalid number(s)");
return false;
}
//--- prepare output arrays
if(x.Size()<count)
if(!x.Resize(count))
return(false);
if(cdf.Size()<count)
if(!cdf.Resize(count))
return(false);
//--- search for min,max and range
double minv=vec.Min();
double maxv=vec.Max();
double range=maxv-minv;
if(range==0)
return(false);
//--- calculate probability density of the empirical distribution
vector<T> pdf;
if(!pdf.Resize(count))
return(false);
for(ulong i=0; i<count; ++i)
{
x[i]=minv+i*range/(count-1);
pdf[i]=0;
}
for(ulong i=0; i<size; ++i)
{
double v=(vec[i]-minv)/range;
ulong ind=ulong((v*(count-1)));
++pdf[ind];
}
//--- normalize values
double dx=range/count;
double sum=0;
for(ulong i=0; i<count; ++i)
sum+=pdf[i]*dx;
if(sum==0)
return(false);
double coef=1.0/sum;
for(ulong i=0; i<count; ++i)
pdf[i]*=coef;
//--- calculate cumulative distribution function
sum=0.0;
for(ulong i=0; i<count; ++i)
{
sum+=pdf[i]*dx;
cdf[i]=sum;
}
//---
return(true);
}
//+------------------------------------------------------------------+
//| MathFactorial |
//+------------------------------------------------------------------+
double factorial(const uint n)
{
double ft[21]=
{
1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600,
6227020800,87178291200,1307674368000,20922789888000,355687428096000,
6402373705728000,121645100408832000,2432902008176640000
};
if(n<=20)
//--- use values from the factorials table
return(ft[n]);
else
{
//--- calculate product starting from 20th element of factorials table
double val=ft[20];
for(uint i=21; i<=n; ++i)
val*=i;
//--- return result
return(val);
}
}
//+------------------------------------------------------------------+
//| MathCorrelationSpearman's RHO |
//+------------------------------------------------------------------+
template <typename T>
double spearmanRho(vector<T> &vec1,vector<T> &vec2)
{
ulong size=vec1.Size();
if(size<1 || vec2.Size()!=size)
{
Print(__FUNCTION__, " vector lengths donot match ");
return(EMPTY_VALUE);
}
//--- calculate ranks
vector rank_x = rankVector(vec1);
vector rank_y = rankVector(vec2);
//---
return rank_x.CorrCoef(rank_y);
}
//+------------------------------------------------------------------+
//| MathCorrelationKendall |
//+------------------------------------------------------------------+
template <typename T>
double kendalTau(const vector<T> &vec1,const vector<T> &vec2)
{
double tau = double("nan");
ulong size=vec1.Size();
if(size==0 || vec2.Size()!=size)
{
Print(__FUNCTION__, " size of input vectors donot match ");
return(tau);
}
//---
long cnt1=0,cnt2=0,cnt=0;
//---
for(long i=0; i<long(size); ++i)
{
for(long j=i+1; j<long(size); ++j)
{
T delta1=vec1[i]-vec1[j];
T delta2=vec2[i]-vec2[j];
T delta=delta1*delta2;
if(delta==0)
{
if(delta1!=0)
++cnt1;
if(delta2!=0)
++cnt2;
}
else
{
++cnt1;
++cnt2;
if(delta>0.0)
++cnt;
else
cnt--;
}
}
}
//--- calculate Kendall tau
long den=cnt1*cnt2;
if(den==0)
{
Print(__FUNCTION__, " failed zero check at line 3139 ");
return tau;
}
tau=double(cnt)/MathSqrt(den);
//---
return(tau);
}
//+------------------------------------------------------------------+
//| MathRank |
//+-------------------------------------------------------------------+
template<typename T>
vector<T> rankVector(vector<T> &in_vec)
{
vector<T> rank = vector<T>::Ones(in_vec.Size());
ulong size=in_vec.Size();
if(size<1)
return rank;
if(size==1)
{
rank[0]=1;
return rank;
}
//--- prepare arrays
vector<T>values = in_vec;
ulong indices[];
//---
if(!arange(indices,int(size)))
return rank;
//---
ulong i,j,k,t,tmpi;
double tmp;
//--- sort
if(size!=1)
{
i=2;
do
{
t=i;
while(t!=1)
{
k=t/2;
if(values[k-1]>=values[t-1])
t=1;
else
{
//--- swap
tmp=values[k-1];
values[k-1]=values[t-1];
values[t-1]=tmp;
tmpi=indices[k-1];
indices[k-1]=indices[t-1];
indices[t-1]=tmpi;
t=k;
}
}
i=i+1;
}
while(i<=size);
i=size-1;
do
{
//--- swap
tmp=values[i];
values[i]=values[0];
values[0]=tmp;
tmpi=indices[i];
indices[i]=indices[0];
indices[0]=tmpi;
t=1;
while(t!=0)
{
k=2*t;
if(k>i)
t=0;
else
{
if(k<i)
if(values[k]>values[k-1])
++k;
if(values[t-1]>=values[k-1])
t=0;
else
{
//--- swap
tmp=values[k-1];
values[k-1]=values[t-1];
values[t-1]=tmp;
tmpi=indices[k-1];
indices[k-1]=indices[t-1];
indices[t-1]=tmpi;
t=k;
}
}
}
i=i-1;
}
while(i>=1);
}
//--- compute tied ranks
i=0;
while(i<size)
{
j=i+1;
while(j<size)
{
//--- check
if(values[j]!=values[i])
break;
j=j+1;
}
for(k=i; k<j; ++k)
values[k]=1+(i+j-1)*0.5;
i=j;
}
//--- set output values
for(i=0; i<size; ++i)
rank[indices[i]]=values[i];
//---
return rank;
}
//+------------------------------------------------------------------+
//| masking vector based on < condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsLt(vector<T>&in, const T value)
{
vector<T>out = vector<T>::Zeros(in.Size());
for(ulong i = 0; i<in.Size(); ++i)
{
if(in[i]<value)
out[i] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking vector based on <= condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsLte(vector<T>&in, const T value)
{
vector<T>out = vector<T>::Zeros(in.Size());
for(ulong i = 0; i<in.Size(); ++i)
{
if(in[i]<=value)
out[i] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking vector based on > condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsGt(vector<T>&in, const T value)
{
vector<T>out = vector<T>::Zeros(in.Size());
for(ulong i = 0; i<in.Size(); ++i)
{
if(in[i]>value)
out[i] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking vector based on >= condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsGte(vector<T>&in, const T value)
{
vector<T>out = vector<T>::Zeros(in.Size());
for(ulong i = 0; i<in.Size(); ++i)
{
if(in[i]>=value)
out[i] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking matrix based on < condition |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsLt(matrix<T>&in, const T value)
{
matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(in[i][j]<value)
out[i][j] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking matrix based on <= condition |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsLte(matrix<T>&in, const T value)
{
matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(in[i][j]<=value)
out[i][j] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking matrix based on > condition |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsGt(matrix<T>&in, const T value)
{
matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(in[i][j]>value)
out[i][j] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking matrix based on >= condition |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsGte(matrix<T>&in, const T value)
{
matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(in[i][j]>=value)
out[i][j] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking vector based on == condition by digits |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsEq(vector<T>&in, const T value, int _digits = 8)
{
vector<T>out = vector<T>::Zeros(in.Size());
for(ulong i = 0; i<in.Size(); ++i)
{
if(NormalizeDouble(in[i],_digits)==value)
out[i] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| masking matrix based on == condition by digits |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsEq(matrix<T>&in, const T value, int _digits = 8)
{
matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(NormalizeDouble(in[i][j],_digits)>=value)
out[i][j] = 1.0;
}
return out;
}
//+------------------------------------------------------------------+
//| replace vector members that are not numbers |
//+------------------------------------------------------------------+
template<typename T>
vector<T> replace_whereVectorIsNaN(vector<T>&in,const T fill = 0.0)
{
vector<T>out = in;
if(in.HasNan()==in.Size())
{
out.Fill(fill);
return out;
}
for(ulong i = 0; i<in.Size(); ++i)
{
if(MathClassify(in[i])==FP_NAN)
out[i] = fill;
}
return out;
}
//+------------------------------------------------------------------+
//| replace matrix members that are not numbers |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> replace_whereMatrixIsNaN(matrix<T>&in,const T fill = 0.0)
{
matrix<T>out = in;
if(in.HasNan()==(in.Rows()*in.Cols()))
{
out.Fill(fill);
return out;
}
for(ulong i = 0; i<in.Rows(); ++i)
{
for(ulong j = 0; j<in.Cols(); ++j)
if(MathClassify(in[i][j])==FP_NAN)
out[i][j] = fill;
}
return out;
}
//+------------------------------------------------------------------+
//| Copy to select columns |
//+------------------------------------------------------------------+
template<typename T>
bool copyToSelectCols(matrix<T> &copyto, matrix<T> &from, vector &v_mask)
{
if(from.Cols()!=copyto.Cols() || v_mask.Size()!=copyto.Cols())
{
Print(__FUNCTION__, " invalid parameters ");
return false;
}
for(ulong i = 0; i<copyto.Cols(); ++i)
{
if(v_mask[i])
{
if(!copyto.Col(from.Col(i),i))
{
Print(__FUNCTION__, " failed column assignment ", GetLastError());
return false;
}
}
}
return true;
}
//+------------------------------------------------------------------+
//| Copy to select rows |
//+------------------------------------------------------------------+
template<typename T>
bool copyToSelectRows(matrix<T> &copyto, matrix<T> from, vector &v_mask)
{
if(from.Rows()!=copyto.Rows() || v_mask.Size()!=copyto.Rows())
{
Print(__FUNCTION__, " invalid parameters ");
return false;
}
for(ulong i = 0; i<copyto.Rows(); ++i)
{
if(v_mask[i])
{
if(!copyto.Row(from.Row(i),i))
{
Print(__FUNCTION__, " failed column assignment ", GetLastError());
return false;
}
}
}
return true;
}
//+------------------------------------------------------------------+
//|Copy select elements |
//+------------------------------------------------------------------+
template<typename T>
bool copySelectElements(vector<T> &copyto, vector<T> &from, vector &v_mask)
{
if(from.Size()!=copyto.Size() || v_mask.Size()!=copyto.Size())
{
Print(__FUNCTION__, " invalid parameters ");
return false;
}
for(ulong i = 0; i<copyto.Size(); ++i)
{
if(v_mask[i])
copyto[i] = from[i];
}
return true;
}
//+------------------------------------------------------------------+
//| sign |
//+------------------------------------------------------------------+
template<typename T>
T sign(const T a)
{
if(a<T(0))
return T(-1);
else
if(a>T(0))
return T(1);
else
return T(0);
}
//+------------------------------------------------------------------+
//| calculation of raw moment about specified center |
//+------------------------------------------------------------------+
double moment(vector& a, ulong order, double mean=NULL)
{
if(!a.Size())
return a.Mean();
if(order == 0 || (order == 1 && !mean))
return order==0?1.0:0.;
ulong list[];
ArrayResize(list,list.Size()+1,100);
list[list.Size()-1]=order;
ulong current = order;
while(current>2)
{
if(MathMod(current,2))
current=(current-1)/2;
else
current/=2;
ArrayResize(list,list.Size()+1,100);
list[list.Size()-1]= current;
}
double mn = (mean == 0.0)?a.Mean():mean;
vector azeromean = a - mn;
vector s;
if(list[list.Size()-1] == 1)
s = azeromean;
else
s = pow(azeromean,2.);
for(int i = int(list.Size())-2; i>=0; --i)
{
s = pow(s,2.0);
if(MathMod(i,2))
s*=azeromean;
}
return s.Mean();
}
//+------------------------------------------------------------------+
//| sin(pi*x) |
//+------------------------------------------------------------------+
template <typename T>
T sinpi(T x)
{
T s = 1.0;
if(x < 0.0)
{
x = -x;
s = -1.0;
}
T r = fmod(x, 2.0);
if(r < 0.5)
{
return s * sin(M_PI * r);
}
else
if(r > 1.5)
{
return s * sin(M_PI * (r - 2.0));
}
else
{
return -s * sin(M_PI * (r - 1.0));
}
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
double lgam_large_x(double x)
{
const double LS2PI = 0.91893853320467274178;
double q = (x - 0.5) * log(x) - x + LS2PI;
if(x > 1.0e8)
{
return (q);
}
double p = 1.0 / (x * x);
p = ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x;
return q + p;
}
//+------------------------------------------------------------------+
//| helper |
//+------------------------------------------------------------------+
double polevl(double x, double &coef[], int N)
{
double ans;
int i;
int p;
p = 0;
ans = coef[p++];
i = N;
do
{
ans = ans * x + coef[p++];
}
while(--i != 0);
return (ans);
}
//+------------------------------------------------------------------+
//| helper |
//+------------------------------------------------------------------+
double p1evl(double x, double &coef[], int N)
{
double ans;
uint p;
int i;
p = 0;
ans = x + coef[p++];
i = N - 1;
do
ans = ans * x + coef[p++];
while(--i != 0);
return (ans);
}
//+------------------------------------------------------------------+
//| sign of log gamma |
//+------------------------------------------------------------------+
double lgam_sgn(double x, int &sign)
{
double p, q, u, w, z;
int i;
const double LS2PI = 0.91893853320467274178;
const double MAXLGM = 2.556348e305;
double gamma_A[] =
{
8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
-2.77777777730099687205E-3, 8.33333333333331927722E-2
};
double gamma_B[] =
{
-1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
-1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5
};
double gamma_C[] =
{
-3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
-1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6
};
sign = 1;
if(MathClassify(x)!=FP_NORMAL)
{
return x;
}
if(x < -34.0)
{
q = -x;
w = lgam_sgn(q, sign);
p = floor(q);
if(p == q)
{
Print(__FUNCTION__, " singular error ");
return (double("inf"));
}
i = (int)p;
if((i & 1) == 0)
{
sign = -1;
}
else
{
sign = 1;
}
z = q - p;
if(z > 0.5)
{
p += 1.0;
z = p - q;
}
z = q * sinpi(z);
if(z == 0.0)
{
Print(__FUNCTION__, " singular error ");
return (double("inf"));
}
z = log(M_PI) - log(z) - w;
return (z);
}
if(x < 13.0)
{
z = 1.0;
p = 0.0;
u = x;
while(u >= 3.0)
{
p -= 1.0;
u = x + p;
z *= u;
}
while(u < 2.0)
{
if(u == 0.0)
{
Print(__FUNCTION__, " singular error ");
return (double("inf"));
}
z /= u;
p += 1.0;
u = x + p;
}
if(z < 0.0)
{
sign = -1;
z = -z;
}
else
{
sign = 1;
}
if(u == 2.0)
{
return (log(z));
}
p -= 2.0;
x = x + p;
p = x * polevl(x, gamma_B, 5) / p1evl(x, gamma_C, 6);
return (log(z) + p);
}
if(x>MAXLGM)
return sign*double("inf");
if(x>=1000.)
return lgam_large_x(x);
q = (x-0.5)*log(x) - x +LS2PI;
p = 1./(x*x);
return q+polevl(p,gamma_A,4)/x;
}
//+-------------------------------------------------------------------------+
//| Asymptotic expansion for ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1)|
//+-------------------------------------------------------------------------+
double lbeta_asymp(double a, double b, int &sgn)
{
double r = lgam_sgn(b, sgn);
r -= b * log(a);
r += b * (1 - b) / (2 * a);
r += b * (1 - b) * (1 - 2 * b) / (12 * a * a);
r += -b * b * (1 - b) * (1 - b) / (12 * a * a * a);
return r;
}
//+------------------------------------------------------------------+
//| Special case for a negative integer argument |
//+------------------------------------------------------------------+
double beta_negint(int a, double b)
{
int sgn;
if(b == int(b) && 1 - a - b > 0)
{
sgn = (MathMod(b, 2) == 0) ? 1 : -1;
return sgn * CAlglib::Beta(1. - a - b, b);
}
else
{
Print(__FUNCTION__, " overflow error ");
return double("inf");
}
}
//+------------------------------------------------------------------+
//| helper |
//+------------------------------------------------------------------+
double lbeta_negint(int a, double b)
{
double r;
if(b == int(b) && 1 - a - b > 0)
{
r = lbeta(1 - a - b, b);
return r;
}
else
{
Print(__FUNCTION__, " overflow error ");
return double("inf");
}
}
//+------------------------------------------------------------------+
//| helper |
//+------------------------------------------------------------------+
double lbeta(double a, double b)
{
double y;
int sign;
const double MAXGAM = 171.624376956302725;
const double beta_ASYMP_FACTOR = 1e6;
sign = 1;
if(a <= 0.0)
{
if(a == floor(a))
{
if(a == int(a))
{
return lbeta_negint(int(a), b);
}
else
{
Print(__FUNCTION__, " overflow error ");
return (sign * double("inf"));
}
}
}
if(b <= 0.0)
{
if(b == floor(b))
{
if(b == int(b))
{
return lbeta_negint(int(b), a);
}
else
{
Print(__FUNCTION__, " overflow error ");
return (sign * double("inf"));
}
}
}
if(MathAbs(a) < MathAbs(b))
{
y = a;
a = b;
b = y;
}
if(MathAbs(a) > beta_ASYMP_FACTOR * MathAbs(b) && a > beta_ASYMP_FACTOR)
{
/* Avoid loss of precision in lgam(a + b) - lgam(a) */
y = lbeta_asymp(a, b, sign);
return y;
}
y = a + b;
if(MathAbs(y) > MAXGAM || MathAbs(a) > MAXGAM || MathAbs(b) > MAXGAM)
{
int sgngam;
y = lgam_sgn(y, sgngam);
sign *= sgngam; /* keep track of the sign */
y = lgam_sgn(b, sgngam) - y;
sign *= sgngam;
y = lgam_sgn(a, sgngam) + y;
sign *= sgngam;
return (y);
}
y = CAlglib::GammaFunction(y);
a = CAlglib::GammaFunction(a);
b = CAlglib::GammaFunction(b);
if(y == 0.0)
{
Print(__FUNCTION__, " overflow error ");
return (sign * double("inf"));
}
if(MathAbs(MathAbs(a) - MathAbs(y)) > MathAbs(MathAbs(b) - MathAbs(y)))
{
y = b / y;
y *= a;
}
else
{
y = a / y;
y *= b;
}
if(y < 0)
{
y = -y;
}
return (log(y));
}
//+------------------------------------------------------------------+
//| Binomial coefficient |
//+------------------------------------------------------------------+
double binom(double n, double k)
{
double kx,nx,num, den, dk, sgn;
if(n < 0)
{
nx = floor(n);
if(n == nx)
{
// Undefined
return double("nan");
}
}
kx = floor(k);
if(k == kx && (MathAbs(n) > 1e-8 || n == 0))
{
nx = floor(n);
if(nx == n && kx > nx / 2 && nx > 0)
{
// Reduce kx by symmetry
kx = nx - kx;
}
if(kx >= 0 && kx < 20)
{
num = 1.0;
den = 1.0;
for(int i = 1; i < 1 + int(kx); i++)
{
num *= i + n - kx;
den *= i;
if(MathAbs(num) > 1e50)
{
num /= den;
den = 1.0;
}
}
return num / den;
}
}
// general case
if(n >= 1e10 * k && k > 0)
{
// avoid under/overflows intermediate results
return exp(-lbeta(1 + n - k, 1 + k) - log(n + 1));
}
if(k > 1e8 * MathAbs(n))
{
// avoid loss of precision
num = CAlglib::GammaFunction(1 + n) / MathAbs(k) + CAlglib::GammaFunction(1 + n) * n / (2 * k * k); // + ...
num /= M_PI * pow(MathAbs(k), n);
if(k > 0)
{
kx = floor(k);
if(int(kx) == kx)
{
dk = k - kx;
sgn = (MathMod(int(kx), 2) == 0) ? 1 : -1;
}
else
{
dk = k;
sgn = 1;
}
return num * sin((dk - n) * M_PI) * sgn;
}
kx = floor(k);
if(int(kx) == kx)
{
return 0;
}
return num * sin(k * M_PI);
}
return 1 / (n + 1) / CAlglib::Beta(1 + n - k, 1 + k);
}
//+------------------------------------------------------------------+
//| N choose K |
//+------------------------------------------------------------------+
template<typename T>
T comb(ulong N, T k, bool exact=false, bool repetition = false)
{
if(repetition)
return comb(N+ulong(k)-1,(T)k,exact);
if(exact)
return comb(N,k);
else
{
bool cond = false;
cond = ((k<=T(N)) && (N>=0) && (k>=0));
T val = (T)binom(double(N),double(k));
if(!cond)
val = 0;
return val;
}
}
//+------------------------------------------------------------------+
//| Compute the kurtosis (Fisher or Pearson) of a dataset. |
//+------------------------------------------------------------------+
double kurtosis(vector& in,bool fisher=true, bool bias=true)
{
ulong n = in.Size();
double out;
double mean = in.Mean();
double m2 = moment(in,2,mean);
double m4 = moment(in,4,mean);
bool zero = (m2<= pow(mean*2.e-16,2.));
out = (!zero)?m4/pow(m2,2.):double("nan");
if(!bias)
{
bool cancorrect = (~zero&(n>3));
double nval = 1.0/(n-2)/(n-3) * ((pow(n,2)-1.0)*m4/pow(m2,2.0) - 3*pow((n-1),2.0));
if(cancorrect)
out = nval + 3.;
else
out = double("nan");
}
if(fisher)
out = out - 3.;
return out;
}
//+------------------------------------------------------------------+
//| get element wise maximum |
//+------------------------------------------------------------------+
vector maximum(vector &x, vector &y)
{
if(x.Size()!=y.Size())
{
Print(__FUNCTION__, " vector size mismatch ");
return vector::Zeros(0);
}
matrix xy(x.Size(),2);
if(!xy.Col(x,0) || !xy.Col(y,1))
{
Print(__FUNCTION__, " column insertion error ", GetLastError());
return vector::Zeros(0);
}
return xy.Max(1);
}
//+------------------------------------------------------------------+
//| get element wise minimum |
//+------------------------------------------------------------------+
vector minimum(vector &x, vector &y)
{
if(x.Size()!=y.Size())
{
Print(__FUNCTION__, " vector size mismatch ");
return vector::Zeros(0);
}
matrix xy(x.Size(),2);
if(!xy.Col(x,0) || !xy.Col(y,1))
{
Print(__FUNCTION__, " column insertion error ", GetLastError());
return vector::Zeros(0);
}
return xy.Min(1);
}
//+------------------------------------------------------------------+
//| get element wise maximum |
//+------------------------------------------------------------------+
vector maximum(vector &x, double v)
{
vector y(x.Size());
y.Fill(v);
matrix xy(x.Size(),2);
if(!xy.Col(x,0) || !xy.Col(y,1))
{
Print(__FUNCTION__, " column insertion error ", GetLastError());
return vector::Zeros(0);
}
return xy.Max(1);
}
//+------------------------------------------------------------------+
//| get element wise minimum |
//+------------------------------------------------------------------+
vector minimum(vector &x, double v)
{
vector y(x.Size());
y.Fill(v);
matrix xy(x.Size(),2);
if(!xy.Col(x,0) || !xy.Col(y,1))
{
Print(__FUNCTION__, " column insertion error ", GetLastError());
return vector::Zeros(0);
}
return xy.Min(1);
}
//+------------------------------------------------------------------+
//| return vector with elements that meet certain condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereIsNan(vector<T>& in)
{
ulong any = in.HasNan();
if(any)
{
vector out = vector::Zeros(any)
ulong count = 0;
for(ulong i = 0; i<in.Size(); ++i)
if(MathClassify(in[i]) == FP_NAN)
out[count++] = double(i);
return out;
}
else
return vector::Zeros(0);
}
//+------------------------------------------------------------------+
//| return vector with elements that meet certain condition |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereIsNotNan(vector<T>& in)
{
ulong any = in.HasNan();
if(any)
{
vector out = vector::Zeros(in.Size() - any);
ulong count = 0;
for(ulong i = 0; i<in.Size(); ++i)
if(MathClassify(in[i]) != FP_NAN)
out[count++] = in[i];
return out;
}
else
return in;
}
//+------------------------------------------------------------------+
//| matrix rank |
//+------------------------------------------------------------------+
ulong matrix_rank(matrix& in)
{
matrix A = in;
matrix u,v;
vector s;
if(!A.SVD(u,v,s))
{
Print(__FUNCTION__, " svd failure ", GetLastError());
return 0;
}
double rtol = MathMax(A.Rows(),A.Cols())*DBL_EPSILON;
double tol = s.Max()*rtol;
vector count = whereVectorIsGt(s,tol);
return ulong(count.Sum());
}
//+------------------------------------------------------------------+
//| toeplitz matrix construction |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> toeplitz(vector<T>& x)
{
ulong n = x.Size();
matrix<T> out(n,n);
for(ulong i = 0; i<n; ++i)
for(ulong j = 0; j<n; ++j)
if(i<=j)
out[i,j] = x[j-i];
else
out[i,j] = x[i-j];
return out;
}
//+------------------------------------------------------------------+
//| Compute the (Moore-Penrose) pseudo-inverse of a matrix.
//+------------------------------------------------------------------+
matrix pinv(matrix& a)
{
double rcond = 1.e-15;
vector vrcond(a.Cols());
vrcond.Fill(rcond);
if(!a.Rows())
return a;
matrix aa = a.Conjugate();
matrix u,v;
vector s;
if(!aa.SVD(u,v,s))
{
Print(__FUNCTION__," svd failed ",GetLastError());
return matrix::Zeros(0,0);
}
CMatrixDouble m = aa;
double singular_values[];
CMatrixDouble U,Vt;
if(!CAlglib::RMatrixSVD(m,m.Rows(),m.Cols(),2,2,2,singular_values,U,Vt))
{
Print(__FUNCTION__," error ", GetLastError());
return matrix::Zeros(0,0);;
}
else
{
s.Assign(singular_values);
u = U.ToMatrix();
u = sliceMatrixCols(u,0,long(s.Size()));
v = Vt.ToMatrix();
}
vector cuttoff = vrcond*s.Max();
for(ulong i = 0; i<s.Size(); ++i)
if(s[i]>cuttoff[i])
s[i] = 1./s[i];
else
s[i] = 0.0;
matrix mats = matrix::Zeros(s.Size(),1);
mats.Col(s,0);
matrix vm = multiply(u.Transpose(),s);
return v.Transpose().MatMul(vm);
}
}
//+------------------------------------------------------------------+
Reference in a new issue View git blame Copy permalink