5050 行
313 KiB
MQL5
5050 行
313 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| np.mqh |
|
|
//| Copyright 2024, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2024, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include<Math\Stat\Normal.mqh>
|
|
#include<Math\Alglib\specialfunctions.mqh>
|
|
#include<Math\Alglib\linalg.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> ©to,matrix ©from, 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> ©to, matrix ©from, 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 vector of arbitrary elements specified as collection |
|
|
//| of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
vector<T> select(vector<T> &in,vector<T> &indices)
|
|
{
|
|
ulong size = (ulong)MathMin(in.Size(),indices.Size());
|
|
|
|
if(size<1 || indices.ArgMax()>=in.Size())
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameter ");
|
|
return vector<T>::Zeros(1);
|
|
}
|
|
|
|
vector<T> out(size);
|
|
|
|
for(ulong i=0; i<size; ++i)
|
|
out[i] = in[ulong(indices[i])];
|
|
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return vector of arbitrary elements specified as collection |
|
|
//| of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
vector<T> select(vector<T> &in,const ulong &indices[])
|
|
{
|
|
ulong size = (ulong)MathMin(in.Size(),indices.Size());
|
|
|
|
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=ulong(in.Size()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameter ");
|
|
return vector<T>::Zeros(1);
|
|
}
|
|
|
|
vector<T> out(size);
|
|
|
|
for(ulong i=0; i<size; ++i)
|
|
out[i] = in[indices[i]];
|
|
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return vector of arbitrary elements specified as collection |
|
|
//| of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
vector<T> select(vector<T> &in,const int &indices[])
|
|
{
|
|
int size = (int)MathMin(in.Size(),indices.Size());
|
|
|
|
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=int(in.Size()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameter ");
|
|
return vector<T>::Zeros(1);
|
|
}
|
|
|
|
vector<T> out(size);
|
|
|
|
for(int i=0; i<size; ++i)
|
|
out[i] = in[indices[i]];
|
|
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return vector of arbitrary elements specified as collection |
|
|
//| of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
vector<T> select(vector<T> &in,const uint &indices[])
|
|
{
|
|
uint size = (uint)MathMin(in.Size(),indices.Size());
|
|
|
|
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=uint(in.Size()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameter ");
|
|
return vector<T>::Zeros(1);
|
|
}
|
|
|
|
vector<T> out(size);
|
|
|
|
for(uint 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;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return arbitrary columns of a matrix specified as |
|
|
//|collection of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
matrix<T> selectMatrixCols(matrix<T> &in,const ulong &indices[])
|
|
{
|
|
ulong size = (ulong)MathMin(in.Cols(),indices.Size());
|
|
|
|
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=ulong(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 ulong &indices[])
|
|
{
|
|
ulong size = (ulong)MathMin(in.Rows(),indices.Size());
|
|
|
|
if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=ulong(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;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return arbitrary columns of a matrix specified as |
|
|
//|collection of indices |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
matrix<T> selectMatrixCols(matrix<T> &in,vector& indices)
|
|
{
|
|
ulong size = (ulong)MathMin(in.Cols(),indices.Size());
|
|
|
|
if(size<1 || indices.Max()>=double(in.Cols()) || indices.Min()<0.0)
|
|
{
|
|
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(ulong(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,vector& indices)
|
|
{
|
|
ulong size = (ulong)MathMin(in.Rows(),indices.Size());
|
|
|
|
if(size<1 || indices.Max()>=double(in.Rows()) || indices.Min()<0.0)
|
|
{
|
|
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(ulong(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> ©to,vector<T> ©from,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> ©to, 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> ©to,matrix ©from, 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> ©to,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, double step_value=1)
|
|
{
|
|
vector out(size);
|
|
|
|
for(ulong i = 0; i<size; ++i, value+=step_value)
|
|
out[i] = value;
|
|
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| returns vector |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> arange(T start_value, T stop_value, T step_value)
|
|
{
|
|
T numerator = (stop_value - start_value);
|
|
ulong size = ulong((numerator)/step_value);
|
|
if(MathMod(numerator,step_value)>T(0))
|
|
size+=1;
|
|
|
|
vector<T> out = vector<T>::Zeros(size);
|
|
|
|
for(ulong i = 0; i<size; ++i, start_value+=step_value)
|
|
out[i] = start_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);
|
|
|
|
if(!num)
|
|
return vector<T>::Zeros(num);
|
|
if(num==1)
|
|
{
|
|
out[0] = start;
|
|
return out;
|
|
}
|
|
|
|
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 ©_from[], matrix ©_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>
|
|
bool vecAsArray(const vector<T> &in, double &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] = double(in[i]);
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic vector to array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool vecAsArray(const vector<T> &in, int &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] = int(in[i]);
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic vector to array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool vecAsArray(const vector<T> &in, uint &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] = uint(in[i]);
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic vector to array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool vecAsArray(const vector<T> &in, long &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] = long(in[i]);
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic vector to array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool vecAsArray(const vector<T> &in, ulong &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] = ulong(in[i]);
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic matrix to 1-D array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, double &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]=double(in.Flat(i));
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic matrix to 1-D array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, int &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]=int(in.Flat(i));
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic matrix to 1-D array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, uint &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]=uint(in.Flat(i));
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic matrix to 1-D array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, long &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]=long(in.Flat(i));
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|generic matrix to 1-D array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, ulong &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]=ulong(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));
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| The function extracts the unique columns or rows from the matrix |
|
|
//| |
|
|
//| Arguments: |
|
|
//| in : matrix with values |
|
|
//| axis : int by row = 1, column = 0 |
|
|
//| |
|
|
//| Return value: out matrix of unique rows or columns |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> unique(matrix<T> &in,int axis = 0)
|
|
{
|
|
vector sum = in.Sum(axis);
|
|
vector spec;
|
|
long indices[],counts[];
|
|
if(!unique(sum,spec,indices,counts))
|
|
{
|
|
Print(__FUNCTION__,": error ");
|
|
return matrix<T>::Zeros(0,0);
|
|
}
|
|
if(axis)
|
|
return selectMatrixRows(in,indices);
|
|
else
|
|
return selectMatrixCols(in,indices);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| 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);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| read csv file into matrix |
|
|
//+------------------------------------------------------------------+
|
|
matrix readcsv_from_string(string csv_string_data,bool common=false,string delimiter=",",bool handle_date_column=false,int date_column_index=-1,bool skip_header=true)
|
|
{
|
|
string Arr[],csv_header[],lines[];
|
|
int all_size = 0;
|
|
int rows = 0;
|
|
int cols_total=0;
|
|
|
|
int handle = StringSplit(csv_string_data,StringGetCharacter("\n",0),lines);
|
|
|
|
if(handle <= 0)
|
|
{
|
|
Print(__FUNCTION__, " no lines found in the string data ");
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
else
|
|
{
|
|
int column = 0;
|
|
int shift = 0;
|
|
cols_total = StringSplit(lines[shift],StringGetCharacter(delimiter,0),csv_header);
|
|
if(StringLen(csv_header[cols_total-1]) == 0)
|
|
--cols_total;
|
|
all_size = handle - 1;
|
|
if(StringLen(lines[handle-1]) == 0)
|
|
--all_size;
|
|
rows = all_size;
|
|
all_size*=cols_total;
|
|
ArrayResize(Arr,all_size);
|
|
for(int i = 0; i<(all_size); i+=cols_total)
|
|
{
|
|
StringSplit(lines[skip_header?((i+cols_total)/cols_total):(i/cols_total)],StringGetCharacter(delimiter,0),csv_header);
|
|
ArrayCopy(Arr,csv_header,i,0,cols_total);
|
|
}
|
|
}
|
|
|
|
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;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Generates a random sample from a given 1-D array a |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> choice(vector<T>& a, vector& p, ulong size = 0, bool replace = true)
|
|
{
|
|
ulong pop_size = a.Size();
|
|
vector idx;
|
|
if(!pop_size && size)
|
|
{
|
|
Print(__FUNCTION__, " a cannot be empty unless no samples are taken ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
if(p.Size())
|
|
{
|
|
ulong d = p.Size();
|
|
double atol = DBL_EPSILON;
|
|
if(d!=pop_size)
|
|
{
|
|
Print(__FUNCTION__, " a and p must have same size ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
double psum = p.Sum();
|
|
if(p.Min()<0.0)
|
|
{
|
|
Print(__FUNCTION__, " probabilities cannot be negative ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
if(psum-1.>atol)
|
|
{
|
|
Print(__FUNCTION__, " probabilities must sum to 1 ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
}
|
|
CHighQualityRandStateShell rstate;
|
|
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
|
|
bool is_scalar = (size==0);
|
|
ulong shape = 0;
|
|
vector cdf;
|
|
if(!is_scalar)
|
|
shape = size;
|
|
else
|
|
size = 1;
|
|
if(replace)
|
|
{
|
|
if(p.Size())
|
|
{
|
|
cdf = p.CumSum();
|
|
cdf/=cdf[cdf.Size()-1];
|
|
vector uniform_samples = vector::Zeros(shape);
|
|
for(ulong i = 0; i<shape; uniform_samples[i] = CHighQualityRand::HQRndUniformR(rstate.GetInnerObj()), ++i);
|
|
if(!searchsorted(cdf,uniform_samples,true,idx))
|
|
{
|
|
Print(__FUNCTION__, " searchsorted failed ", __LINE__);
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
idx = vector::Zeros(size);
|
|
for(ulong i = 0; i<shape; idx[i] = (double)CHighQualityRand::HQRndUniformI(rstate.GetInnerObj(),int(pop_size)), ++i);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(size>pop_size)
|
|
{
|
|
Print(__FUNCTION__, " cannot take a larger sample thatn population when replace = false ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
if(p.Size())
|
|
{
|
|
ulong n_uniq = 0;
|
|
vector pp = p;
|
|
vector found = vector::Zeros(shape);
|
|
vector cdf,_new,temp;
|
|
long uindices[],counted[];
|
|
while(n_uniq<size)
|
|
{
|
|
vector x = vector::Zeros(size - n_uniq);
|
|
for(ulong i = 0; i<x.Size(); x[i] = CHighQualityRand::HQRndUniformR(rstate.GetInnerObj()), ++i);
|
|
if(n_uniq>0)
|
|
{
|
|
for(ulong i = 0; i<n_uniq; ++i)
|
|
pp[ulong(found[i])] = 0.0;
|
|
}
|
|
cdf = pp.CumSum();
|
|
cdf /= cdf[cdf.Size()-1];
|
|
if(!searchsorted(cdf,x,true,_new))
|
|
{
|
|
Print(__FUNCTION__, " searchsorted failed ", __LINE__);
|
|
return vector::Zeros(0);
|
|
}
|
|
if(!unique(_new,temp,uindices,counted) || !ArraySort(uindices))
|
|
{
|
|
Print(__FUNCTION__, " unique failed ", __LINE__);
|
|
return vector::Zeros(0);
|
|
}
|
|
_new = select(_new,uindices);
|
|
for(ulong i = n_uniq, j = 0; i<(n_uniq+_new.Size()); ++i)
|
|
found[i] = _new[j++];
|
|
n_uniq+=_new.Size();
|
|
}
|
|
idx = found;
|
|
}
|
|
else
|
|
{
|
|
idx = permutation(pop_size);
|
|
idx = sliceVector(idx,0,long(size));
|
|
idx.Resize(shape);
|
|
}
|
|
}
|
|
return select(a,idx);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| 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_NORMAL)
|
|
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> ©to, 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> ©to, 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> ©to, 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 * CBetaF::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 = CGammaFunc::GammaFunc(y);
|
|
a = CGammaFunc::GammaFunc(a);
|
|
b = CGammaFunc::GammaFunc(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 = CGammaFunc::GammaFunc(1 + n) / MathAbs(k) + CGammaFunc::GammaFunc(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) / CBetaF::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(!CSingValueDecompose::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);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| 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)(CHighQualityRand::HQRndUniformR(rstate.GetInnerObj())*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)(CHighQualityRand::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;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| folds structure |
|
|
//+------------------------------------------------------------------+
|
|
|
|
struct Folds
|
|
{
|
|
vector test_indices;
|
|
vector train_indices;
|
|
|
|
Folds(void)
|
|
{
|
|
test_indices = train_indices = vector::Zeros(0);
|
|
}
|
|
|
|
Folds(Folds& other)
|
|
{
|
|
test_indices = other.test_indices;
|
|
train_indices = other.train_indices;
|
|
}
|
|
void operator=(Folds& other)
|
|
{
|
|
test_indices = other.test_indices;
|
|
train_indices = other.train_indices;
|
|
}
|
|
};
|
|
|
|
//+------------------------------------------------------------------+
|
|
//|KFold cross validation: returns the indices for the test/train |
|
|
//|split. |
|
|
//+------------------------------------------------------------------+
|
|
bool kfold(ulong data_length, Folds& folds[],ulong n_splits = 2, bool shuffle = true, int random_seed = 0)
|
|
{
|
|
if(n_splits==0 || data_length<=n_splits)
|
|
{
|
|
Print(__FUNCTION__, " invalid inputs ");
|
|
return false;
|
|
}
|
|
|
|
vector indices = arange(data_length);
|
|
|
|
if(!indices.Size())
|
|
return false;
|
|
|
|
if(shuffle)
|
|
{
|
|
CHighQualityRandStateShell rstate;
|
|
if(random_seed)
|
|
CHighQualityRand::HQRndSeed(random_seed,random_seed+1,rstate.GetInnerObj());
|
|
else
|
|
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
|
|
|
|
indices = shuffleVector(indices,(CHighQualityRandStateShell*)GetPointer(rstate));
|
|
}
|
|
|
|
long start,stop,current,size;
|
|
vector fold_sizes = vector::Zeros(n_splits);
|
|
ArrayResize(folds,(int)n_splits);
|
|
fold_sizes.Fill(floor(data_length/n_splits));
|
|
ulong to = ulong(MathMod(data_length,n_splits));
|
|
vector temp;
|
|
for(ulong i = 0; i<to; ++i)
|
|
fold_sizes[i]+=1.;
|
|
|
|
current = start = stop = size = 0;
|
|
|
|
for(ulong i = 0; i<fold_sizes.Size(); ++i)
|
|
{
|
|
start = current;
|
|
size = long(fold_sizes[i]);
|
|
stop = current + size;
|
|
folds[i].test_indices = sliceVector(indices,start,stop);
|
|
size = start+(long(data_length)-stop);
|
|
if(size)
|
|
{
|
|
folds[i].train_indices = vector::Zeros(size);
|
|
if(start)
|
|
{
|
|
temp = sliceVector(indices,0,start);
|
|
if(!vectorCopy(folds[i].train_indices,temp,0,start))
|
|
return false;
|
|
}
|
|
if(stop<long(data_length))
|
|
{
|
|
temp = sliceVector(indices,stop);
|
|
if(!vectorCopy(folds[i].train_indices,temp,start))
|
|
return false;
|
|
}
|
|
}
|
|
current = stop;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| expanding windows |
|
|
//+------------------------------------------------------------------+
|
|
bool exp_windows(ulong data_length, Folds& folds[],ulong n_splits = 2)
|
|
{
|
|
vector indices = arange(data_length);
|
|
|
|
long test_size = long(data_length/(n_splits+1));
|
|
|
|
ArrayResize(folds,int(n_splits));
|
|
|
|
long test_end, test_start;
|
|
for(ulong i = 0; i<n_splits; ++i)
|
|
{
|
|
test_end = long(data_length)-long(i*test_size);
|
|
test_start = test_end - test_size;
|
|
|
|
folds[n_splits-i-1].train_indices=sliceVector(indices,0,test_start);
|
|
folds[n_splits-i-1].test_indices=sliceVector(indices,test_start,test_end);
|
|
}
|
|
|
|
return true;
|
|
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| The Rosenbrock function. |
|
|
//+-------------------------------------------------------------------+
|
|
double rosen(vector &x)
|
|
{
|
|
vector x1 = sliceVector(x,1);
|
|
vector x2 = sliceVector(x,0,long(x.Size()-1));
|
|
vector r = 100.0*pow(x1-pow(x2,2.0),2.0)+pow(1.-x2,2.0);
|
|
return r.Sum();
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|Gradient of the Rosenbrock function |
|
|
//+------------------------------------------------------------------+
|
|
vector rosen_gradient(vector &x)
|
|
{
|
|
vector grad = vector::Zeros(x.Size());
|
|
if(x.Size()>2)
|
|
{
|
|
vector xm = sliceVector(x,1,long(x.Size()-1));
|
|
vector xm_m1 = sliceVector(x,0,long(x.Size()-2));
|
|
vector xm_p1 = sliceVector(x,2);
|
|
vector temp = (200.0 * (xm - pow(xm_m1,2.0)) - 400.0 * (xm_p1 - pow(xm,2.0)) * xm - 2.0 * (1.0 - xm));
|
|
if(!vectorCopy(grad,temp,1,long(grad.Size()-1)))
|
|
return vector::Zeros(0);
|
|
}
|
|
grad[0] = -400.0 * x[0] * (x[1] - pow(x[0],2.0)) - 2.0 * (1.0 - x[0]);
|
|
grad[grad.Size()-1] = 200.0 * (x[x.Size()-1] - pow(x[x.Size()-2],2.0));
|
|
return grad;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Hessian matrix of the Rosenbrock function. |
|
|
//+------------------------------------------------------------------+
|
|
matrix rosen_hess(vector &x)
|
|
{
|
|
vector x1 = sliceVector(x,0,long(x.Size()-1));
|
|
vector x2 = sliceVector(x,1,long(x.Size()-1));
|
|
vector x3 = vector::Zeros(x.Size());
|
|
matrix h1,h2,h3,H;
|
|
if(!h1.Diag(-400.0*x1,1) || !h2.Diag(400.0*x1,-1))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
H = h1 - h2;
|
|
|
|
x3[0] = 1200.0 * pow(x[0],2.0) - 400.0 * x[1] + 2.0;
|
|
vector temp = 202.0 + 1200.0 * pow(x2,2) - 400.0 * sliceVector(x,2);
|
|
x3[x3.Size()-1] = 200.0;
|
|
|
|
if(!vectorCopy(x3,temp,1,long(x3.Size()-1)))
|
|
return matrix::Zeros(0,0);
|
|
|
|
if(!h3.Diag(x3))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return matrix::Zeros(0,0);
|
|
}
|
|
|
|
return H+h3;
|
|
}
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|