5633 lines
372 KiB
MQL5
5633 lines
372 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); // Key exceeds maximum value; return out-of-bounds upper limit
|
|
}
|
|
else
|
|
if(key < arr[0])
|
|
{
|
|
return -1; // Key is lower than minimum value; return out-of-bounds lower limit
|
|
}
|
|
|
|
// Standard binary search convergence loop
|
|
while(imin<=imax)
|
|
{
|
|
mid = (imin+imax)/2; // Safe midpoint extraction (using integer truncation)
|
|
// Check if key sits perfectly between the current midpoint index and the next index
|
|
if(mid<int(arr.Size()-1) && arr[mid]<=key && arr[mid+1]>key)
|
|
return mid;
|
|
else
|
|
{
|
|
if(key<arr[mid])
|
|
imax = mid - 1; // Narrow search range to left partition
|
|
else
|
|
imin = mid + 1; // Narrow search range to right partition
|
|
}
|
|
}
|
|
return ArraySize(arr) - 2; // Fallback bound adjustment boundary indicator
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| 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); // Pre-allocate return MQL5 vector container
|
|
//---
|
|
if(handle_dates && date_column == column_index)
|
|
{
|
|
// Convert standard datetime string format elements into raw Unix epoch timestamp doubles
|
|
for(int i=0; i<size; ++i)
|
|
ret[i] = (double)StringToTime(Arr[i]);
|
|
}
|
|
else
|
|
{
|
|
// Parse numeric string values straight into doubles
|
|
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; // Infer implicit matrix row count via total size and column count
|
|
ArrayResize(Col,rows);
|
|
|
|
int start = 0;
|
|
for(int i=0; i<cols; ++i)
|
|
{
|
|
start = i;
|
|
|
|
if(i != column-1) // Standardize 1-based index input parameter against internal 0-based iteration
|
|
continue;
|
|
else
|
|
for(int j=0; j<rows; ++j)
|
|
{
|
|
Col[j] = Matrix[start]; // Extract element values sequentially
|
|
start += cols; // Advance flat 1D array offset memory address pointer by row stride
|
|
}
|
|
}
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| product for arrays |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
T prodArray(T &in_array[], uint from=0)
|
|
{
|
|
T result = T(1); // Set multiplicative identity base
|
|
for(uint i = from; i<in_array.Size(); ++i)
|
|
result*=in_array[i]; // Cumulatively compute total array product elements
|
|
return result;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//|calculate size of container |
|
|
//+------------------------------------------------------------------+
|
|
long calculatesize(long initial,long stop,long step)
|
|
{
|
|
//---calculate initial mathematical step lengths
|
|
long len = fabs(initial-stop)/fabs(step);
|
|
//---
|
|
if(len*fabs(step) < fabs(initial-stop))
|
|
len += (fabs(initial-stop) - len*fabs(step))/fabs(step); // Ceiling adjustments for uneven remainder segments
|
|
//---
|
|
return len;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| process start parameter |
|
|
//+------------------------------------------------------------------+
|
|
long normalizestart(long start,ulong size)
|
|
{
|
|
// Implements Python-style wrapping rules mapping negative indices onto end arrays safely
|
|
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)
|
|
{
|
|
// Standardizes terminal slicing markers allowing negative wrap limits or hard clamp size thresholds
|
|
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)
|
|
{
|
|
// Bounds safety verification checks protecting slice loops from infinite ranges, empty outputs, or overflows
|
|
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());
|
|
}
|
|
|
|
// Compute localized continuous positional interpolation mapping ranges
|
|
vector aleph = (double(n)*probs+abprobs);
|
|
|
|
vector k = aleph;
|
|
if(!k.Clip(1.0, double(n-1))) // Enforce boundary clamp constraints within valid indexing limits
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(probs.Size());
|
|
}
|
|
|
|
k = floor(k); // Isolate base array structural address segments
|
|
vector gamma = aleph-k; // Isolate remaining fractional weight offset scalars
|
|
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); // Select values at structural lower-bound thresholds
|
|
|
|
ArrayFree(index);
|
|
|
|
if(!np::vecAsArray(k,index))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(probs.Size());
|
|
}
|
|
|
|
vector x2 = np::select(k,index); // Select values at structural upper-bound thresholds
|
|
|
|
return (1.0-gamma)*x1+gamma*x2; // Linearly blend the values using calculated interpolation weights
|
|
}
|
|
|
|
} // namespace internal
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Return arrays representing indices of a grid |
|
|
//+------------------------------------------------------------------+
|
|
void indices(ulong rows, ulong cols, matrix &outRC[])
|
|
{
|
|
ulong N = 2; // Expected dual coordinate dimension tracking layers count (row mapping matrix, col mapping matrix)
|
|
|
|
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); // Broaden row axis coordinates across columns array grid
|
|
v = np::arange(cols);
|
|
outRC[1] = np::repeat_vector_as_rows_cols(v,rows,true); // Broaden column axis coordinates across rows array grid
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| QuickSort core tracking driver algorithm implementation |
|
|
//+------------------------------------------------------------------+
|
|
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 bitwise division by 2 to extract the midpoint index safely
|
|
p_double=in[(first+last)>>1];
|
|
while(i<j)
|
|
{
|
|
// Evaluate partitions matching sort direction directives
|
|
while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
|
|
{
|
|
if(i==in.Size()-1)
|
|
break;
|
|
i++;
|
|
}
|
|
while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
|
|
{
|
|
if(j==0)
|
|
break;
|
|
j--;
|
|
}
|
|
// Swap structural items across left and right sub-regions when elements are misplaced
|
|
if(i<=j)
|
|
{
|
|
t_double=in[i];
|
|
in[i]=in[j];
|
|
in[j]=t_double;
|
|
//--
|
|
i++;
|
|
if(j==0)
|
|
break;
|
|
j--;
|
|
}
|
|
}
|
|
if(first<j)
|
|
quickSort(in,ascending,first,j); // Recursively slice and sort left sub-partitions
|
|
first=i;
|
|
j=last;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| QuickSort routing tracker linking indices mutations synchronously|
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
bool quickSortIndices(vector<T>& _in,bool ascending,long &indices[],long first,long last)
|
|
{
|
|
vector<T> in = _in; // Fork internal local proxy replica tracking arrays
|
|
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 bitwise division by 2 to extract the pivot index
|
|
p_double=in[(first+last)>>1];
|
|
while(i<j)
|
|
{
|
|
while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
|
|
{
|
|
if(i==in.Size()-1)
|
|
break;
|
|
i++;
|
|
}
|
|
while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
|
|
{
|
|
if(j==0)
|
|
break;
|
|
j--;
|
|
}
|
|
if(i<=j)
|
|
{
|
|
//-- swap value elements inside working arrays
|
|
t_double=in[i];
|
|
in[i]=in[j];
|
|
in[j]=t_double;
|
|
//-- mirror layout swaps tracking original structural addresses
|
|
t_int=indices[i];
|
|
indices[i]=indices[j];
|
|
indices[j]=t_int;
|
|
|
|
i++;
|
|
if(j==0)
|
|
break;
|
|
j--;
|
|
}
|
|
}
|
|
if(first<j)
|
|
quickSortIndices(in,ascending,indices,first,j); // Recursively slice left partitions
|
|
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); // Redirect processing straight to the vector-based QuickSort engine
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| sort matrix by columns or rows |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool sort(matrix<T>&in, bool by_rows, bool ascending=true)
|
|
{
|
|
if(by_rows)
|
|
{
|
|
// Iterate sequentially over rows, isolating vector subsets to sort in-place
|
|
for(ulong i = 0; i<in.Rows(); ++i)
|
|
{
|
|
vector row = in.Row(i);
|
|
if(!sort(row,ascending))
|
|
return false;
|
|
in.Row(row,i);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Iterate sequentially over columns, isolating vector subsets to sort in-place
|
|
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 |
|
|
//+------------------------------------------------------------------+
|
|
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 directly to internal working targets in-place
|
|
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; // Shift bits down rightwards by 1 to isolate midpoint boundaries safely
|
|
for(ulong i = 0; i<half; ++i)
|
|
{
|
|
T temp = in[i];
|
|
in[i] = in[size - i]; // Mirror opposite array value positions inwards
|
|
in[size - i] = temp;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| extract a portion a matrix within a range at specific intervals |
|
|
//+------------------------------------------------------------------+
|
|
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 column 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 grid using row/column combinations
|
|
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 |
|
|
//+------------------------------------------------------------------+
|
|
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 column 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 directly inside working inputs in-place
|
|
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 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> sliceMatrixCols(matrix<T> &in, const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
|
|
{
|
|
// Leverages base sliceMatrix function while keeping full range settings intact along row tracks
|
|
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; // Derive column processing boundary limits safely via right shift bits
|
|
|
|
// Swaps columns inwards sequentially toward structural midlines
|
|
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 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> sliceMatrixRows(matrix<T> &in, const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP)
|
|
{
|
|
// Leverages base sliceMatrix function while keeping full range settings intact along column tracks
|
|
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;
|
|
|
|
// Swaps rows inwards sequentially toward structural midlines
|
|
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 cols |
|
|
//+------------------------------------------------------------------+
|
|
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 rows |
|
|
//+------------------------------------------------------------------+
|
|
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 (Generic Index Type Array) |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T, typename P>
|
|
vector<T> select(vector<T> &in,const P &indices[])
|
|
{
|
|
P size = (P)MathMin(in.Size(),indices.Size()); // Safe clamp checking allocation targets
|
|
|
|
// Guard checks protecting structural lookups against out of bounds array values
|
|
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);
|
|
|
|
// Populate return vectors via targeted indexes arrays mapping
|
|
for(P i=0; i<size; ++i)
|
|
out[i] = in[indices[i]];
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//|return vector of arbitrary elements (MQL5 Vector Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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])]; // Explicit scalar type casting to extract coordinate values safely
|
|
|
|
return out;
|
|
}
|
|
//+----------------------------------------------------------------------------------+
|
|
//|return vector with shifted contents and vacated positions filled with custom value|
|
|
//+----------------------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> shiftVector(vector<T> &_in, long shift, T placeholder)
|
|
{
|
|
if(shift == 0)
|
|
return _in;
|
|
//---
|
|
vector<T> shifted = vector<T>::Zeros(_in.Size());
|
|
//--- Fill vector with place holder value
|
|
if(!shifted.Fill(placeholder))
|
|
{
|
|
Print(__FUNCTION__,": Vector fill failed ", GetLastError());
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
//---Copy values
|
|
if(shift>0)
|
|
for(long i = shift; i<long(shifted.Size()); shifted[i] = _in[i-shift], ++i);
|
|
else
|
|
for(long i = 0; i<long(_in.Size())+shift; shifted[i] = _in[i+fabs(shift)], ++i);
|
|
return shifted;
|
|
}
|
|
//+----------------------------------------------------------------------------------+
|
|
//|return vector with shifted contents and vacated positions filled with custom value|
|
|
//+----------------------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> shiftMatrix(matrix<T> &_in, long shift, T placeholder)
|
|
{
|
|
if(shift == 0)
|
|
return _in;
|
|
//---
|
|
matrix<T> shifted = matrix<T>::Zeros(_in.Rows(),_in.Cols());
|
|
//--- Fill vector with place holder value
|
|
if(!shifted.Fill(placeholder))
|
|
{
|
|
Print(__FUNCTION__,": Matrix fill failed ", GetLastError());
|
|
return matrix<T>::Zeros(0,0);
|
|
}
|
|
//---Copy values
|
|
if(shift>0)
|
|
for(long i = shift; i<shifted.Rows(); shifted.Row(_in.Row(i-shift),i), ++i);
|
|
else
|
|
for(long i = 0; i<_in.Rows(); shifted.Row(_in.Row(i+fabs(shift)),i), ++i);
|
|
return shifted;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|return vector of arbitrary elements (ulong Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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 (int Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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 (uint Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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 (long Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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);
|
|
|
|
// Extract specific source column layouts sequentially into target positions
|
|
for(ulong i=0; i<size; ++i)
|
|
out.Col(in.Col(indices[i]),i);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//|return arbitrary rows of a matrix (long Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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());
|
|
|
|
// Extract specific source row layouts sequentially into target positions
|
|
for(ulong i=0; i<size; ++i)
|
|
out.Row(in.Row(indices[i]),i);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//|return arbitrary columns of a matrix (ulong Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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 (ulong Array Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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 (MQL5 Vector Overload) |
|
|
//+------------------------------------------------------------------+
|
|
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); // Convert index values stored inside vector doubles to ulong sizes
|
|
|
|
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 to store normalized indices and size
|
|
long st, sp, stp, size;
|
|
|
|
//--- Convert potentially negative/default boundaries to concrete vector indices
|
|
st = internal::normalizestart(dest_start, copyto.Size());
|
|
sp = internal::normalizestop(dest_stop, copyto.Size());
|
|
stp = dest_step;
|
|
|
|
//--- Ensure parameters conform to limits and step logic
|
|
if(internal::invalidparams(st, sp, stp, copyto.Size()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters ");
|
|
return false;
|
|
}
|
|
|
|
//--- Calculate total number of elements to copy based on step and bounds
|
|
size = internal::calculatesize(st, sp, stp);
|
|
|
|
//--- Ensure source vector has enough elements to satisfy the slice size
|
|
if(copyfrom.Size() < ulong(size))
|
|
{
|
|
Print(__FUNCTION__, " source vector is of insufficient size ");
|
|
return false;
|
|
}
|
|
|
|
//--- Map elements from source to destination using the slicing sequence
|
|
for(long pos = st, index = 0; index < size; ++index, pos += stp)
|
|
copyto[pos] = copyfrom[index];
|
|
|
|
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 for boundaries and sizing
|
|
long st, sp, stp, size;
|
|
|
|
//--- Normalize start, stop and step values for the destination vector
|
|
st = internal::normalizestart(dest_start, copyto.Size());
|
|
sp = internal::normalizestop(dest_stop, copyto.Size());
|
|
stp = dest_step;
|
|
|
|
//--- Validate normalized index bounds
|
|
if(internal::invalidparams(st, sp, stp, copyto.Size()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters ");
|
|
return false;
|
|
}
|
|
|
|
//--- Compute how many targeted index points match the slice criteria
|
|
size = internal::calculatesize(st, sp, stp);
|
|
|
|
//--- Sequentially write the fixed fill_value to the specified step indices
|
|
for(long pos = st, index = 0; index < size; ++index, pos += stp)
|
|
copyto[pos] = fill_value;
|
|
|
|
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)
|
|
{
|
|
//--- Variables tracking row/column properties: bounds, steps, and sizes
|
|
long rst, rsp, rstp, rsize, cst, csp, cstp, csize;
|
|
|
|
//--- Process and normalize row dimensions
|
|
rst = internal::normalizestart(row_dest_start, copyto.Rows());
|
|
rsp = internal::normalizestop(row_dest_stop, copyto.Rows());
|
|
rstp = row_dest_step;
|
|
|
|
if(internal::invalidparams(rst, rsp, rstp, copyto.Rows()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters for row specification ");
|
|
return false;
|
|
}
|
|
rsize = internal::calculatesize(rst, rsp, rstp);
|
|
|
|
//--- Process and normalize column dimensions
|
|
cst = internal::normalizestart(column_dest_start, copyto.Cols());
|
|
csp = internal::normalizestop(column_dest_stop, copyto.Cols());
|
|
cstp = column_dest_step;
|
|
|
|
if(internal::invalidparams(cst, csp, cstp, copyto.Cols()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters for column specification ");
|
|
return false;
|
|
}
|
|
csize = internal::calculatesize(cst, csp, cstp);
|
|
|
|
//--- Verify the incoming matrix is large enough to supply data for calculated bounds
|
|
if(copyfrom.Cols() < ulong(csize) || copyfrom.Rows() < ulong(rsize))
|
|
{
|
|
Print(__FUNCTION__, " source matrix has insufficient number of columns and/or rows ");
|
|
return false;
|
|
}
|
|
|
|
//--- Multi-dimensional iteration over row and column targets to perform submatrix copying
|
|
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];
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Assign selected elements of matrix 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)
|
|
{
|
|
//--- Variables tracking row/column properties: bounds, steps, and sizes
|
|
long rst, rsp, rstp, rsize, cst, csp, cstp, csize;
|
|
|
|
//--- Process and validate row ranges
|
|
rst = internal::normalizestart(row_dest_start, copyto.Rows());
|
|
rsp = internal::normalizestop(row_dest_stop, copyto.Rows());
|
|
rstp = row_dest_step;
|
|
|
|
if(internal::invalidparams(rst, rsp, rstp, copyto.Rows()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters for row specification ");
|
|
return false;
|
|
}
|
|
rsize = internal::calculatesize(rst, rsp, rstp);
|
|
|
|
//--- Process and validate column ranges
|
|
cst = internal::normalizestart(column_dest_start, copyto.Cols());
|
|
csp = internal::normalizestop(column_dest_stop, copyto.Cols());
|
|
cstp = column_dest_step;
|
|
|
|
if(internal::invalidparams(cst, csp, cstp, copyto.Cols()))
|
|
{
|
|
Print(__FUNCTION__, " invalid function parameters for column specification ");
|
|
return false;
|
|
}
|
|
csize = internal::calculatesize(cst, csp, cstp);
|
|
|
|
//--- Run nested loops to flood targeted coordinate blocks with fill_value
|
|
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;
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns double vector generated incrementally over a fixed size |
|
|
//+------------------------------------------------------------------+
|
|
vector arange(const ulong size, double value=0, double step_value=1)
|
|
{
|
|
//--- Initialize output sequence container
|
|
vector out(size);
|
|
|
|
//--- Compute linearly increasing elements starting from base value
|
|
for(ulong i = 0; i < size; ++i, value += step_value)
|
|
out[i] = value;
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns generic type vector bounded by start and stop conditions |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> arange(T start_value, T stop_value, T step_value)
|
|
{
|
|
//--- Calculate the distance between stop and start points
|
|
T numerator = (stop_value - start_value);
|
|
ulong size = ulong((numerator) / step_value);
|
|
|
|
//--- Include trailing fractional segment step if remainder exists
|
|
if(MathMod(numerator, step_value) > T(0))
|
|
size += 1;
|
|
|
|
//--- Instantiate a zeroed tracking vector
|
|
vector<T> out = vector<T>::Zeros(size);
|
|
|
|
//--- Loop-populate sequential steps onto the template typed container
|
|
for(ulong i = 0; i < size; ++i, start_value += step_value)
|
|
out[i] = start_value;
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns vector covering inclusive integer limits up to stop value|
|
|
//+------------------------------------------------------------------+
|
|
vector range(int stop_value, int start_value = 0, int step_value = 1)
|
|
{
|
|
//--- Round up mathematical size calculation to catch non-divisible ranges
|
|
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);
|
|
|
|
//--- Linearly step up double tracking sequence
|
|
for(ulong i = 0; i < ulong(size); ++i, value += step)
|
|
out[i] = value;
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Outputs range sequence via mutable passing of classic array structures |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool arange(T &in_out[], const int size, T value = 0, T step_value = 1)
|
|
{
|
|
//--- Fail immediately if array size is specified as zero
|
|
if(!size)
|
|
return false;
|
|
|
|
//--- Try resizing array reference, checking for allocation failures
|
|
if(ArrayResize(in_out, fabs(size)) != fabs(size))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
//--- Populate step values across linear dynamic allocations
|
|
for(int i = 0; i < size; ++i, value += step_value)
|
|
in_out[i] = value;
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Equivalent to NumPy linspace(), outputs evenly spaced intervals |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> linspace(T start, T stop, ulong num=50)
|
|
{
|
|
vector<T> out(num);
|
|
|
|
//--- Edge Case handling: Return empty vector if no points requested
|
|
if(!num)
|
|
return vector<T>::Zeros(num);
|
|
|
|
//--- Edge Case handling: If 1 point, matching NumPy convention, use start value
|
|
if(num == 1)
|
|
{
|
|
out[0] = start;
|
|
return out;
|
|
}
|
|
|
|
//--- Identify physical step size (chunk delta) per element gap
|
|
T chunk = (stop - start) / T(num - 1);
|
|
|
|
//--- Construct discrete segments relative to baseline position
|
|
for(ulong i = 0; i < out.Size(); ++i)
|
|
out[i] = start + (T(i) * chunk);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Equivalent to NumPy logspace(), logs linear space out to custom bases|
|
|
//+------------------------------------------------------------------+
|
|
vector<double> logspace(double start, double stop, ulong num=50, double base = 10.0)
|
|
{
|
|
//--- Generate base exponent ranges linearly using internal linspace
|
|
vector<double> out = linspace(start, stop, num);
|
|
|
|
//--- Elevate logarithmic scales per value index transformation
|
|
for(ulong i = 0; i < out.Size(); ++i)
|
|
out[i] = pow(base, out[i]);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Fill main diagonal of square matrix with standard value |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool fillDiagonal(matrix<T> &in_out, T fill_value)
|
|
{
|
|
//--- Enforce square matrix condition constraints
|
|
if(in_out.Rows() != in_out.Cols())
|
|
{
|
|
Print(__FUNCTION__ "Invalid Input. Matrix is not square ");
|
|
return false;
|
|
}
|
|
|
|
//--- Traverse uniform index coordinate paths [i][i]
|
|
for(ulong i = 0; i < in_out.Rows(); ++i)
|
|
in_out[i][i] = fill_value;
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Copy an array of matrices (simulating a 3D structural container) |
|
|
//+------------------------------------------------------------------+
|
|
bool copy3D(matrix ©_from[], matrix ©_to[])
|
|
{
|
|
//--- Realign destination tracking array allocations to match donor length
|
|
if(ArraySize(copy_to) != ArraySize(copy_from))
|
|
ArrayResize(copy_to, (int)copy_from.Size());
|
|
|
|
//--- Individually pass reference matrix blocks over the layer arrays
|
|
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)
|
|
{
|
|
//--- Default column sizing relies entirely on incoming reference limits
|
|
ulong n = num_columns ? num_columns : in.Size();
|
|
|
|
matrix<T> out(in.Size(), n);
|
|
|
|
//--- Populate matrix column vectors with exponential powers of input terms
|
|
for(ulong i = 0; i < n; ++i)
|
|
out.Col(pow(in, ascending ? i : n - 1 - i), i);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generic vector extraction to a raw typed data array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool vecAsArray(const vector<T> &in, T &out[], uint srce_start = 0, uint count = 0)
|
|
{
|
|
//--- Validate input data presence
|
|
if(in.Size() < 1)
|
|
{
|
|
Print(__FUNCTION__, " Empty vector");
|
|
return false;
|
|
}
|
|
|
|
//--- Calculate target footprint slice
|
|
if(!count)
|
|
count = uint(in.Size());
|
|
|
|
//--- Allocate backing arrays to accept target conversion sequences
|
|
if(out.Size() != count && ArrayResize(out, int(count)) < int(count))
|
|
{
|
|
Print(__FUNCTION__, " resize error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
//--- Iteratively populate matching element indices via explicit cast conversions
|
|
for(uint i = srce_start; i < (srce_start + count); ++i)
|
|
out[i - srce_start] = T(in[i]);
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Type-specific specialization: Vector to Double 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Type-specific specialization: Vector to Integer 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Type-specific specialization: Vector to Unsigned Int 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Type-specific specialization: Vector to Long Integer 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Type-specific specialization: Vector to Unsigned Long 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 flattening wrapper to 1-D Double Array |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T, typename S>
|
|
bool matAsArray(const matrix<T> &in, double &out[])
|
|
{
|
|
//--- Confirm dimensional sizing viability across matrix fields
|
|
if(in.Cols() + in.Rows() < 1)
|
|
{
|
|
Print(__FUNCTION__, " Empty vector");
|
|
return false;
|
|
}
|
|
|
|
//--- Validate flat linear destination scaling space matches rows * columns
|
|
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;
|
|
}
|
|
|
|
//--- Stream sequentially flattened structural offsets directly using .Flat() indices
|
|
for(uint i = 0; i < out.Size(); ++i)
|
|
out[i] = double(in.Flat(i));
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generic matrix flattening wrapper to 1-D Integer 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 flattening wrapper to 1-D Unsigned Int 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 flattening wrapper to 1-D Long Integer 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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Least squares polynomial fit (Curve Fitting Functionality) |
|
|
//+------------------------------------------------------------------+
|
|
vector polyfit(vector& x, vector& y, int deg)
|
|
{
|
|
int order = deg + 1;
|
|
|
|
//--- Verify error validation criteria constraints
|
|
if(deg < 0)
|
|
{
|
|
Print(__FUNCTION__, ": Param deg should be >= 0.");
|
|
return vector::Zeros(0);
|
|
}
|
|
if(!x.Size())
|
|
{
|
|
Print(__FUNCTION__, ": Param x is empty.");
|
|
return vector::Zeros(0);
|
|
}
|
|
if(x.Size() != y.Size())
|
|
{
|
|
Print(__FUNCTION__, ": Params x and y not of equal size");
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
//--- Generate left-hand side values utilizing a Vandermonde structural matrix
|
|
matrix lhs = vander(x, ulong(order));
|
|
vector rhs = y;
|
|
|
|
//--- Normalize scale parameters to guard algorithmic stability against floating anomalies
|
|
vector scale = sqrt((lhs * lhs).Sum(0));
|
|
lhs = divide(lhs, scale);
|
|
|
|
//--- Solve regularized expressions applying native Least Squares algorithm (.LstSq)
|
|
vector c = lhs.LstSq(rhs);
|
|
c /= scale; // Reverse scaling multiplier offset adjustment factors
|
|
return c;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generic matrix flattening wrapper to 1-D Unsigned Long 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 calculation order degree |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> diff(vector<T> &in, ulong degree = 1)
|
|
{
|
|
//--- Out of bounds/invalid param check
|
|
if(in.Size() < 1 || degree < 1 || in.Size() < degree + 1)
|
|
{
|
|
Print(__FUNCTION__, " Invalid parameter ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
|
|
vector<T> yy, zz;
|
|
|
|
//--- Generate matching slices offset by degree to evaluate consecutive deltas
|
|
yy = sliceVector(in, long(degree));
|
|
zz = sliceVector(in, BEGIN, long(in.Size() - degree));
|
|
|
|
//--- Return differential subtraction matrices directly
|
|
return yy - zz;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Difference matrix computed independently relative to rows/columns|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> diff(matrix<T> &in, ulong degree = 1, bool by_row = true)
|
|
{
|
|
//--- Parameter bounds structural check validation
|
|
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;
|
|
|
|
//--- Slice configuration mapped relative to active evaluation pathways (row vs col)
|
|
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 expanding copies of original internal elements |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> repeat_vector_elements(vector<T> &vect, ulong num_repetitions)
|
|
{
|
|
//--- Allocate output matrix sized exactly to handle expanded tracking allocations
|
|
vector<T> out(num_repetitions * vect.Size());
|
|
|
|
//--- Repeatedly clone structural sets using core vector copy subroutines
|
|
for(ulong i = 0; i < out.Size(); i += vect.Size())
|
|
vectorCopy(out, vect, i, i + vect.Size());
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generate vector matching multi-reps of internal 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());
|
|
|
|
//--- Sequentially step through grid layers flattening individual row segments
|
|
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 broadcasting vectors along matching grid spaces |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> repeat_vector_as_rows_cols(vector<T> &in, ulong num_repetitions, bool as_rows = true)
|
|
{
|
|
//--- Calculate inverted dimensions mapping configuration targets
|
|
matrix<T> out(as_rows ? num_repetitions : in.Size(), as_rows ? in.Size() : num_repetitions);
|
|
|
|
//--- Direct matrix dimension line assignment mapping loop
|
|
for(ulong i = 0; i < num_repetitions; ++i)
|
|
as_rows ? out.Row(in, i) : out.Col(in, i);
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Generate matrix mirroring structured lines of raw input matrixes|
|
|
//+------------------------------------------------------------------+
|
|
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);
|
|
}
|
|
|
|
//--- Pre-allocate memory reflecting global row/col repetition factor dimensions
|
|
matrix<T> out((by_rows) ? num_repetitions * in.Rows() : in.Rows(), (!by_rows) ? num_repetitions * in.Cols() : in.Cols());
|
|
|
|
//--- Block map assignments matching indexed layout coordinates
|
|
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 variable scaled matrix sequences matching discrete array maps|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> repeat_matrix_rows_cols(matrix<T> &in, ulong &reps[], bool by_rows)
|
|
{
|
|
ulong max = 0;
|
|
|
|
//--- Aggregate compound structural element spans
|
|
for(uint i = 0; i < reps.Size(); ++i)
|
|
max += reps[i];
|
|
|
|
//--- Enforce input size conformance constraints against raw processing limits
|
|
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());
|
|
|
|
//--- Segment-slice elements processing arbitrary variable dimension arrays
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Helper wrapper converting vector into an expanded column matrix |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> vectorAsColumnMatrix(vector<T> &in, ulong num_cols)
|
|
{
|
|
return repeat_vector_as_rows_cols(in, num_cols, false);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Helper wrapper converting vector into an expanded row matrix |
|
|
//+------------------------------------------------------------------+
|
|
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 tracker array to flag processed indices
|
|
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 vector with maximum possible unique size
|
|
vector<T>out(size);
|
|
//--- find unique elements
|
|
ulong unique_count=0;
|
|
T value=0;
|
|
while(true)
|
|
{
|
|
bool flag=false;
|
|
// Loop through the vector starting from the current unique item baseline
|
|
for(ulong i=unique_count; i<size; ++i)
|
|
{
|
|
// If a new unseen element is found, record it as a unique candidate
|
|
if(!flag && !checked[i])
|
|
{
|
|
value=in[i];
|
|
out[unique_count]=in[i];
|
|
++unique_count;
|
|
checked[i]=true;
|
|
flag=true; // Set flag to mark that we are scanning for duplicates of this value
|
|
}
|
|
else
|
|
// If the current element matches our active unique value, mark it checked
|
|
if(flag && value==in[i])
|
|
checked[i]=true;
|
|
}
|
|
// If no new unique item was encountered during the pass, extraction is complete
|
|
if(!flag)
|
|
break;
|
|
}
|
|
//--- resize target vector to the exact number of unique elements found
|
|
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 occurrences 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 matching maximum possible allocation sizes
|
|
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 along with their first occurrences and frequency counts
|
|
ulong unique_count=0;
|
|
T value=0;
|
|
while(true)
|
|
{
|
|
bool flag=false;
|
|
for(ulong i=unique_count; i<size; ++i)
|
|
{
|
|
// Extract the next unvisited item
|
|
if(!flag && !checked[i])
|
|
{
|
|
value=in[i];
|
|
out[unique_count]=in[i];
|
|
indices[unique_count]=long(i); // Trace back the index position of original item
|
|
++unique_count;
|
|
checked[i]=true;
|
|
flag=true;
|
|
}
|
|
else
|
|
// If a duplicate is encountered later in the vector, update its frequency counter
|
|
if(flag && value==in[i])
|
|
{
|
|
checked[i]=true;
|
|
++counts[unique_count-1];
|
|
}
|
|
}
|
|
if(!flag)
|
|
break;
|
|
}
|
|
//--- Downsize and fit all resulting arrays to match exact unique count
|
|
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)
|
|
{
|
|
// Sum cross-sections along the axis to act as a footprint signature proxy for rows/columns
|
|
vector sum = in.Sum(axis);
|
|
vector spec;
|
|
long indices[],counts[];
|
|
// Find uniqueness positions based on structural proxy data signatures
|
|
if(!unique(sum,spec,indices,counts))
|
|
{
|
|
Print(__FUNCTION__,": error ");
|
|
return matrix<T>::Zeros(0,0);
|
|
}
|
|
// Select and form the unique output sub-matrix based on extracted valid axis indices
|
|
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 specific coordinates are not given, generate evenly distributed steps using dx
|
|
if(x.Size()==0)
|
|
{
|
|
d = vector::Zeros(y.Size());
|
|
d.Fill(dx);
|
|
}
|
|
else
|
|
d = diff(x); // Calculate step increments between points
|
|
|
|
// Compute area boundaries: step * average heights of neighboring trapezoid segments
|
|
vector<T> ret = (d*(sliceVector(y,1) + sliceVector(y,0,-1)))/2.0;
|
|
|
|
// Sum up individual slices to return total definite integration value
|
|
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());
|
|
|
|
// Check if a periodic wrap (e.g. angle wrap-around interpolation) is specified
|
|
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));
|
|
|
|
// Wrap input values inside the bounded range using modulo
|
|
x = MathMod(xv,vp);
|
|
xp = MathMod(xpoints,vp);
|
|
|
|
long indices[];
|
|
|
|
// Generate sequence range and sort internal parameters on period wrapper
|
|
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);
|
|
|
|
// Scale buffers up to handle repeating boundary edge conditions
|
|
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[];
|
|
|
|
// Copy point references onto a basic array structure for binary exploration lookup
|
|
if(!vecAsArray(xpoints,xvals) || !ArraySort(xvals))
|
|
{
|
|
Print(__FUNCTION__, " ", __LINE__," error ", GetLastError());
|
|
return out;
|
|
}
|
|
|
|
// Define edge fallback values if input elements fall entirely out of range boundaries
|
|
lval=(left==DBL_MIN)?ypoints[0]:left;
|
|
rval=(right==DBL_MAX)?ypoints[ypoints.Size()-1]:right;
|
|
|
|
// Handle special edge configuration where only a single data point mapping exists
|
|
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
|
|
{
|
|
// Calculate slopes beforehand if points structural size accommodates them
|
|
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; // Pre-calculated rise over run variations
|
|
}
|
|
|
|
int j = 0;
|
|
// Loop over every target interpolation coordinate element
|
|
for(ulong i = 0; i<xv.Size(); ++i)
|
|
{
|
|
if(MathClassify(xv[i])==FP_NAN)
|
|
out[i] = xv[i];
|
|
|
|
// Use binary search to locate neighborhood bounds inside known map intervals
|
|
j = internal::binary_search(xv[i],xvals);
|
|
if(j == -1)
|
|
out[i] = lval; // Falls short of left edge bounds
|
|
else
|
|
if(j == int(xpoints.Size()))
|
|
out[i] = rval; // Extends past right edge bounds
|
|
else
|
|
if(j == int(xpoints.Size()-1))
|
|
out[i] = ypoints[j];
|
|
else
|
|
if(xpoints[j] == xv[i])
|
|
out[i] = ypoints[j]; // Exact match found
|
|
else
|
|
{
|
|
// General continuous point-to-point step slope computation
|
|
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];
|
|
|
|
// Edge validation fallback verification check on NaN conditions
|
|
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);
|
|
|
|
// Validation rule check: Columns count of factor matrix A must match row count of factor matrix B
|
|
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);
|
|
|
|
// Perform manual standard triple nested loop cross-dot-product calculations
|
|
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)
|
|
{
|
|
// Adapt the vector input format structurally into a single column matrix layout
|
|
matrix<T> matrix_b_mod = matrix<T>::Zeros(vector_b.Size(),1);
|
|
matrix_b_mod.Col(vector_b,0);
|
|
|
|
// Call standard matrix-to-matrix multiplication handler
|
|
matrix<T>out = matmul(matrix_a,matrix_b_mod);
|
|
return out.Col(0); // Return result converted back as a flattened flat vector array
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| 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;
|
|
|
|
// Verify array profile data orientation before launching comparative scans
|
|
if(!np::sort(samplesorted,sorted[0]<=sorted[sorted.Size()-1]))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return false;
|
|
}
|
|
|
|
// Search insertions tracking boundaries item by item
|
|
for(ulong j = 0; j<sample.Size(); ++j)
|
|
{
|
|
for(ulong i = 1; i<sorted.Size(); ++i)
|
|
{
|
|
// Verify bounded constraints based on right-side insertion parameter directives
|
|
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
|
|
{
|
|
// Boundary fallback edge cases: value belongs lower than the absolute minimum
|
|
if((!right_side && sample[j]<=sorted.Min()) || (right_side && sample[j]<sorted.Min()))
|
|
{
|
|
out[j] = 0;
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
// Boundary fallback edge cases: value belongs higher than the absolute maximum
|
|
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)
|
|
{
|
|
// Broadcaster layout safety constraint check
|
|
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());
|
|
|
|
// Apply vector addition row-by-row if vector dimensions line up with columns count
|
|
if(vec.Size() == mat.Cols())
|
|
for(ulong i = 0; i<mat.Rows(); ++i)
|
|
out.Row(mat.Row(i)+vec,i);
|
|
else
|
|
{
|
|
// Apply vector addition column-by-column if matching up against row count dimensions
|
|
out=matrix::Zeros(mat.Rows(),vec.Size());
|
|
for(ulong i = 0; i<out.Cols(); ++i)
|
|
out.Col(mat.Col(0)+vec[i],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;
|
|
// Subtract vector values element-by-element across matrix row profiles
|
|
if(vec.Size() == mat.Cols())
|
|
{
|
|
out = matrix::Zeros(mat.Rows(),mat.Cols());
|
|
for(ulong i = 0; i<mat.Rows(); ++i)
|
|
out.Row(mat.Row(i)-vec,i);
|
|
}
|
|
else
|
|
{
|
|
// Subtract vector values element-by-element across matrix column profiles
|
|
out = matrix::Zeros(mat.Rows(),vec.Size());
|
|
for(ulong i = 0; i<out.Cols(); ++i)
|
|
out.Col(mat.Col(0)-vec[i],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());
|
|
// Inverse dimension order subtraction (Vector - Matrix elements) along rows orientation
|
|
if(vec.Size()==mat.Cols())
|
|
for(ulong i = 0; i<mat.Rows(); ++i)
|
|
out.Row(vec-mat.Row(i),i);
|
|
else
|
|
{
|
|
// Inverse dimension order subtraction along columns orientation
|
|
out=matrix::Zeros(mat.Rows(),vec.Size());
|
|
for(ulong i = 0; i<out.Cols(); ++i)
|
|
out.Col(vec[i]-mat.Col(0),i);
|
|
}
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| vector subtraction with broadcasting |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> minus(vector<T>& va, vector<T>& vb)
|
|
{
|
|
// Direct subtraction shortcut if structural sizes are identical
|
|
if(va.Size()==vb.Size())
|
|
return va-vb;
|
|
else
|
|
{
|
|
// If sizes differ, one must strictly be a scalar singleton (Size = 1) to enable broadcasting
|
|
if(va.Size()!=1 && vb.Size()!=1)
|
|
{
|
|
Print(__FUNCTION__, " operands could not be broadcast together ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
else
|
|
{
|
|
// Subtract scalar va from array vb or scalar vb from array va depending on sizes layout
|
|
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 multiplication 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());
|
|
// Scalar-by-vector multiplication applied uniformly across structural rows or columns layers
|
|
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());
|
|
// Divide structural matrix values uniformly by the broadcaster vector values profile
|
|
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());
|
|
// Inverted dimension order division (Vector elements divided by Matrix cells) row/col loops
|
|
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;
|
|
|
|
// Convert raw dimension limits into manageable vectors tracking format lengths
|
|
if(!vdims.Assign(dims))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(1);
|
|
}
|
|
|
|
// Extract stride steps configuration specifications
|
|
vector<T> vec = np::sliceVector(vdims,1);
|
|
|
|
if(!np::reverseVector(vec))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(1);
|
|
}
|
|
|
|
// Compute cumulative products of shapes to define sub-dimensional stride offsets lengths
|
|
vec = vec.CumProd();
|
|
|
|
if(!np::reverseVector(vec))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(1);
|
|
}
|
|
|
|
// Expand layout mapping tracker buffer to account for the leading dimensionality level
|
|
vector<T> nvec = vector::Ones(vec.Size()+1);
|
|
|
|
if(!vectorCopy(nvec,vec,1))
|
|
{
|
|
Print(__FUNCTION__, " error ", GetLastError());
|
|
return vector::Zeros(1);
|
|
}
|
|
|
|
// Initialize data conversion matrices targeting multi-index parameters dimensions
|
|
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);
|
|
}
|
|
|
|
// Use Matrix Multiply (Dot Product) to calculate total linear base offsets values
|
|
return nvec.MatMul(out);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| get the bin edges |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> binEdges(vector<T>&in, ulong numbins=3)
|
|
{
|
|
// Determine range limits boundaries across the entire input cluster vector
|
|
T first_edge = in.Min();
|
|
T last_edge = in.Max();
|
|
|
|
// Construct linear evenly-spaced coordinate boundaries arrays matching numbins scale
|
|
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 |
|
|
//| Computes the bitwise NOT operation on elements of a vector. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseNot(vector<T>&left)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Input validation: ensure the vector contains elements
|
|
if(left.Size()<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function input ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Convert the MQL vector object into a standard array for Math bitwise operations
|
|
if(!vecAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Execute the internal MQL mathematical bitwise NOT operation
|
|
if(!MathBitwiseNot(aleft,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Assign the resulting native array back into the native vector object
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseAnd |
|
|
//| Computes the bitwise AND operation between two vectors. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseAnd(vector<T>&left,vector<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Validation: Vectors must be of identical size and not empty
|
|
if(left.Size()!=right.Size() || right.Size()<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Convert both input vectors into native array types
|
|
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Perform element-wise bitwise AND
|
|
if(!MathBitwiseAnd(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseAnd failed ", GetLastError()); // Fixed log message to match operation
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Re-populate and return the vector object
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseOr |
|
|
//| Computes the bitwise OR operation between two vectors. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseOr(vector<T>&left,vector<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Validation: Dimension matching and size boundary check
|
|
if(left.Size()!=right.Size() || right.Size()<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Convert vector objects into internal raw arrays
|
|
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Execute the mathematical bitwise OR function
|
|
if(!MathBitwiseOr(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseOr failed ", GetLastError()); // Fixed log string
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Re-assign raw array data back onto the generic vector
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseXor |
|
|
//| Computes the bitwise XOR operation between two vectors. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseXor(vector<T>&left,vector<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Validation: Match dimensions and check for empty elements
|
|
if(left.Size()!=right.Size() || right.Size()<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Convert matrix vector wrappers to system arrays
|
|
if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size()); // Added missing semicolon
|
|
}
|
|
|
|
// Run the mathematical bitwise Exclusive OR (XOR)
|
|
if(!MathBitwiseXor(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseXor failed ", GetLastError()); // Fixed log string
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Repack results array into a generic MQL vector
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseShiftL |
|
|
//| Shifts bits to the left for all vector elements. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseShiftL(vector<T>&left,int shift)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Validation: Vector must have contents and shift amount must be positive
|
|
if(!left.Size() || shift<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Map vector to temporary calculations array
|
|
if(!vecAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size()); // Added missing semicolon
|
|
}
|
|
|
|
// Call internal bit shifting calculation
|
|
if(!MathBitwiseShiftL(aleft,shift,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseShiftL failed ", GetLastError()); // Fixed log string
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Assign result buffer back to structural vector
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseShiftR |
|
|
//| Shifts bits to the right for all vector elements. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T>bitwiseShiftR(vector<T>&left,int shift)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Validation: Check array bounds and enforce positive shift count
|
|
if(!left.Size() || shift<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
// Extract flat layout array from vector template object
|
|
if(!vecAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " vector conversion failed ");
|
|
return vector::Zeros(left.Size()); // Added missing semicolon
|
|
}
|
|
|
|
// Compute systemic right-shift math logic
|
|
if(!MathBitwiseShiftR(aleft,shift,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseShiftR failed ", GetLastError()); // Fixed log string
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
vector<T> vec;
|
|
|
|
// Map array primitives back into your vector struct object
|
|
if(!vec.Assign(res))
|
|
{
|
|
Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
|
|
return vector::Zeros(left.Size());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| BitwiseNot (Matrix Overload) |
|
|
//| Computes the bitwise NOT operation element-wise on a matrix. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseNot(matrix<T>&left)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Structure Integrity Verification: Reject flat matrices
|
|
if(!left.Rows() || !left.Cols())
|
|
{
|
|
Print(__FUNCTION__, " invalid function input ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Flatten the 2D matrix into a 1D continuous array format
|
|
if(!matAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Run bitwise processing on array representation
|
|
if(!MathBitwiseNot(aleft,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Reconstruct the array output back into matrix dimensionality structures
|
|
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 (Matrix Overload) |
|
|
//| Computes the bitwise AND operation element-wise on two matrices. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseAnd(matrix<T>&left,matrix<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Bound check validations across both matrices axes
|
|
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Convert multi-dimensional templates down to processing raw arrays
|
|
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Compute element-by-element intersection bits
|
|
if(!MathBitwiseAnd(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseAnd failed ", GetLastError()); // Fixed log string
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Map flattened operational logs into target output dimensions
|
|
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 (Matrix Overload) |
|
|
//| Computes the bitwise OR operation element-wise on two matrices. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseOr(matrix<T>&left,matrix<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Validation checks for size bounds
|
|
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Flatten processing structures
|
|
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Run MQL bit math block logic
|
|
if(!MathBitwiseOr(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseOr failed ", GetLastError()); // Fixed log string
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Copy layout matrix back onto instance object
|
|
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 (Matrix Overload) |
|
|
//| Computes the bitwise XOR operation element-wise on two matrices. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseXor(matrix<T>&left,matrix<T>&right)
|
|
{
|
|
int aleft[];
|
|
int aright[];
|
|
int res[];
|
|
|
|
// Matrix non-empty structure assertions
|
|
if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Extract raw values
|
|
if(!matAsArray(left,aleft) || !matAsArray(right,aright))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Perform XOR math processing calculation
|
|
if(!MathBitwiseXor(aleft,aright,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseXor failed ", GetLastError()); // Fixed log string
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Build spatial coordinate tracking fields back into the result container
|
|
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 (Matrix Overload) |
|
|
//| Shifts bits left for all matrix elements. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseShiftL(matrix<T>&left,int shift)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Ensure valid boundaries and positive values
|
|
if(!left.Rows() || !left.Cols() || shift<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Convert matrix format down into array buffers
|
|
if(!matAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Execute left-shift execution pipeline
|
|
if(!MathBitwiseShiftL(aleft,shift,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseShiftL failed ", GetLastError()); // Fixed log string
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Reshape into original grid profile constraints
|
|
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 (Matrix Overload) |
|
|
//| Shifts bits right for all matrix elements. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T>bitwiseShiftR(matrix<T>&left,int shift)
|
|
{
|
|
int aleft[];
|
|
int res[];
|
|
|
|
// Structural parameter sanitization check
|
|
if(!left.Rows() || !left.Cols() || shift<1)
|
|
{
|
|
Print(__FUNCTION__, " invalid function inputs ");
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
// Convert matrix wrapper instance properties into native array blocks
|
|
if(!matAsArray(left,aleft))
|
|
{
|
|
Print(__FUNCTION__, " matrix conversion failed ");
|
|
return matrix::Zeros(left.Rows(),left.Cols()); // Added missing semicolon
|
|
}
|
|
|
|
// Perform binary structural shift calculation
|
|
if(!MathBitwiseShiftR(aleft,shift,res))
|
|
{
|
|
Print(__FUNCTION__ " MathBitwiseShiftR failed ", GetLastError()); // Fixed log string
|
|
return matrix::Zeros(left.Rows(),left.Cols());
|
|
}
|
|
|
|
matrix<T> mat;
|
|
|
|
// Reconstruct geometric configuration layouts matching the parent element dimensions
|
|
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 |
|
|
//| Outputs content from a generic matrix instance to filesystem log |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
bool matrix2csv(string csv_name, matrix<T> &matrix_, string header_string="", bool common=false, int digits=8)
|
|
{
|
|
// Erase pre-existing file configurations safely before writing new assets
|
|
FileDelete(csv_name);
|
|
int handle = FileOpen(csv_name,FILE_WRITE|FILE_CSV|FILE_ANSI|(common?FILE_COMMON:FILE_ANSI),",",CP_UTF8);
|
|
|
|
// Generate default programmatic numeric indexing string headers if left blank
|
|
if(header_string == "")
|
|
for(ulong i=0; i<matrix_.Cols(); ++i)
|
|
header_string += IntegerToString(long(i)) + (i==matrix_.Cols()-1?"":",");
|
|
|
|
// File access error checking routines
|
|
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 = StringGetCharacter(",",0);
|
|
|
|
// Tokenize structural headers using selected formatting character split configurations
|
|
StringSplit(header_string,u_sep, header);
|
|
|
|
vector<T> colsinrows = matrix_.Row(0);
|
|
|
|
// Verify parsed configuration strings match explicit geometric column widths
|
|
if(ArraySize(header) != (int)colsinrows.Size())
|
|
{
|
|
printf("headers=%d and columns=%d from the matrix vary is size ",ArraySize(header),colsinrows.Size());
|
|
FileClose(handle); // Added tracking file close protection execution
|
|
return false;
|
|
}
|
|
|
|
// Construct clean file header string outputs
|
|
string header_str = "";
|
|
for(int i=0; i<ArraySize(header); ++i)
|
|
header_str += header[i] + (i+1 == colsinrows.Size() ? "" : ",");
|
|
|
|
// Write string metadata configuration headers directly onto output file streams
|
|
FileWrite(handle,header_str);
|
|
FileSeek(handle,0,SEEK_SET);
|
|
|
|
// Parse matrix element details block rows safely avoiding infinite lockups
|
|
for(ulong i=0; i<matrix_.Rows() && !IsStopped(); ++i)
|
|
{
|
|
ZeroMemory(concstring);
|
|
row = matrix_.Row(i);
|
|
|
|
// Loop individual matrix columns element locations within tracked row subsets
|
|
for(ulong j=0, cols =1; j<row.Size() && !IsStopped(); ++j, ++cols)
|
|
{
|
|
current_time = GetTickCount();
|
|
// Sanitize and serialize numeric metrics according to maximum floating-precision parameter controls
|
|
concstring += (string)NormalizeDouble(row[j],digits) + (cols == matrix_.Cols() ? "" : ",");
|
|
}
|
|
|
|
// Append row records safely directly matching file boundary end locations
|
|
FileSeek(handle,0,SEEK_END);
|
|
FileWrite(handle,concstring);
|
|
}
|
|
|
|
FileClose(handle);
|
|
return (true);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Read CSV file into matrix |
|
|
//| Parses local flat tables and structures inside structured metrics|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Open files ensuring data-sharing reading protection is allowed concurrently
|
|
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;
|
|
|
|
// Read target strings to terminal locations cleanly safely checking state conditions
|
|
while(!FileIsEnding(handle) && !IsStopped())
|
|
{
|
|
string data = FileReadString(handle);
|
|
|
|
// Intercept first-line elements to dynamically read text schema mapping definitions
|
|
if((rows == 0 && skip_header) || (rows == 1 && !skip_header))
|
|
{
|
|
ArrayResize(csv_header,column+1);
|
|
csv_header[column] = skip_header?data:string(column);
|
|
}
|
|
|
|
++column;
|
|
|
|
// Keep tracked item positions outside basic header boundary contexts
|
|
if(rows>0)
|
|
{
|
|
++all_size;
|
|
ArrayResize(Arr,all_size);
|
|
Arr[all_size-1] = data;
|
|
}
|
|
|
|
// Monitor tracking checks across string character transitions points
|
|
if(FileIsLineEnding(handle))
|
|
{
|
|
cols_total=column;
|
|
++rows;
|
|
column = 0;
|
|
}
|
|
}
|
|
|
|
FileClose(handle);
|
|
}
|
|
|
|
// Compute total matrix properties matching calculations mapping arrays
|
|
int rows = all_size/cols_total;
|
|
Comment("");
|
|
|
|
matrix mat(rows, cols_total);
|
|
string Col[];
|
|
|
|
// Distribute parsed cell text data back across vertical structural column segments
|
|
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 data from string into matrix |
|
|
//| Converts raw raw multi-line strings into formatted matrices |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Segment parent structural block layouts directly at endline boundary nodes
|
|
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;
|
|
|
|
// Inspect base layout patterns using character split configurations
|
|
cols_total = StringSplit(lines[shift],StringGetCharacter(delimiter,0),csv_header);
|
|
|
|
// Clean blank parameters inside calculations arrays safely
|
|
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);
|
|
|
|
// Re-populate inner arrays with matching parameters
|
|
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[];
|
|
|
|
// Map elements vertically to structure target matrix layouts
|
|
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);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sample Vector |
|
|
//| Selects a uniform random sample from the vector elements |
|
|
//| with replacement using the basic MathRand generator. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool sampleVector(vector<T> &in,const ulong count,vector<T> &result)
|
|
{
|
|
if(!in.Size())
|
|
return false;
|
|
|
|
// Ensure target allocation space is sized adequately
|
|
if(result.Size()<count)
|
|
if(!result.Resize(count))
|
|
return(false);
|
|
|
|
// Populate target records using pseudo-random indices
|
|
for(ulong i=0; i<count; ++i)
|
|
{
|
|
int ind=(int)(in.Size()*MathRand()/32768);
|
|
result[i]=in[ind];
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Choice |
|
|
//| Generates a random sample from a given 1-D vector with options |
|
|
//| for custom probability distributions and non-replacement sampling|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> choice(vector<T>& a, vector& p, ulong size = 0, bool replace = true)
|
|
{
|
|
ulong pop_size = a.Size();
|
|
vector idx;
|
|
|
|
// Structural boundaries validation
|
|
if(!pop_size && size)
|
|
{
|
|
Print(__FUNCTION__, " a cannot be empty unless no samples are taken ");
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
|
|
// Custom distribution setup and validation routines
|
|
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);
|
|
}
|
|
}
|
|
|
|
// Initialize the high-quality crypto-grade pseudo-random state structure objects
|
|
CHighQualityRandStateShell rstate;
|
|
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
|
|
|
|
bool is_scalar = (size==0);
|
|
ulong shape = 0;
|
|
vector cdf;
|
|
|
|
if(!is_scalar)
|
|
shape = size;
|
|
else
|
|
size = 1;
|
|
|
|
// CASE 1: Sampling with Replacement allowed
|
|
if(replace)
|
|
{
|
|
if(p.Size())
|
|
{
|
|
// Create a Cumulative Distribution Function array track matching probabilities
|
|
cdf = p.CumSum();
|
|
cdf/=cdf[cdf.Size()-1]; // Normalize cumulative tracks
|
|
vector uniform_samples = vector::Zeros(shape);
|
|
|
|
// Gather uniform random values between 0.0 and 1.0
|
|
for(ulong i = 0; i<shape; uniform_samples[i] = CHighQualityRand::HQRndUniformR(rstate.GetInnerObj()), ++i);
|
|
|
|
// Map random values inside calculated ranges using structural binary search lookups
|
|
if(!searchsorted(cdf,uniform_samples,true,idx))
|
|
{
|
|
Print(__FUNCTION__, " searchsorted failed ", __LINE__);
|
|
return vector<T>::Zeros(0);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Default path: Generate uniform completely unweighted random index maps
|
|
idx = vector::Zeros(size);
|
|
for(ulong i = 0; i<shape; idx[i] = (double)CHighQualityRand::HQRndUniformI(rstate.GetInnerObj(),int(pop_size)), ++i);
|
|
}
|
|
}
|
|
// CASE 2: Sampling WITHOUT replacement (unique values forced)
|
|
else
|
|
{
|
|
if(size>pop_size)
|
|
{
|
|
Print(__FUNCTION__, " cannot take a larger sample than 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 _new,temp;
|
|
long uindices[],counted[];
|
|
|
|
// Loop structure selection blocks until we hit unique requested limit records
|
|
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);
|
|
|
|
// Mask out and invalidate previously selected entries by wiping their weight variables
|
|
if(n_uniq>0)
|
|
{
|
|
for(ulong i = 0; i<n_uniq; ++i)
|
|
pp[ulong(found[i])] = 0.0;
|
|
}
|
|
|
|
// Refresh normalized CDF tracks dynamically over available subsets
|
|
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
|
|
{
|
|
// Shuffling shortcut configuration logic without weighting distribution fields
|
|
idx = permutation(pop_size);
|
|
idx = sliceVector(idx,0,long(size));
|
|
idx.Resize(shape);
|
|
}
|
|
}
|
|
|
|
// Pull final vector configurations out matching parsed target positional indices mapping flags
|
|
return select(a,idx);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sample Vector (With Overloaded Replace Controller) |
|
|
//| Handles vector selection utilizing custom flags conditions flags |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool sampleVector(vector<T> &in,const ulong count,const bool replace,vector<T> &result)
|
|
{
|
|
// Shortcut directly back onto default handler if replacement paths are chosen
|
|
if(replace)
|
|
return sampleVector(in,count,result);
|
|
|
|
if(!in.Size())
|
|
return(false);
|
|
|
|
// Resize target objects validation check structures
|
|
if(result.Size()<count)
|
|
if(!result.Resize(count))
|
|
return(false);
|
|
|
|
ulong indices[];
|
|
|
|
// Generate continuous array maps mirroring data configurations dimensions indices
|
|
if(!arange(indices,int(in.Size())))
|
|
return(false);
|
|
|
|
// Fisher-Yates Shuffling Variant Logic pipeline setup processing block
|
|
for(ulong i=0; i<count; ++i)
|
|
{
|
|
// Select random position from elements remaining inside active target regions
|
|
ulong j=i+ulong((ulong)(in.Size()-i)*MathRand()/32768);
|
|
|
|
// Perform position element swap track steps directly
|
|
if(j!=i)
|
|
{
|
|
ulong t=indices[i];
|
|
indices[i]=indices[j];
|
|
indices[j]=t;
|
|
}
|
|
}
|
|
|
|
// Extract matching values using randomized uniquely shuffled target mapping lists
|
|
for(ulong i=0; i<count; ++i)
|
|
result[i]=in[indices[i]];
|
|
|
|
return(true);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sample Array |
|
|
//| Selects values uniformly from primitive arrays with replacement |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool sampleArray(T &array[],const ulong count,T &result[])
|
|
{
|
|
if(!array.Size())
|
|
return false;
|
|
|
|
// Allocation scaling checks for array buffers structures types parameters
|
|
if(result.Size()<uint(count))
|
|
if(ArrayResize(result,int(count))!=int(count))
|
|
return(false);
|
|
|
|
// Fill positions safely with pseudo-random math elements index maps
|
|
for(ulong i=0; i<count; ++i)
|
|
{
|
|
int ind=(int)(array.Size()*MathRand()/32768);
|
|
result[i]=array[ind];
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sample Array (With Replacement Controller Overload) |
|
|
//| Handles unique selection options on standard primitive variables |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool sampleArray(T &array[],const ulong count,const bool replace,T &result[])
|
|
{
|
|
// Route over tracking patterns back directly to baseline execution states
|
|
if(replace)
|
|
return sampleArray(array,count,result);
|
|
|
|
if(!array.Size())
|
|
return(false);
|
|
|
|
// Resize system buffer tracks sizes validations checks controls
|
|
if(result.Size()<uint(count))
|
|
if(ArrayResize(result,int(count))!=int(count))
|
|
return(false);
|
|
|
|
ulong indices[];
|
|
|
|
// Construct native indexes coordinates tracks lists configurations
|
|
if(!arange(indices,int(array.Size())))
|
|
return(false);
|
|
|
|
// Perform primitive layout array content structure randomized shuffle mapping operations
|
|
for(ulong i=0; i<count; ++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;
|
|
}
|
|
}
|
|
|
|
// Write metrics collections onto destination tracking parameters indices definitions arrays
|
|
for(ulong i=0; i<count; ++i)
|
|
result[i]=array[indices[i]];
|
|
|
|
return(true);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| EmpiricalProbabilityDensityEmpirical |
|
|
//| Computes the Empirical Probability Density Function (EPDF). |
|
|
//| divides the range of data into 'count' bins and normalizes. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool epdf(const vector<T> &vec, const ulong count, vector<T> &x, vector<T> &pdf)
|
|
{
|
|
// Bins must be greater than 1 to establish intervals
|
|
if(count<=1)
|
|
return(false);
|
|
|
|
ulong size=vec.Size();
|
|
if(size==0)
|
|
return(false);
|
|
|
|
//--- check NaN values
|
|
// Verify all elements are valid numbers before proceeding
|
|
for(ulong i=0; i<size; ++i)
|
|
{
|
|
if(!MathIsValidNumber(vec[i]))
|
|
return(false);
|
|
}
|
|
|
|
//--- prepare output arrays
|
|
// Ensure output arrays have enough capacity to hold the data
|
|
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) // Avoid division by zero if all elements are identical
|
|
return(false);
|
|
|
|
//--- calculate probability density of the empirical distribution
|
|
// Initialize intervals (x-axis coordinates for each bin)
|
|
for(ulong i=0; i<count; ++i)
|
|
{
|
|
x[i]=minv+i*range/(count-1);
|
|
pdf[i]=0;
|
|
}
|
|
|
|
// Distribute data points into their respective bins (histogram counting)
|
|
for(ulong i=0; i<size; ++i)
|
|
{
|
|
double v=(vec[i]-minv)/range;
|
|
ulong ind=ulong((v*(count-1)));
|
|
++pdf[ind];
|
|
}
|
|
|
|
//--- normalize values
|
|
// Normalize the frequency values so the integral under the curve equals 1.0
|
|
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 |
|
|
//| Computes the Empirical Cumulative Distribution Function (ECDF). |
|
|
//| Integrates the generated EPDF to calculate cumulative values. |
|
|
//+------------------------------------------------------------------+
|
|
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
|
|
// Step 1: Generate an internal PDF first
|
|
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
|
|
// Step 2: Accumulate the area under the PDF curve for the CDF
|
|
sum=0.0;
|
|
for(ulong i=0; i<count; ++i)
|
|
{
|
|
sum+=pdf[i]*dx;
|
|
cdf[i]=sum;
|
|
}
|
|
|
|
return(true);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| MathFactorial |
|
|
//| Calculates the factorial of an integer 'n'. |
|
|
//| Uses a lookup table for n <= 20 to improve performance. |
|
|
//+------------------------------------------------------------------+
|
|
double factorial(const uint n)
|
|
{
|
|
// Precomputed factorials from 0! to 20! to save memory/processing cycles
|
|
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 |
|
|
//| Computes Spearman's rank correlation coefficient between |
|
|
//| two data vectors. |
|
|
//+------------------------------------------------------------------+
|
|
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
|
|
// Converts raw values into rank orders to perform non-parametric correlation
|
|
vector rank_x = rankVector(vec1);
|
|
vector rank_y = rankVector(vec2);
|
|
|
|
//---
|
|
// Spearman correlation is equivalent to Pearson correlation on ranked variables
|
|
return rank_x.CorrCoef(rank_y);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| MathCorrelationKendall |
|
|
//| Computes Kendall's Tau-b rank correlation coefficient |
|
|
//| measuring concordant and discordant pairs. |
|
|
//+------------------------------------------------------------------+
|
|
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);
|
|
}
|
|
|
|
// cnt1/cnt2 track non-tied elements in each vector, cnt tracks (concordant - discordant)
|
|
long cnt1=0, cnt2=0, cnt=0;
|
|
|
|
//--- Compare every possible pair combination (i, j)
|
|
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) // At least one tie exists in the pair
|
|
{
|
|
if(delta1!=0)
|
|
++cnt1;
|
|
if(delta2!=0)
|
|
++cnt2;
|
|
}
|
|
else // No ties in either vector
|
|
{
|
|
++cnt1;
|
|
++cnt2;
|
|
if(delta>0.0)
|
|
++cnt; // Concordant pair
|
|
else
|
|
cnt--; // Discordant pair
|
|
}
|
|
}
|
|
}
|
|
|
|
//--- calculate Kendall tau
|
|
long den=cnt1*cnt2;
|
|
if(den==0)
|
|
{
|
|
Print(__FUNCTION__, " failed zero check at line 3139 ");
|
|
return tau;
|
|
}
|
|
|
|
// Calculate Tau formula adjusting for ties
|
|
tau=double(cnt)/MathSqrt(den);
|
|
|
|
return(tau);
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| MathRank |
|
|
//| Ranks the elements of a vector, assigning average fractional |
|
|
//| ranks to tied elements. Uses an in-place Heap Sort routine. |
|
|
//+------------------------------------------------------------------+
|
|
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[];
|
|
|
|
// Create an index sequence tracking original positions [0, 1, 2... size-1]
|
|
if(!arange(indices, int(size)))
|
|
return rank;
|
|
|
|
ulong i, j, k, t, tmpi;
|
|
double tmp;
|
|
|
|
//--- sort (Heap Sort Implementation keeping track of data indices)
|
|
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 values and tracking indices
|
|
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 root with current boundary element
|
|
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
|
|
// Finds contiguous repeating values and updates them to their fractional rank average
|
|
i=0;
|
|
while(i<size)
|
|
{
|
|
j=i+1;
|
|
while(j<size)
|
|
{
|
|
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
|
|
// Maps the computed ranks back into their original spatial arrangement order
|
|
for(i=0; i<size; ++i)
|
|
rank[indices[i]]=values[i];
|
|
|
|
return rank;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Modifies vector in place: elements matching 'value' (within |
|
|
//| maxdiff tolerance) become 1.0, others become 0.0. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsEqual(vector<T>&in, const T value, const T maxdiff)
|
|
{
|
|
in -= value;
|
|
whereIsLessEqual(fabs(in), fabs(maxdiff));
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Modifies matrix in place: elements matching 'value' (within |
|
|
//| maxdiff tolerance) become 1.0, others become 0.0. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsEqual(matrix<T>&in, const T value, const T maxdiff)
|
|
{
|
|
in -= value;
|
|
whereIsLessEqual(fabs(in), fabs(maxdiff));
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets vector elements to 1.0 if less than value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsLess(vector<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Size(); ++i)
|
|
in[i] = in[i]<value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets vector elements to 1.0 if greater than value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsGreater(vector<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Size(); ++i)
|
|
in[i] = in[i]>value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets vector elements to 1.0 if less or equal to value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsLessEqual(vector<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Size(); ++i)
|
|
in[i] = in[i]<=value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets vector elements to 1.0 if greater or equal to value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsGreaterEqual(vector<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Size(); ++i)
|
|
in[i] = in[i]>=value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets matrix elements to 1.0 if less than value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsLess(matrix<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Rows(); ++i)
|
|
for(ulong j = 0; j<in.Cols(); ++j)
|
|
in[i,j] = in[i,j]<value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets matrix elements to 1.0 if less or equal to value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsLessEqual(matrix<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Rows(); ++i)
|
|
for(ulong j = 0; j<in.Cols(); ++j)
|
|
in[i,j] = in[i,j]<=value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets matrix elements to 1.0 if greater than value, else 0.0 |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsGreater(matrix<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Rows(); ++i)
|
|
for(ulong j = 0; j<in.Cols(); ++j)
|
|
in[i,j] = in[i,j]>value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Sets matrix elements to 1.0 if greater or equal to value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
void whereIsGreaterEqual(matrix<T>&in, const T value)
|
|
{
|
|
for(ulong i = 0; i<in.Rows(); ++i)
|
|
for(ulong j = 0; j<in.Cols(); ++j)
|
|
in[i,j] = in[i,j]>=value?1.0:0.0;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask vector: 1.0 if element < value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask vector: 1.0 if element <= value, else 0|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask vector: 1.0 if element > value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask vector: 1.0 if element >= value, else 0|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask matrix: 1.0 if element < value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask matrix: 1.0 if element <= value, else 0|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask matrix: 1.0 if element > value, else 0.|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a new binary mask matrix: 1.0 if element >= value, else 0|
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a vector mask: 1.0 if rounded values equal target value. |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a matrix mask based on rounding precision. |
|
|
//| FIX?: Note that condition checks >= value instead of == value. |
|
|
//+------------------------------------------------------------------+
|
|
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) // Double check if this should be ==
|
|
out[i][j] = 1.0;
|
|
}
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a copy of the vector with NaN / irregular items replaced |
|
|
//| by a fallback filler value. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector<T> replace_whereVectorIsNaN(vector<T>&in, const T fill = 0.0)
|
|
{
|
|
vector<T>out = in;
|
|
|
|
// Quick exit check if the entire vector is completely made of NaN values
|
|
if(in.HasNan()==in.Size())
|
|
{
|
|
out.Fill(fill);
|
|
return out;
|
|
}
|
|
|
|
for(ulong i = 0; i<in.Size(); ++i)
|
|
{
|
|
// Catch NaN, Subnormal, or Inf configurations
|
|
if(MathClassify(in[i])!=FP_NORMAL)
|
|
out[i] = fill;
|
|
}
|
|
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Returns a copy of the matrix with NaN components replaced |
|
|
//| by a fallback filler value. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> replace_whereMatrixIsNaN(matrix<T>&in, const T fill = 0.0)
|
|
{
|
|
matrix<T>out = in;
|
|
|
|
// Quick exit check if every single matrix location is NaN
|
|
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;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|
|
//| Selective column copier: Transfer columns from 'from' matrix |
|
|
//| into 'copyto' matrix ONLY if v_mask element is non-zero (true). |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool copyToSelectCols(matrix<T> ©to, matrix<T> &from, vector &v_mask)
|
|
{
|
|
// Bounds verification checking that shapes/dimensions line up correctly
|
|
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 mask evaluated index is active, update targeted column
|
|
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 |
|
|
//| Copies rows from 'from' matrix to 'copyto' matrix where the |
|
|
//| corresponding index in 'v_mask' evaluates to true. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool copyToSelectRows(matrix<T> ©to, matrix<T> from, vector &v_mask)
|
|
{
|
|
// Validate dimension compatibility: source, destination, and mask must match row-wise
|
|
if(from.Rows()!=copyto.Rows() || v_mask.Size()!=copyto.Rows())
|
|
{
|
|
Print(__FUNCTION__, " invalid parameters ");
|
|
return false;
|
|
}
|
|
|
|
// Iterate through each row
|
|
for(ulong i = 0; i<copyto.Rows(); ++i)
|
|
{
|
|
// If the mask element is non-zero, perform the row copy
|
|
if(v_mask[i])
|
|
{
|
|
if(!copyto.Row(from.Row(i),i))
|
|
{
|
|
Print(__FUNCTION__, " failed column assignment ", GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Copy select elements |
|
|
//| Copies individual elements from the 'from' vector to 'copyto' |
|
|
//| vector using 'v_mask' as an element-wise conditional switch. |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
bool copySelectElements(vector<T> ©to, vector<T> &from, vector &v_mask)
|
|
{
|
|
// Ensure all vectors share the same length
|
|
if(from.Size()!=copyto.Size() || v_mask.Size()!=copyto.Size())
|
|
{
|
|
Print(__FUNCTION__, " invalid parameters ");
|
|
return false;
|
|
}
|
|
|
|
// Transfer elements conditionally based on mask values
|
|
for(ulong i = 0; i<copyto.Size(); ++i)
|
|
{
|
|
if(v_mask[i])
|
|
copyto[i] = from[i];
|
|
}
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Sign function |
|
|
//| Returns -1 for negative values, 1 for positive, and 0 for zero. |
|
|
//+------------------------------------------------------------------+
|
|
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 |
|
|
//| Computes statistical moments efficiently using a binary |
|
|
//| exponentiation breakdown to calculate power transformations. |
|
|
//+------------------------------------------------------------------+
|
|
double moment(vector& a, ulong order, double mean=NULL)
|
|
{
|
|
// Edge case: Empty vector defaults to statistical mean
|
|
if(!a.Size())
|
|
return a.Mean();
|
|
|
|
// Direct results for 0th or 1st-order basic moments
|
|
if(order == 0 || (order == 1 && !mean))
|
|
return order==0?1.0:0.;
|
|
|
|
// Build an execution chain for dynamic power processing (binary exponentiation)
|
|
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;
|
|
}
|
|
|
|
// Identify the distribution center (mean)
|
|
double mn = (mean == 0.0)?a.Mean():mean;
|
|
|
|
// Center the data values
|
|
vector azeromean = a - mn;
|
|
vector s;
|
|
if(list[list.Size()-1] == 1)
|
|
s = azeromean;
|
|
else
|
|
s = pow(azeromean,2.);
|
|
|
|
// Unwind the exponentiation tree to compute final elevated matrix/vector power states
|
|
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) implementation |
|
|
//| Computes sin(pi*x) reliably by wrapping the argument inside the |
|
|
//| standard period range [0, 2) to eliminate large float errors. |
|
|
//+------------------------------------------------------------------+
|
|
template <typename T>
|
|
T sinpi(T x)
|
|
{
|
|
T s = 1.0;
|
|
|
|
// Force input positivity and track original quadrant sign
|
|
if(x < 0.0)
|
|
{
|
|
x = -x;
|
|
s = -1.0;
|
|
}
|
|
|
|
// Map value into a standard [0.0, 2.0) periodicity window
|
|
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));
|
|
}
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Log Gamma asymptotic expansion for large values of x |
|
|
//| Uses Stirling's approximation to prevent overflow on large input |
|
|
//+------------------------------------------------------------------+
|
|
double lgam_large_x(double x)
|
|
{
|
|
const double LS2PI = 0.91893853320467274178; // Logarithm of sqrt(2*PI)
|
|
double q = (x - 0.5) * log(x) - x + LS2PI;
|
|
|
|
// If x is massively large, additional fraction expansions drop below precision limits
|
|
if(x > 1.0e8)
|
|
{
|
|
return (q);
|
|
}
|
|
// Polynomial correction factor optimization
|
|
double p = 1.0 / (x * x);
|
|
p = ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x;
|
|
return q + p;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Horner's Method Polynomial Evaluation |
|
|
//| Evaluates standard algebraic polynomials given standard coefficients|
|
|
//+------------------------------------------------------------------+
|
|
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);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Horner's Method Variant for Monic Polynomials |
|
|
//| Polynomial evaluation with assumptions that the highest order |
|
|
//| coefficient equals 1. |
|
|
//+------------------------------------------------------------------+
|
|
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-aware Log Gamma Calculation |
|
|
//| Computes ln(|Gamma(x)|) and returns the sign of Gamma(x) via ref.|
|
|
//+------------------------------------------------------------------+
|
|
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; // Guard rails for double numeric limitations
|
|
|
|
// Coefficients for minimax approximations over specific scalar domains
|
|
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;
|
|
|
|
// Reject NaNs or Infinite structures immediately
|
|
if(MathClassify(x)!=FP_NORMAL)
|
|
{
|
|
return x;
|
|
}
|
|
|
|
// Use Reflection Formula transformation principles for highly negative inputs
|
|
if(x < -34.0)
|
|
{
|
|
q = -x;
|
|
w = lgam_sgn(q, sign);
|
|
p = floor(q);
|
|
if(p == q)
|
|
{
|
|
Print(__FUNCTION__, " singular error ");
|
|
return (double("inf")); // Gamma undefined for negative integers
|
|
}
|
|
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);
|
|
}
|
|
|
|
// For small to mid ranges, recurse and rely on rational evaluations
|
|
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);
|
|
}
|
|
|
|
// Handle large inputs via Stirling scaling variations
|
|
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) |
|
|
//| Prevents serious cancellation precision losses when parameters differ |
|
|
//| significantly by magnitude scale. |
|
|
//+-------------------------------------------------------------------------+
|
|
double lbeta_asymp(double a, double b, int &sgn)
|
|
{
|
|
double r = lgam_sgn(b, sgn);
|
|
r -= b * log(a);
|
|
|
|
// Infinite series correction structures
|
|
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 Beta handling for negative integer inputs |
|
|
//+------------------------------------------------------------------+
|
|
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");
|
|
}
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Special case Log Beta handling for negative integers |
|
|
//+------------------------------------------------------------------+
|
|
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");
|
|
}
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Log Beta Function: ln(B(a,b)) |
|
|
//| Handles reflection limits, variable sorting, and asymptotic scaling|
|
|
//+------------------------------------------------------------------+
|
|
double lbeta(double a, double b)
|
|
{
|
|
double y;
|
|
int sign;
|
|
const double MAXGAM = 171.624376956302725; // Standard double precision cutoff limit for basic gamma
|
|
const double beta_ASYMP_FACTOR = 1e6;
|
|
sign = 1;
|
|
|
|
// Process negative or integer domain transitions for 'a'
|
|
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"));
|
|
}
|
|
}
|
|
}
|
|
|
|
// Process negative or integer domain transitions for 'b'
|
|
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"));
|
|
}
|
|
}
|
|
}
|
|
|
|
// Enforce optimization symmetry mapping (a >= b)
|
|
if(MathAbs(a) < MathAbs(b))
|
|
{
|
|
y = a;
|
|
a = b;
|
|
b = y;
|
|
}
|
|
|
|
// Deploy specialized fast asymptotic formulas if a is dramatically larger than b
|
|
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;
|
|
}
|
|
|
|
// Use Log Gamma components directly if standard values bypass safe range thresholds
|
|
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);
|
|
}
|
|
|
|
// Direct native calculation if arguments fit inside reliable evaluation limits
|
|
y = CGammaFunc::GammaFunc(y);
|
|
a = CGammaFunc::GammaFunc(a);
|
|
b = CGammaFunc::GammaFunc(b);
|
|
if(y == 0.0)
|
|
{
|
|
Print(__FUNCTION__, " overflow error ");
|
|
return (sign * double("inf"));
|
|
}
|
|
|
|
// Dynamic ratio factorization balancing to maximize float precision accuracy
|
|
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));
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Generalized Binomial Coefficient calculation |
|
|
//| Computes combinations allowing real parameters (n choose k) |
|
|
//+------------------------------------------------------------------+
|
|
double binom(double n, double k)
|
|
{
|
|
double kx,nx,num, den, dk, sgn;
|
|
// Undefined conditions validation
|
|
if(n < 0)
|
|
{
|
|
nx = floor(n);
|
|
if(n == nx)
|
|
{
|
|
return double("nan"); // Undefined for negative integer 'n'
|
|
}
|
|
}
|
|
|
|
kx = floor(k);
|
|
// Optimization layer targeting standard small integer fields
|
|
if(k == kx && (MathAbs(n) > 1e-8 || n == 0))
|
|
{
|
|
nx = floor(n);
|
|
if(nx == n && kx > nx / 2 && nx > 0)
|
|
{
|
|
// Reduce kx using identity properties of symmetry
|
|
kx = nx - kx;
|
|
}
|
|
|
|
// Explicit iterative loops for low ranges
|
|
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;
|
|
// Balance scale to check overflow mid-loop
|
|
if(MathAbs(num) > 1e50)
|
|
{
|
|
num /= den;
|
|
den = 1.0;
|
|
}
|
|
}
|
|
return num / den;
|
|
}
|
|
}
|
|
|
|
// General Log-Beta path logic for large scales to minimize overflow
|
|
if(n >= 1e10 * k && k > 0)
|
|
{
|
|
return exp(-lbeta(1 + n - k, 1 + k) - log(n + 1));
|
|
}
|
|
|
|
// Approximation mechanisms for extremely large indices using Sine reflection properties
|
|
if(k > 1e8 * MathAbs(n))
|
|
{
|
|
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);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Combination Count interface: N choose K |
|
|
//| Supports variations containing exact matching or repetition paths|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
T comb(ulong N, T k, bool exact=false, bool repetition = false)
|
|
{
|
|
// Handle combinations with replacements allowed
|
|
if(repetition)
|
|
return comb(N+ulong(k)-1,(T)k,exact);
|
|
if(exact)
|
|
return comb(N,k);
|
|
else
|
|
{
|
|
// Boundary error protection evaluation
|
|
bool cond = false;
|
|
cond = ((k<=T(N)) && (N>=0) && (k>=0));
|
|
T val = (T)binom(double(N),double(k));
|
|
if(!cond)
|
|
val = 0; // Return zero for mathematically invalid constraints
|
|
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); // 2nd Moment (Variance)
|
|
double m4 = moment(in,4,mean); // 4th Moment
|
|
|
|
// Guard logic against zero-variance inputs
|
|
bool zero = (m2<= pow(mean*2.e-16,2.));
|
|
out = (!zero)?m4/pow(m2,2.):double("nan");
|
|
|
|
// Strip biases for small sample profiles if flagged
|
|
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");
|
|
}
|
|
|
|
// Convert from Pearson default format to Fisher excess metrics if flagged
|
|
if(fisher)
|
|
out = out - 3.;
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Get Element-wise Vector Maximum |
|
|
//| Returns a vector containing the largest item at matching steps |
|
|
//+------------------------------------------------------------------+
|
|
vector maximum(vector &x, vector &y)
|
|
{
|
|
if(x.Size()!=y.Size())
|
|
{
|
|
Print(__FUNCTION__, " vector size mismatch ");
|
|
return vector::Zeros(0);
|
|
}
|
|
|
|
// Bundle both inputs horizontally to evaluate row-wise metrics
|
|
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); // Column axis evaluation
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Get Element-wise Vector Minimum |
|
|
//| Returns a vector containing the smallest item at matching steps |
|
|
//+------------------------------------------------------------------+
|
|
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 Vector/Scalar 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 Vector/Scalar 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 mask highlighting matching NaN conditions |
|
|
//+------------------------------------------------------------------+
|
|
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++] = 1.0;
|
|
return out;
|
|
}
|
|
else
|
|
return vector::Zeros(0);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Return mask tracking locations without NaN objects |
|
|
//+------------------------------------------------------------------+
|
|
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++] = 1.0;
|
|
return out;
|
|
}
|
|
else
|
|
return vector::Ones(in.Size());
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Matrix Rank Determination |
|
|
//| Computes rank by identifying structural Singular Value |
|
|
//| Decomposition counts exceeding tolerance calculations. |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
}
|
|
|
|
// Compute a machine epsilon threshold tolerance scale relative to maximum dimensions
|
|
double rtol = MathMax(A.Rows(),A.Cols())*DBL_EPSILON;
|
|
double tol = s.Max()*rtol;
|
|
|
|
// Filter values above computational noise bands
|
|
vector count = whereVectorIsGt(s,tol);
|
|
return ulong(count.Sum());
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Toeplitz Matrix Construction |
|
|
//| Builds symmetric diagonal banded matrices where each row is shifted|
|
|
//| one step relative to the row above it. |
|
|
//+------------------------------------------------------------------+
|
|
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]; // Upper triangle assignments
|
|
else
|
|
out[i,j] = x[i-j]; // Lower triangle assignments
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| Compute the (Moore-Penrose) pseudo-inverse of a matrix |
|
|
//| Implements SVD matrix conversions to generate stable inversions |
|
|
//| for non-square or ill-conditioned matrices. |
|
|
//+------------------------------------------------------------------+
|
|
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);
|
|
}
|
|
|
|
// Fallback handling to explicitly ensure calculation safety
|
|
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();
|
|
}
|
|
|
|
// Invert valid non-zero singular coefficients
|
|
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;
|
|
|
|
// Reconstruct inversion states via matrix composition maps (V * S_inv * U^T)
|
|
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 a vector |
|
|
//| Implements a random index permutation using the Fisher-Yates method|
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
vector shuffleVector(vector<T> &in,CHighQualityRandStateShell* rstate=NULL)
|
|
{
|
|
vector<T>out(in.Size());
|
|
bool destroy = false;
|
|
|
|
// Instantiate local random state if caller did not supply one
|
|
if(rstate == NULL)
|
|
{
|
|
destroy = true;
|
|
rstate = new CHighQualityRandStateShell();
|
|
CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
|
|
}
|
|
|
|
ulong k=0;
|
|
ulong size=in.Size();
|
|
out=in;
|
|
T tmp;
|
|
|
|
// Perform element position swaps across uniform random distributions
|
|
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 an array |
|
|
//+------------------------------------------------------------------+
|
|
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();
|
|
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 matrix components (either row-wise or column-wise) |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Process and re-inject vectors individually
|
|
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 |
|
|
//| Storage container encapsulating training and validation testing |
|
|
//| index splits for Cross-Validation routines. |
|
|
//+------------------------------------------------------------------+
|
|
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 split mapping function |
|
|
//| Divides dataset arrays into training and testing partitions. |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Scramble the base sequence indices before splitting if requested
|
|
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));
|
|
|
|
// Distribute remainder elements evenly across initial chunks
|
|
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;
|
|
|
|
// Segment dataset ranges into consecutive testing blocks and aggregate remainders as training data
|
|
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 Splitter (Time Series Split) |
|
|
//| Generates forward-chaining rolling folds where historical spaces |
|
|
//| accumulate over time to predict subsequent intervals. |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Everything leading up to the evaluation slice becomes training context
|
|
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 (Banana Function) |
|
|
//| Popular non-convex mathematical optimization test 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 |
|
|
//| Calculates the first-order partial derivatives vector array. |
|
|
//+------------------------------------------------------------------+
|
|
vector rosen_gradient(vector &x)
|
|
{
|
|
vector grad = vector::Zeros(x.Size());
|
|
if(x.Size()>2)
|
|
{
|
|
// Slice vectors to compute partial values efficiently in blocks
|
|
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);
|
|
}
|
|
// Explicit derivative boundary conditions for head and tail variables
|
|
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 |
|
|
//| Constructs the square matrix of second-order partial derivatives. |
|
|
//+------------------------------------------------------------------+
|
|
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;
|
|
|
|
// Set up off-diagonal element bands
|
|
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;
|
|
|
|
// Calculate primary diagonal values
|
|
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);
|
|
}
|
|
|
|
// Combine off-diagonal and main diagonal matrices to complete the Hessian
|
|
return H+h3;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|Convert vector into either a column or row matrix |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> asMatrix(vector<T>& in , bool as_column_matrix = true)
|
|
{
|
|
// Prepare output matrix
|
|
matrix<T> mat = as_column_matrix?matrix<T>::Zeros(in.Size(),1):matrix<T>::Zeros(1,in.Size());
|
|
// Add vector to new matrix
|
|
if((as_column_matrix && !mat.Col(in,0)) ||
|
|
(!as_column_matrix && !mat.Row(in,0)))
|
|
{
|
|
Print(__FUNCTION__," Row insertion failure ", GetLastError());
|
|
return matrix<T>::Zeros(0,0);
|
|
}
|
|
// Return matrix
|
|
return mat;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|Mahalanobis distance calculation |
|
|
//+------------------------------------------------------------------+
|
|
double mahalanobis(vector& u,vector& v, matrix& VI)
|
|
{
|
|
//---
|
|
vector delta = u - v;
|
|
//---
|
|
return sqrt((delta.MatMul(VI)).Dot(delta));
|
|
//---
|
|
}
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|