//+------------------------------------------------------------------+
//|                                                           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> &copyto,matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
   return matrixCopy(copyto,copyfrom,BEGIN,END,STEP,dest_start,dest_stop,dest_step);
  }

//+------------------------------------------------------------------+
//|  assign values to a portion of matrix selected within range rows |
//+------------------------------------------------------------------+
template <typename T>
bool matrixCopyRows(matrix<T> &copyto, matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
   return matrixCopy(copyto,copyfrom,dest_start,dest_stop,dest_step);
  }

//+------------------------------------------------------------------+
//|return vector of arbitrary elements (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> &copyto, vector<T> &copyfrom, 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> &copyto, 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> &copyto, matrix &copyfrom, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP, const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
  {
//--- 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> &copyto, T fill_value, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP, const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
  {
//--- 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 &copy_from[], matrix &copy_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> &copyto, 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> &copyto, 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> &copyto, 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));
//---
  }
}

//+------------------------------------------------------------------+