//+------------------------------------------------------------------+
//|                                                           np.mqh |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Math\Stat\Normal.mqh>
#include<Math\Alglib\specialfunctions.mqh>
#include<Math\Alglib\linalg.mqh>
#define _NP_
#define BEGIN 0
#define STEP 1
#define END LONG_MAX
#define BEGIN_REVERSE -1
#define STEP_REVERSE -1
#define END_REVERSE LONG_MIN
//+------------------------------------------------------------------+
//| library constants                                                |
//+------------------------------------------------------------------+
namespace np
{
namespace internal
{

//+------------------------------------------------------------------+
//|find index of a sorted array such that arr[i] <= key < arr[i + 1] |
//+------------------------------------------------------------------+
int binary_search(const double key, const double &arr[])
  {

   int imin = 0;
   int imax = ArraySize(arr)-1;
   int mid = -1;

   /* Handle keys outside of the arr range first */
   if(key > arr[arr.Size()-1])
     {
      return ArraySize(arr);
     }
   else
      if(key < arr[0])
        {
         return -1;
        }

   while(imin<=imax)
     {
      mid = (imin+imax)/2;
      if(mid<int(arr.Size()-1) && arr[mid]<=key && arr[mid+1]>key)
         return mid;
      else
        {
         if(key<arr[mid])
            imax = mid - 1;
         else
            imin = mid + 1;
        }
     }
   return ArraySize(arr) - 2;


  }
//+------------------------------------------------------------------+
//| cast string data to vector format doubles and dates              |
//+------------------------------------------------------------------+
vector formatCol(string &Arr[],bool handle_dates,int date_column, int column_index)
  {
   int size = ArraySize(Arr);
   vector ret(size);
//---
   if(handle_dates && date_column == column_index)
     {
      for(int i=0; i<size; ++i)
         ret[i] = (double)StringToTime(Arr[i]);
     }
   else
     {
      for(int i=0; i<size; ++i)
         ret[i] = StringToDouble(Arr[i]);
     }
//---
   return ret;
  }
//+------------------------------------------------------------------+
//|get column                                                        |
//+------------------------------------------------------------------+
template<typename T>
void col(const T &Matrix[], T &Col[], int column, int cols)
  {
   int rows = ArraySize(Matrix)/cols;
   ArrayResize(Col,rows);

   int start = 0;
   for(int i=0; i<cols; ++i)
     {
      start = i;

      if(i != column-1)
         continue;
      else
         for(int j=0; j<rows; ++j)
           {
            //printf("ColMatrix[%d} Matrix{%d]",j,start);
            Col[j] = Matrix[start];

            start += cols;
           }
     }
  }
//+------------------------------------------------------------------+
//| product for arrays                                               |
//+------------------------------------------------------------------+
template<typename T>
T prodArray(T &in_array[], uint from=0)
  {
   T result = T(1);
   for(uint i = from; i<in_array.Size(); ++i)
      result*=in_array[i];
   return result;

  }
//+------------------------------------------------------------------+
//|calculate size of container                                       |
//+------------------------------------------------------------------+
long calculatesize(long initial,long stop,long step)
  {
//---calculate len
   long len = fabs(initial-stop)/fabs(step);
//---
   if(len*fabs(step) < fabs(initial-stop))
      len += (fabs(initial-stop) - len*fabs(step))/fabs(step);
//---
   return len;
  }
//+------------------------------------------------------------------+
//| process start parameter                                          |
//+------------------------------------------------------------------+
long normalizestart(long start,ulong size)
  {
   long s = (start<0 && ulong(fabs(start))<size)?long(size)+start:(start<0 && ulong(fabs(start))>=size)?0:start;
   return s;
  }
//+------------------------------------------------------------------+
//|process stop parameter                                            |
//+------------------------------------------------------------------+
long normalizestop(long stop, ulong size)
  {
   long s = (stop>0 && ulong(stop)>=size)?long(size):(stop<0 && ulong(fabs(stop))>size)?-1:(stop<0 && ulong(fabs(stop))<=size)?long(size)+stop:stop;
   return s;
  }
//+------------------------------------------------------------------+
//| validate all calculated params                                   |
//+------------------------------------------------------------------+
bool invalidparams(long start,long stop, long step, ulong size)
  {
   if(step==0 || start==stop || ulong(start)>=size || start<0 || ulong(fabs(stop))>size || (step<0 && start<stop) || (step>0 && start>stop))
      return true;
   return false;
  }
//+------------------------------------------------------------------+
//| quantiles helper function                                        |
//+------------------------------------------------------------------+
vector quantile(vector &in,vector &abprobs,vector &probs)
  {
   vector x = in;
   if(!sort(x))
     {
      Print(__FUNCTION__, " sort failure ");
      vector::Zeros(probs.Size());
     }

   ulong n = x.Size();

   if(!x.Size())
     {
      Print(__FUNCTION__, " invalid parameter, empty vector ");
      vector::Zeros(probs.Size());
     }

   vector aleph = (double(n)*probs+abprobs);

   vector k = aleph;
   if(!k.Clip(1.0, double(n-1)))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(probs.Size());
     }

   k = floor(k);
   vector gamma =  aleph-k;
   if(!gamma.Clip(0.0, 1.0))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(probs.Size());
     }

   vector k1 = k-1.0;
   ulong index[];
   if(!np::vecAsArray(k1,index))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(probs.Size());
     }

   vector x1 = np::select(k1,index);

   ArrayFree(index);

   if(!np::vecAsArray(k,index))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(probs.Size());
     }

   vector x2 = np::select(k,index);

   return (1.0-gamma)*x1+gamma*x2;

  }

}
//+------------------------------------------------------------------+
//| Return arrays representing indices of a grid                     |
//+------------------------------------------------------------------+
void indices(ulong rows, ulong cols, matrix &outRC[])
  {
   ulong N = 2;

   if(ArrayResize(outRC,2)<0 || !outRC[0].Resize(rows,cols) || !outRC[1].Resize(rows,cols))
     {
      Print(__FUNCTION__, " resize error ");
      return;
     }

   vector v = np::arange(rows);
   outRC[0] = np::repeat_vector_as_rows_cols(v,cols,false);
   v = np::arange(cols);
   outRC[1] = np::repeat_vector_as_rows_cols(v,rows,true);

  }
//+------------------------------------------------------------------+
//| The function sorts (in) vector and indices[] simultaneously      |
//| using QuickSort algorithm. After sorting the initial ordering    |
//| of the elements is located in the indices[] array.               |
//|                                                                  |
//| Arguments:                                                       |
//| in          : vector with values to sort                         |
//| ascending   : direction of sort (true:ascending)                 |
//| (false:descending)                                               |
//| first       : First element index                                |
//| last        : Last element index                                 |
//|                                                                  |
//| Return value: None                                               |
//+------------------------------------------------------------------+
template <typename T>
bool quickSort(vector<T> &in,bool ascending,long first,long last)
  {
   long    i,j;
   T p_double,t_double;
//--- check
   if(first<0 || last<0)
     {
      return false;
     }
//--- sort
   i=first;
   j=last;
   while(i<last)
     {
      //--- ">>1" is quick division by 2
      p_double=in[(first+last)>>1];
      while(i<j)
        {
         while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
           {
            //--- control the output of the in bounds
            if(i==in.Size()-1)
               break;
            i++;
           }
         while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
           {
            //--- control the output of the in bounds
            if(j==0)
               break;
            j--;
           }
         if(i<=j)
           {
            //-- swap elements i and j
            t_double=in[i];
            in[i]=in[j];
            in[j]=t_double;
            //--
            i++;
            //--- control the output of the in bounds
            if(j==0)
               break;
            j--;
           }
        }
      if(first<j)
         quickSort(in,ascending,first,j);
      first=i;
      j=last;
     }

   return true;
  }
//+------------------------------------------------------------------+
//| The function sorts (in) vector and indices[] simultaneously      |
//| using QuickSort algorithm. After sorting the initial ordering    |
//| of the elements is located in the indices[] array.               |
//|                                                                  |
//| Arguments:                                                       |
//| in          : vector with values to sort                         |
//| ascending   : direction of sort (true:ascending)                 |
//| (false:descending)                                               |
//| indices     : array of indexes                                   |
//| first       : First element index                                |
//| last        : Last element index                                 |
//|                                                                  |
//| Return value: None                                               |
//+------------------------------------------------------------------+
template <typename T>
bool quickSortIndices(vector<T>& _in,bool ascending,long &indices[],long first,long last)
  {
   vector<T> in = _in;
   long    i,j,t_int;
   T p_double,t_double;
//--- check
   if(first<0 || last<0)
     {
      return false;
     }
//--- sort
   i=first;
   j=last;
   while(i<last)
     {
      //--- ">>1" is quick division by 2
      p_double=in[(first+last)>>1];
      while(i<j)
        {
         while((ascending && in[i]<p_double) || (!ascending && in[i]>p_double))
           {
            //--- control the output of the in bounds
            if(i==in.Size()-1)
               break;
            i++;
           }
         while((ascending && in[j]>p_double) || (!ascending && in[j]<p_double))
           {
            //--- control the output of the in bounds
            if(j==0)
               break;
            j--;
           }
         if(i<=j)
           {
            //-- swap elements i and j
            t_double=in[i];
            in[i]=in[j];
            in[j]=t_double;
            //-- swap indices i and j
            t_int=indices[i];
            indices[i]=indices[j];
            indices[j]=t_int;
            i++;
            //--- control the output of the in bounds
            if(j==0)
               break;
            j--;
           }
        }
      if(first<j)
         quickSortIndices(in,ascending,indices,first,j);
      first=i;
      j=last;
     }
   return true;
  }

//+------------------------------------------------------------------+
//| sort a vector                                                    |
//+------------------------------------------------------------------+
template<typename T>
bool sort(vector<T>&in, bool ascending=true)
  {
   if(in.Size()<1)
      return false;

   return quickSort(in,ascending,0,in.Size()-1);
  }
//+------------------------------------------------------------------+
//| sort matrix by columns or rows                                   |
//+------------------------------------------------------------------+
template<typename T>
bool sort(matrix<T>&in, bool by_rows, bool ascending=true)
  {
   if(by_rows)
     {
      for(ulong i = 0; i<in.Rows(); ++i)
        {
         vector row = in.Row(i);
         if(!sort(row,ascending))
            return false;
         in.Row(row,i);
        }
     }
   else
     {
      for(ulong i = 0; i<in.Cols(); ++i)
        {
         vector col = in.Col(i);
         if(!sort(col,ascending))
            return false;
         in.Col(col,i);
        }
     }
   return true;
  }
//+------------------------------------------------------------------+
//|  extract a portion of vector within a range at specific interval |
//|  with support for negative indexing python style                 |
//| returns new vector from start to stop-1 at step intervals        |
//| start - index of first element (inclusive)                       |
//| stop  - index past last element (exclusive)                      |
//| step  - step size                                                |
//+------------------------------------------------------------------+
template<typename T>
vector<T> sliceVector(vector<T> &in, const long index_start=BEGIN, const long index_stop = END, const long index_step = STEP)
  {
//---local variables
   long st,sp,stp,size;
//---check function parameters
   st = internal::normalizestart(index_start,in.Size());
   sp = internal::normalizestop(index_stop,in.Size());
   stp = index_step;
//---error handling
   if(internal::invalidparams(st,sp,stp,in.Size()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return vector<T>::Zeros(0);
     }
//---calculate size of new vector
   size = internal::calculatesize(st,sp,stp);
//---initialize output vector
   vector<T> out(ulong(size));
//---assign elements to output
   for(long pos = st, index=0; index<size; ++index,pos+=stp)
      out[index] = in[pos];
//---done
   return out;
  }
//+------------------------------------------------------------------+
//| set subset of contiguous container elements to a certain value   |
//+------------------------------------------------------------------+
template<typename T>
bool fillVector(vector<T> &in, T value, const long index_start=BEGIN, const long index_stop = END, const long index_step = STEP)
  {
//---local variables
   long st,sp,stp,size;
//---check function parameters
   st = internal::normalizestart(index_start,in.Size());
   sp = internal::normalizestop(index_stop,in.Size());
   stp = index_step;
//---error handling
   if(internal::invalidparams(st,sp,stp,in.Size()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return false;
     }
//---calculate size of new vector
   size = internal::calculatesize(st,sp,stp);
//---assign elements to output
   for(long pos = st, index=0; index<size; ++index,pos+=stp)
      in[pos] = value;
//---done
   return true;
  }
//+------------------------------------------------------------------+
//| reverse a vector                                                 |
//+------------------------------------------------------------------+
template<typename T>
bool reverseVector(vector<T> &in)
  {
   if(in.Size()<2)
      return in.Size()!=0;

   ulong size = in.Size()-1;
   ulong half = (in.Size())>>1;
   for(ulong i = 0; i<half; ++i)
     {
      T temp = in[i];
      in[i] = in[size - i];
      in[size - i] = temp;
     }
   return true;
  }
//+------------------------------------------------------------------+
//|  extract a portion a matrix within a range at specific intervals |
//|  with support for negative indexing python style                 |
//| returns new matrix                                               |
//| row_start - row index of first row (inclusive)                   |
//| row_stop  - row index past last row (exclusive)                  |
//| row_step  - row step size                                        |
//| column_start - columnindex of first column (inclusive)           |
//| column_stop  - column index past last columun (exclusive)        |
//| column_step  - column step size                                  |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrix(matrix<T> &in, const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP,const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
  {
//---local variables
   long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
   rst =internal::normalizestart(row_start,in.Rows());
   rsp =internal::normalizestop(row_stop,in.Rows());
   rstp = row_step;
//---error handling
   if(internal::invalidparams(rst,rsp,rstp,in.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters for row specification ");
      return matrix<T>::Zeros(0,0);
     }
//---calculate size of new matrix rows
   rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
   cst = internal::normalizestart(column_start,in.Cols());
   csp = internal::normalizestop(column_stop,in.Cols());
   cstp = column_step;
//---error handling
   if(internal::invalidparams(cst,csp,cstp,in.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters for column specification ");
      return matrix<T>::Zeros(0,0);
     }
//---calculate size of new matrix columns
   csize = internal::calculatesize(cst,csp,cstp);
//---initialize output matrix
   matrix<T> out(ulong(rsize),ulong(csize));
//---assign elements to output
   for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
      for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
         out[row_index][column_index] = in[rpos][cpos];
//---done
   return out;
  }
//+------------------------------------------------------------------+
//| fill a portion a matrix within a range at specific intervals     |
//|  with support for negative indexing python style                 |
//| edits container in place                                         |
//| row_start - row index of first row (inclusive)                   |
//| row_stop  - row index past last row (exclusive)                  |
//| row_step  - row step size                                        |
//| column_start - columnindex of first column (inclusive)           |
//| column_stop  - column index past last columun (exclusive)        |
//| column_step  - column step size                                  |
//+------------------------------------------------------------------+
template<typename T>
bool fillMatrix(matrix<T> &in, const T value,const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP,const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
  {
//---local variables
   long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
   rst =internal::normalizestart(row_start,in.Rows());
   rsp =internal::normalizestop(row_stop,in.Rows());
   rstp = row_step;
//---error handling
   if(internal::invalidparams(rst,rsp,rstp,in.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters for row specification ");
      return false;
     }
//---calculate size of new matrix rows
   rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
   cst = internal::normalizestart(column_start,in.Cols());
   csp = internal::normalizestop(column_stop,in.Cols());
   cstp = column_step;
//---error handling
   if(internal::invalidparams(cst,csp,cstp,in.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters for column specification ");
      return false;
     }
//---calculate size of new matrix columns
   csize = internal::calculatesize(cst,csp,cstp);
//---assign elements to output
   for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
      for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
         in[rpos][cpos] = value;
//---done
   return true;
  }
//+------------------------------------------------------------------+
//|  extract a portion of matrix by selecting range of columns and   |
//|  optional interval                                               |
//|  with support for negative indexing python style                 |
//| returns new matrix                                               |
//| column_start - columnindex of first column (inclusive)           |
//| column_stop  - column index past last columun (exclusive)        |
//| column_step  - column step size                                  |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrixCols(matrix<T> &in, const long column_start=BEGIN, const long column_stop = END, const long column_step = STEP)
  {
   return sliceMatrix(in,BEGIN,END,STEP,column_start,column_stop,column_step);
  }
//+------------------------------------------------------------------+
//|reverse columns of matrix                                         |
//+------------------------------------------------------------------+
template<typename T>
bool reverseMatrixCols(matrix<T> &in)
  {

   ulong cols = in.Cols()-1;

   if(in.Cols()<2)
      return in.Cols()!=0;

   ulong cols_end = (in.Cols())>>1;

   for(ulong i = 0; i<cols_end; ++i)
      if(!in.SwapCols(i,cols-i))
         return false;

   return true;
  }
//+------------------------------------------------------------------+
//|  extract a portion of matrix by selecting range of rows and      |
//|  optional interval                                               |
//|  with support for negative indexing python style                 |
//| returns new matrix                                               |
//| row_start - row index of first row (inclusive)                   |
//| row_stop  - row index past last row (exclusive)                  |
//| row_step  - row step size                                        |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> sliceMatrixRows(matrix<T> &in, const long row_start=BEGIN, const long row_stop = END, const long row_step = STEP)
  {
   return sliceMatrix(in,row_start,row_stop,row_step);
  }
//+------------------------------------------------------------------+
//|reverse rows of matrix                                            |
//+------------------------------------------------------------------+
template<typename T>
bool reverseMatrixRows(matrix<T> &in)
  {
   ulong rows = in.Rows()-1;

   if(in.Rows()<2)
      return in.Rows()!=0;

   ulong rows_end = (in.Rows())>>1;

   for(ulong i = 0; i<rows_end; ++i)
      if(!in.SwapRows(i,rows-i))
         return false;

   return true;
  }
//+------------------------------------------------------------------+
//|  assign values to a portion of matrix selected within range      |
//|  of columns and                                                  |
//|  optional interval                                               |
//|  with support for negative indexing python style                 |
//| dest_start - column index of first column (inclusive)            |
//| dest_stop  - column index past last columun (exclusive)          |
//| dest_step  - column step size                                    |
//| indexes refer to positions in copyto                             |
//+------------------------------------------------------------------+
template <typename T>
bool matrixCopyCols(matrix<T> &copyto,matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
   return matrixCopy(copyto,copyfrom,BEGIN,END,STEP,dest_start,dest_stop,dest_step);
  }
//+------------------------------------------------------------------+
//|  assign values to a portion of matrix selected within range      |
//|  of rows and                                                     |
//|  optional interval                                               |
//|  with support for negative indexing python style                 |
//| dest_start - row index of first row (inclusive)                  |
//| dest_stop  - row index past last row (exclusive)                 |
//| dest_step  - row step size                                       |
//| indexes refer to positions in copyto                             |
//+------------------------------------------------------------------+
template <typename T>
bool matrixCopyRows(matrix<T> &copyto, matrix &copyfrom, const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
   return matrixCopy(copyto,copyfrom,dest_start,dest_stop,dest_step);
  }
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as  collection      |
//| of indices                                                       |
//+------------------------------------------------------------------+
template <typename T, typename P>
vector<T> select(vector<T> &in,const P &indices[])
  {
   P size = (P)MathMin(in.Size(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=P(in.Size()))
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return  vector<T>::Zeros(1);
     }

   vector<T> out(size);

   for(P i=0; i<size; ++i)
      out[i] = in[indices[i]];

   return out;
  }
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as  collection      |
//| of indices                                                       |
//+------------------------------------------------------------------+
template <typename T>
vector<T> select(vector<T> &in,vector<T> &indices)
  {
   ulong size = (ulong)MathMin(in.Size(),indices.Size());

   if(size<1 || indices.ArgMax()>=in.Size())
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return  vector<T>::Zeros(1);
     }

   vector<T> out(size);

   for(ulong i=0; i<size; ++i)
      out[i] = in[ulong(indices[i])];

   return out;
  }
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as  collection      |
//| of indices                                                       |
//+------------------------------------------------------------------+
template <typename T>
vector<T> select(vector<T> &in,const ulong &indices[])
  {
   ulong size = (ulong)MathMin(in.Size(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=ulong(in.Size()))
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return  vector<T>::Zeros(1);
     }

   vector<T> out(size);

   for(ulong i=0; i<size; ++i)
      out[i] = in[indices[i]];

   return out;
  }
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as  collection      |
//| of indices                                                       |
//+------------------------------------------------------------------+
template <typename T>
vector<T> select(vector<T> &in,const int &indices[])
  {
   int size = (int)MathMin(in.Size(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=int(in.Size()))
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return  vector<T>::Zeros(1);
     }

   vector<T> out(size);

   for(int i=0; i<size; ++i)
      out[i] = in[indices[i]];

   return out;
  }
//+------------------------------------------------------------------+
//|return vector of arbitrary elements specified as  collection      |
//| of indices                                                       |
//+------------------------------------------------------------------+
template <typename T>
vector<T> select(vector<T> &in,const uint &indices[])
  {
   uint size = (uint)MathMin(in.Size(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=uint(in.Size()))
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return  vector<T>::Zeros(1);
     }

   vector<T> out(size);

   for(uint i=0; i<size; ++i)
      out[i] = in[indices[i]];

   return out;
  }

//+------------------------------------------------------------------+
//|return  arbitrary columns of a matrix specified as                |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixCols(matrix<T> &in,const long &indices[])
  {
   ulong size = (ulong)MathMin(in.Cols(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=long(in.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(in.Rows(),size);

   for(ulong i=0; i<size; ++i)
      out.Col(in.Col(indices[i]),i);

   return out;
  }
//+------------------------------------------------------------------+
//|return  arbitrary rows of a matrix specified as                   |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixRows(matrix<T> &in,const long &indices[])
  {
   ulong size = (ulong)MathMin(in.Rows(),indices.Size());

   if(size<1  || indices[ArrayMaximum(indices,0,int(size))]>=long(in.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(size,in.Cols());

   for(ulong i=0; i<size; ++i)
      out.Row(in.Row(indices[i]),i);

   return out;
  }
//+------------------------------------------------------------------+
//|return  arbitrary columns of a matrix specified as                |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixCols(matrix<T> &in,const ulong &indices[])
  {
   ulong size = (ulong)MathMin(in.Cols(),indices.Size());

   if(size<1 || indices[ArrayMaximum(indices,0,int(size))]>=ulong(in.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(in.Rows(),size);

   for(ulong i=0; i<size; ++i)
      out.Col(in.Col(indices[i]),i);

   return out;
  }
//+------------------------------------------------------------------+
//|return  arbitrary rows of a matrix specified as                   |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixRows(matrix<T> &in,const ulong &indices[])
  {
   ulong size = (ulong)MathMin(in.Rows(),indices.Size());

   if(size<1  || indices[ArrayMaximum(indices,0,int(size))]>=ulong(in.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(size,in.Cols());

   for(ulong i=0; i<size; ++i)
      out.Row(in.Row(indices[i]),i);

   return out;
  }
//+------------------------------------------------------------------+
//|return  arbitrary columns of a matrix specified as                |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixCols(matrix<T> &in,vector& indices)
  {
   ulong size = (ulong)MathMin(in.Cols(),indices.Size());

   if(size<1 || indices.Max()>=double(in.Cols()) || indices.Min()<0.0)
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(in.Rows(),size);

   for(ulong i=0; i<size; ++i)
      out.Col(in.Col(ulong(indices[i])),i);

   return out;
  }
//+------------------------------------------------------------------+
//|return  arbitrary rows of a matrix specified as                   |
//|collection of indices                                             |
//+------------------------------------------------------------------+
template <typename T>
matrix<T> selectMatrixRows(matrix<T> &in,vector& indices)
  {
   ulong size = (ulong)MathMin(in.Rows(),indices.Size());

   if(size<1  || indices.Max()>=double(in.Rows()) || indices.Min()<0.0)
     {
      Print(__FUNCTION__," invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out(size,in.Cols());

   for(ulong i=0; i<size; ++i)
      out.Row(in.Row(ulong(indices[i])),i);

   return out;
  }
//+------------------------------------------------------------------+
//|copy vector copyfrom into specified positions of vector copyto    |
//|with support for negative indexing python style                   |
//+------------------------------------------------------------------+
template <typename T>
bool vectorCopy(vector<T> &copyto,vector<T> &copyfrom,const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
//---local variables
   long st,sp,stp,size;
//---check function parameters
   st = internal::normalizestart(dest_start,copyto.Size());
   sp = internal::normalizestop(dest_stop,copyto.Size());
   stp = dest_step;
//---error handling
   if(internal::invalidparams(st,sp,stp,copyto.Size()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return false;
     }
//---calculate size of new vector
   size = internal::calculatesize(st,sp,stp);
//---error handling
   if(copyfrom.Size()<ulong(size))
     {
      Print(__FUNCTION__," source vector is of insufficient size ");
      return false;
     }
//---assign elements
   for(long pos = st,index=0; index<size; ++index,pos+=stp)
      copyto[pos] = copyfrom[index];
//---done
   return true;
  }

//+------------------------------------------------------------------+
//| assign selected elements of vector the value /fill_value/        |
//+------------------------------------------------------------------+
template <typename T>
bool vectorFill(vector<T> &copyto, T fill_value,const long dest_start=BEGIN, const long dest_stop = END, const long dest_step = STEP)
  {
//---local variables
   long st,sp,stp,size;
//---check function parameters
   st = internal::normalizestart(dest_start,copyto.Size());
   sp = internal::normalizestop(dest_stop,copyto.Size());
   stp = dest_step;
//---error handling
   if(internal::invalidparams(st,sp,stp,copyto.Size()))
     {
      Print(__FUNCTION__," invalid function parameters ");
      return false;
     }
//---calculate size of new vector
   size = internal::calculatesize(st,sp,stp);
//---assign elements
   for(long pos = st,index=0; index<size; ++index,pos+=stp)
      copyto[pos] = fill_value;
//---done
   return true;
  }
//+------------------------------------------------------------------+
//|copy matrix copyfrom into specified positions of matrix copyto    |
//|with support for negative indexing python style                   |
//+------------------------------------------------------------------+
template<typename T>
bool matrixCopy(matrix<T> &copyto,matrix &copyfrom, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP,const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
  {
//---local variables
   long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
   rst =internal::normalizestart(row_dest_start,copyto.Rows());
   rsp =internal::normalizestop(row_dest_stop,copyto.Rows());
   rstp = row_dest_step;
//---error handling
   if(internal::invalidparams(rst,rsp,rstp,copyto.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters for row specification ");
      return false;
     }
//---calculate size of new matrix rows
   rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
   cst = internal::normalizestart(column_dest_start,copyto.Cols());
   csp = internal::normalizestop(column_dest_stop,copyto.Cols());
   cstp = column_dest_step;
//---error handling
   if(internal::invalidparams(cst,csp,cstp,copyto.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters for column specification ");
      return false;
     }
//---calculate size of new matrix columns
   csize = internal::calculatesize(cst,csp,cstp);
//---error handling
   if(copyfrom.Cols() < ulong(csize) || copyfrom.Rows() < ulong(rsize))
     {
      Print(__FUNCTION__, " source matrix has insufficient number of columns and/or rows ");
      return false;
     }
//---assign elements
   for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
      for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
         copyto[rpos][cpos]=copyfrom[row_index][column_index];
//---done
   return true;
  }
//+------------------------------------------------------------------+
//| assign selected elements of vector the value /fill_value/        |
//+------------------------------------------------------------------+
template<typename T>
bool matrixFill(matrix<T> &copyto,T fill_value, const long row_dest_start=BEGIN, const long row_dest_stop = END, const long row_dest_step = STEP,const long column_dest_start=BEGIN, const long column_dest_stop = END, const long column_dest_step = STEP)
  {
//---local variables
   long rst,rsp,rstp,rsize,cst,csp,cstp,csize;
//---check function parameters for row manipulation
   rst =internal::normalizestart(row_dest_start,copyto.Rows());
   rsp =internal::normalizestop(row_dest_stop,copyto.Rows());
   rstp = row_dest_step;
//---error handling
   if(internal::invalidparams(rst,rsp,rstp,copyto.Rows()))
     {
      Print(__FUNCTION__," invalid function parameters for row specification ");
      return false;
     }
//---calculate size of new matrix rows
   rsize = internal::calculatesize(rst,rsp,rstp);
//---check function parameters for row manipulation
   cst = internal::normalizestart(column_dest_start,copyto.Cols());
   csp = internal::normalizestop(column_dest_stop,copyto.Cols());
   cstp = column_dest_step;
//---error handling
   if(internal::invalidparams(cst,csp,cstp,copyto.Cols()))
     {
      Print(__FUNCTION__," invalid function parameters for column specification ");
      return false;
     }
//---calculate size of new matrix columns
   csize = internal::calculatesize(cst,csp,cstp);
//---assign elements
   for(long rpos = rst,row_index = 0; row_index<rsize; rpos+=rstp,++row_index)
      for(long cpos = cst,column_index = 0; column_index<csize; ++column_index,cpos+=cstp)
         copyto[rpos][cpos]=fill_value;
//---done
   return true;
  }
//+------------------------------------------------------------------+
//| returns vector                                                   |
//+------------------------------------------------------------------+
vector arange(const ulong size,  double value=0, double step_value=1)
  {
   vector out(size);

   for(ulong i = 0; i<size; ++i, value+=step_value)
      out[i] = value;

   return out;
  }
//+------------------------------------------------------------------+
//| returns vector                                                   |
//+------------------------------------------------------------------+
template<typename T>
vector<T> arange(T start_value, T stop_value, T step_value)
  {
   T numerator = (stop_value - start_value);
   ulong size = ulong((numerator)/step_value);
   if(MathMod(numerator,step_value)>T(0))
      size+=1;

   vector<T> out = vector<T>::Zeros(size);

   for(ulong i = 0; i<size; ++i, start_value+=step_value)
      out[i] = start_value;

   return out;
  }
//+------------------------------------------------------------------+
//| returns vector                                                   |
//+------------------------------------------------------------------+
vector range(int stop_value,  int start_value = 0, int step_value = 1)
  {
   double size = ceil(double(stop_value - start_value)/double(step_value));
   vector out(ulong(size));
   double value = double(start_value);
   double step = double(step_value);

   for(ulong i = 0; i<ulong(size); ++i, value+=step)
      out[i] = value;

   return out;
  }

//+------------------------------------------------------------------+
//| outputs array                                                    |
//+------------------------------------------------------------------+
template<typename T>
bool arange(T &in_out[],const int size,  T value = 0, T step_value = 1)
  {
   if(!size)
      return false;

   if(ArrayResize(in_out,fabs(size))!=fabs(size))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return false;
     }

   for(int i = 0; i<size; ++i, value+=step_value)
      in_out[i] = value;

   return true;
  }
//+------------------------------------------------------------------+
//| equivalent to numpy linspace(), outputs vector                   |
//+------------------------------------------------------------------+
template<typename T>
vector<T> linspace(T start, T stop, ulong num=50)
  {
   vector<T> out(num);

   if(!num)
      return vector<T>::Zeros(num);
   if(num==1)
     {
      out[0] = start;
      return out;
     }

   T chunk = (stop-start)/T(num-1);

   for(ulong i = 0; i<out.Size(); ++i)
      out[i] = start+(T(i)*chunk);

   return out;
  }
//+------------------------------------------------------------------+
//| fill diagonal                                                    |
//+------------------------------------------------------------------+
template<typename T>
bool fillDiagonal(matrix<T> &in_out, T fill_value)
  {
   if(in_out.Rows()!=in_out.Cols())
     {
      Print(__FUNCTION__ "Invalid Input. Matrix is not square ");
      return false;
     }

   for(ulong i = 0; i<in_out.Rows(); ++i)
      in_out[i][i] = fill_value;

   return true;
  }
//+------------------------------------------------------------------+
//| copy 3d container                                                |
//+------------------------------------------------------------------+
bool copy3D(matrix &copy_from[], matrix &copy_to[])
  {
   if(ArraySize(copy_to) != ArraySize(copy_from))
      ArrayResize(copy_to,(int)copy_from.Size());

   for(uint i = 0; i<copy_to.Size(); ++i)
      copy_to[i] = copy_from[i];

   return true;
  }
//+-----------------------------------------------------------------------+
//| Generate a Vandermonde matrix                                         |
//| in : vector                                                           |
//| num_columns : ulong,  default: size of in input vector.               |
//| ascending : bool, (default: false) Order of the powers of the columns.|
//| If true, the powers increase                                          |
//| from left to right, if false (the default) they are reversed          |
//+-----------------------------------------------------------------------+
template<typename T>
matrix<T> vander(vector<T> &in, ulong num_columns = 0,bool ascending = false)
  {
   ulong n = num_columns?num_columns:in.Size();

   matrix<T> out(in.Size(),n);

   for(ulong i = 0; i<n; ++i)
      out.Col(pow(in,ascending?i:n-1-i),i);

   return out;
  }

//+------------------------------------------------------------------+
//|generic vector to array                                           |
//+------------------------------------------------------------------+
template<typename T>
bool vecAsArray(const vector<T> &in, double &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i] = double(in[i]);
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic vector to array                                           |
//+------------------------------------------------------------------+
template<typename T>
bool vecAsArray(const vector<T> &in, int &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i] = int(in[i]);
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic vector to array                                           |
//+------------------------------------------------------------------+
template<typename T>
bool vecAsArray(const vector<T> &in, uint &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i] = uint(in[i]);
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic vector to array                                           |
//+------------------------------------------------------------------+
template<typename T>
bool vecAsArray(const vector<T> &in, long &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i] = long(in[i]);
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic vector to array                                           |
//+------------------------------------------------------------------+
template<typename T>
bool vecAsArray(const vector<T> &in, ulong &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i] = ulong(in[i]);
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic matrix to 1-D array                                       |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, double &out[])
  {
//---
   if(in.Cols()+in.Rows()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i]=double(in.Flat(i));
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic matrix to 1-D array                                       |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, int &out[])
  {
//---
   if(in.Cols()+in.Rows()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i]=int(in.Flat(i));
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic matrix to 1-D array                                       |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, uint &out[])
  {
//---
   if(in.Cols()+in.Rows()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i]=uint(in.Flat(i));
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic matrix to 1-D array                                       |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, long &out[])
  {
//---
   if(in.Cols()+in.Rows()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i]=long(in.Flat(i));
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//|generic matrix to 1-D array                                       |
//+------------------------------------------------------------------+
template<typename T, typename S>
bool matAsArray(const matrix<T> &in, ulong &out[])
  {
//---
   if(in.Cols()+in.Rows()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Cols()*in.Rows() && ArrayResize(out,int(in.Cols()*in.Rows()))!=int(in.Cols()*in.Rows()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); ++i)
      out[i]=ulong(in.Flat(i));
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//| difference vector according to parameter degreee                 |
//+------------------------------------------------------------------+
template<typename T>
vector<T> diff(vector<T> &in, ulong degree = 1)
  {
//---
   if(in.Size()<1 || degree<1 || in.Size()<degree+1)
     {
      Print(__FUNCTION__," Invalid parameter ");
      return vector<T>::Zeros(0);
     }
//---
   vector<T> yy,zz;
//---
   yy =  sliceVector(in,long(degree));
   zz =  sliceVector(in,BEGIN,long(in.Size()-degree));
//---
   return yy-zz;
  }
//+------------------------------------------------------------------+
//| difference matrix about the rows or columns                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> diff(matrix<T> &in, ulong degree = 1, bool by_row = true)
  {
//---
   if(in.Cols()+in.Rows()<1 || degree<1 || (by_row && in.Cols()<degree+1) || (!by_row && in.Rows()<degree+1))
     {
      Print(__FUNCTION__," Invalid function parameter ");
      return matrix::Zeros(0,0);
     }
//---
   matrix<T> yy,zz;
//---
   if(by_row)
     {
      yy =  sliceMatrix(in,BEGIN,END,STEP,long(degree));
      zz =  sliceMatrix(in,BEGIN,END,STEP,BEGIN,long(in.Cols()-degree));
     }
   else
     {
      yy =  sliceMatrix(in,long(degree));
      zz =  sliceMatrix(in,BEGIN,long(in.Rows()-degree));
     }
//---
   return yy-zz;
  }

//+------------------------------------------------------------------+
//| generate vector of num_repetitions of input vector elements      |
//+------------------------------------------------------------------+
template<typename T>
vector<T> repeat_vector_elements(vector<T> &vect,ulong num_repetitions)
  {
   vector<T> out(num_repetitions*vect.Size());

   for(ulong i = 0; i<out.Size(); i+=vect.Size())
      vectorCopy(out,vect,i,i+vect.Size());

   return out;
  }
//+------------------------------------------------------------------+
//| generate vector of num_repetitions of input matrix elements      |
//+------------------------------------------------------------------+
template<typename T>
vector<T> repeat_matrix_elements(matrix<T> &in, ulong num_repetitions)
  {
   vector<T> out(num_repetitions*in.Cols()*in.Rows());

   for(ulong i = 0; i<out.Size(); i+=(in.Cols()*num_repetitions))
     {
      vector copy = repeat_vector_elements(in.Row(i/(in.Cols()*num_repetitions)),num_repetitions);
      vectorCopy(out,copy,i,i+copy.Size());
     }

   return out;
  }
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as rows or columns           |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_vector_as_rows_cols(vector<T>&in, ulong num_repetitions, bool as_rows = true)
  {
   matrix<T> out(as_rows?num_repetitions:in.Size(),as_rows?in.Size():num_repetitions);

   for(ulong i = 0; i<num_repetitions; ++i)
      as_rows?out.Row(in,i):out.Col(in,i);

   return out;
  }
//+------------------------------------------------------------------+
//| generate matrix of num_repetitions of input columns or rows      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_matrix_rows_cols(matrix<T> &in, ulong num_repetitions, bool by_rows)
  {
   if(num_repetitions<1)
     {
      Print(__FUNCTION__, " invalid function parameter ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out((by_rows)?num_repetitions*in.Rows():in.Rows(),(!by_rows)?num_repetitions*in.Cols():in.Cols());

   if(by_rows)
      for(ulong i = 0; i<out.Rows(); ++i)
         out.Row(in.Row(i/num_repetitions),i);
   else
      for(ulong i = 0; i<out.Cols(); ++i)
         out.Col(in.Col(i/num_repetitions),i);

   return out;
  }

//+------------------------------------------------------------------+
//| generate matrix of an arbitrary number of repetitions of each row|
//| or column as specified by elements of input array                |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> repeat_matrix_rows_cols(matrix<T> &in, ulong &reps[],bool by_rows)
  {
   ulong max = 0 ;

   for(uint i = 0; i<reps.Size(); ++i)
      max+=reps[i];


   if((by_rows && in.Rows()!=ulong(reps.Size())) ||
      (!by_rows && in.Cols()!=ulong(reps.Size()))||
      reps[ArrayMinimum(reps)]<1)
     {
      Print(__FUNCTION__, " invalid function parameters ");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T> out((by_rows)?max:in.Rows(),(!by_rows)?max:in.Cols());

   if(by_rows)
      for(ulong i = 0, j = 0; i<out.Rows(); i+=reps[j],++j)
        {
         matrix sliced =sliceRows(in,j,j+1);
         matrix rpl = repeat_matrix_rows_cols(sliced,reps[j],true);
         copyRows(out,rpl,i,i+reps[j]);
        }
   else
      for(ulong i = 0, j = 0; i<out.Cols(); i+=reps[j],++j)
        {
         matrix sliced =sliceCols(in,j,j+1);
         matrix rpl = repeat_matrix_rows_cols(sliced,reps[j],false);
         copyCols(out,rpl,i,i+reps[j]);
        }

   return out;
  }
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as columns                   |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> vectorAsColumnMatrix(vector<T>&in, ulong num_cols)
  {
   return repeat_vector_as_rows_cols(in,num_cols,false);
  }
//+------------------------------------------------------------------+
//| generate matrix by repeating vector as rows                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> vectorAsRowMatrix(vector<T>&in, ulong num_rows)
  {
   return repeat_vector_as_rows_cols(in,num_rows,true);
  }
//+------------------------------------------------------------------+
//| The function extracts the unique values from the vector.         |
//|                                                                  |
//| Arguments:                                                       |
//| in     : vector with values                                      |
//|                                                                  |
//|                                                                  |
//| Return value: out vector of unique values  .                     |
//+------------------------------------------------------------------+
template<typename T>
vector<T> unique(vector<T> &in)
  {
//--- check array size
   ulong size=in.Size();
   if(size==0)
     {
      Print(__FUNCTION__," Invalid function parameter ");
      return vector<T>::Zeros(1);
     }

//--- prepare additional in
   bool checked[];
   if(ArrayResize(checked,int(size))!=int(size))
     {
      Print(__FUNCTION__," error ",GetLastError());
      return vector<T>::Zeros(0) ;
     }
   ArrayFill(checked,0,int(size),false);

//--- prepare out in
   vector<T>out(size);
//--- find unique elements
   ulong unique_count=0;
   T value=0;
   while(true)
     {
      bool flag=false;
      for(ulong i=unique_count; i<size; ++i)
        {
         if(!flag && !checked[i])
           {
            value=in[i];
            out[unique_count]=in[i];
            ++unique_count;
            checked[i]=true;
            flag=true;
           }
         else
            if(flag && value==in[i])
               checked[i]=true;
        }
      if(!flag)
         break;
     }
//--- resize target in
   out.Resize(unique_count);
//---
   return(out);
  }
//+------------------------------------------------------------------+
//| The function extracts the unique values from the vector.         |
//|                                                                  |
//| Arguments:                                                       |
//| in     : vector with values                                      |
//| out    : vector of unique values                                 |
//| indices: index positions in original vector                      |
//| counts : number of occurances of each unique value               |
//| Return bool: true on success and false on error                  |
//+------------------------------------------------------------------+
template<typename T>
bool  unique(vector<T> &in, vector<T> &out, long &indices[],long &counts[])
  {
//--- check vector size
   ulong size=in.Size();
   if(size==0)
     {
      Print(__FUNCTION__," Invalid function parameter ");
      return false;
     }

//--- prepare internal and output buffers
   bool checked[];
   if(ArrayResize(checked,int(size))!=int(size) ||
      ArrayResize(indices,int(size))!=int(size) ||
      ArrayResize(counts,int(size))!=int(size) ||
      !out.Resize(size))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return(false);
     }
   ArrayFill(checked,0,int(size),false);
   ArrayFill(indices,0,int(size),-1);
   ArrayFill(counts,0,int(size),1);

//--- prepare out vector
   out.Resize(size);
//--- find unique elements
   ulong unique_count=0;
   T value=0;
   while(true)
     {
      bool flag=false;
      for(ulong i=unique_count; i<size; ++i)
        {
         if(!flag && !checked[i])
           {
            value=in[i];
            out[unique_count]=in[i];
            indices[unique_count]=long(i);
            ++unique_count;
            checked[i]=true;
            flag=true;
           }
         else
            if(flag && value==in[i])
              {
               checked[i]=true;
               ++counts[unique_count-1];
              }
        }
      if(!flag)
         break;
     }
//---
   return (out.Resize(unique_count) && ArrayResize(indices,int(unique_count))==int(unique_count) && ArrayResize(counts,int(unique_count))==int(unique_count));
//---
  }
//+------------------------------------------------------------------+
//| The function extracts the unique columns or rows from the matrix |
//|                                                                  |
//| Arguments:                                                       |
//| in     : matrix with values                                      |
//| axis   : int by row = 1, column = 0                              |
//|                                                                  |
//| Return value: out matrix of unique rows or columns               |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> unique(matrix<T> &in,int axis = 0)
 {
   vector sum = in.Sum(axis);
   vector spec;
   long indices[],counts[];
   if(!unique(sum,spec,indices,counts))
    {
     Print(__FUNCTION__,": error ");
     return matrix<T>::Zeros(0,0);
    }
   if(axis)
    return selectMatrixRows(in,indices);
   else
    return selectMatrixCols(in,indices);
 }
//+------------------------------------------------------------------+
//|Integrate using the composite trapezoidal rule                    |
//+------------------------------------------------------------------+
template<typename T>
T trapezoid(vector<T>&y, vector<T>&x,T dx=1)
  {
   vector<T>d;

   if(x.Size()==0)
     {
      d = vector::Zeros(y.Size());
      d.Fill(dx);
     }
   else
      d = diff(x);

   vector<T> ret = (d*(sliceVector(y,1) + sliceVector(y,0,-1)))/2.0;

   return ret.Sum();
  }
//+---------------------------------------------------------------------------------+
//| One-dimensional linear interpolation for monotonically increasing sample points.|
//+---------------------------------------------------------------------------------+
template<typename T>
vector<T>interp(vector<T> &xv,vector<T>& xpoints,vector<T> &ypoints,double left=DBL_MIN,double right=DBL_MAX, double period=EMPTY_VALUE)
  {
   vector<T>out = vector<T>::Zeros(xv.Size());

   if(period != EMPTY_VALUE)
     {
      if(period == 0.0)
        {
         Print(__FUNCTION__," period should be non zero");
         return out;
        }

      vector<T> x = xv;
      vector<T> xp = xpoints;
      vector<T> yp = ypoints;

      vector vp(x.Size());
      vp.Fill(MathAbs(period));

      x = MathMod(xv,vp);
      xp = MathMod(xpoints,vp);

      long indices[];

      if(!arange(indices,int(xp.Size())) || !quickSortIndices(xp,true,indices,0,long(xp.Size()-1)))
        {
         Print(__FUNCTION__, " ", __LINE__);
         return out;
        }

      xp = select(xp,indices);
      yp = select(yp,indices);

      if(!xpoints.Resize((xp.Size()*2)+1) || !ypoints.Resize((yp.Size()*2)+1))
        {
         Print(__FUNCTION__, " ", __LINE__, " error ", GetLastError());
         return out;
        }

      if(!vectorCopy(xpoints,xp,xp.Size(),xp.Size()*2) || !vectorCopy(ypoints,yp,yp.Size(),yp.Size()*2))
        {
         Print(__FUNCTION__, " ", __LINE__);
         return out;
        }

      xpoints[xpoints.Size()-1] = xp[0]+MathAbs(period);
      ypoints[ypoints.Size()-1] = yp[0];

      if(!reverseVector(xp))
        {
         Print(__FUNCTION__, " ", __LINE__);
         return out;
        }

      vector temp = xp-vp;
      if(!vectorCopy(xpoints,temp,0,temp.Size()))
        {
         Print(__FUNCTION__, " ", __LINE__);
         return out;
        }

      if(!reverseVector(yp) && !vectorCopy(ypoints,yp,0,yp.Size()))
        {
         Print(__FUNCTION__, " ", __LINE__);
         return out;
        }

      xv = x;
     }

   double lval,rval;
   double xvals[];

   if(!vecAsArray(xpoints,xvals) || !ArraySort(xvals))
     {
      Print(__FUNCTION__, " ", __LINE__," error ", GetLastError());
      return out;
     }

   lval=(left==DBL_MIN)?ypoints[0]:left;
   rval=(right==DBL_MAX)?ypoints[ypoints.Size()-1]:right;

   if(xpoints.Size()==1)
     {
      for(ulong i =0; i<out.Size(); ++i)
         out[i] = xv[i]<xpoints[0]?lval:(xv[i]>xpoints[0])?rval:ypoints[0];
     }
   else
     {
      vector slopes=vector::Zeros(0);
      if(xpoints.Size()<=xv.Size())
         slopes = vector::Zeros(xpoints.Size()-1);

      if(slopes.Size())
        {
         vector ypd = sliceVector(ypoints,1) - sliceVector(ypoints,0,ypoints.Size()-1);
         vector xpd = sliceVector(xpoints,1) - sliceVector(xpoints,0,xpoints.Size()-1);
         slopes = ypd/xpd;
        }

      int j = 0;
      for(ulong i = 0; i<xv.Size(); ++i)
        {
         if(MathClassify(xv[i])==FP_NAN)
            out[i] = xv[i];
         j = internal::binary_search(xv[i],xvals);
         if(j == -1)
            out[i] = lval;
         else
            if(j == int(xpoints.Size()))
               out[i] = rval;
            else
               if(j == int(xpoints.Size()-1))
                  out[i] = ypoints[j];
               else
                  if(xpoints[j] == xv[i])
                     out[i] = ypoints[j];
                  else
                    {
                     double slope = slopes.Size()?slopes[j]:(ypoints[j+1] - ypoints[j]) / (xpoints[j+1] - xpoints[j]);
                     out[i] = slope*(xv[i] - xpoints[j])+ypoints[j];
                     if(MathClassify(out[i])==FP_NAN)
                       {
                        out[i] = slope*(xv[i] - xpoints[j+1]) + ypoints[j+1];
                        if(MathClassify(out[i])==FP_NAN && ypoints[j] == ypoints[j+1])
                           out[i] = ypoints[j];
                       }
                    }
        }
     }

   return out;

  }
//+------------------------------------------------------------------+
//| matrix multiplication                                            |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> matmul(const matrix<T>& matrix_a, const matrix<T>& matrix_b)
  {
   matrix<T> matrix_c = matrix<T>::Zeros(0,0);

   if(matrix_a.Cols()!=matrix_b.Rows())
      return(matrix_c);

   ulong M=matrix_a.Rows();
   ulong K=matrix_a.Cols();
   ulong N=matrix_b.Cols();
   matrix_c=matrix<T>::Zeros(M,N);

   for(ulong m=0; m<M; ++m)
      for(ulong k=0; k<K; ++k)
         for(ulong n=0; n<N; ++n)
            matrix_c[m][n]+=matrix_a[m][k]*matrix_b[k][n];

   return(matrix_c);
  }
//+------------------------------------------------------------------+
//| Matrix multiply with vector                                      |
//+------------------------------------------------------------------+
template<typename T>
vector<T> matmul(const matrix<T>& matrix_a, const vector<T>& vector_b)
  {
   matrix<T> matrix_b = matrix<T>::Zeros(vector_b.Size(),1);
   matrix_b.Col(vector_b,0);
   matrix<T>out = matmul(matrix_a,matrix_b);
   return out.Col(0);
  }
//+------------------------------------------------------------------+
//| Find indices where elements should be inserted to maintain order |
//+------------------------------------------------------------------+
template<typename T>
bool searchsorted(vector<T>&sorted,vector<T>&sample,bool right_side, vector &out)
  {
   out = vector::Zeros(sample.Size());

   vector<T> samplesorted = sample;

   if(!np::sort(samplesorted,sorted[0]<=sorted[sorted.Size()-1]))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return false;
     }


   for(ulong j = 0; j<sample.Size(); ++j)
     {
      for(ulong i = 1; i<sorted.Size(); ++i)
        {
         if((!right_side && sorted[i-1]<sample[j] && sample[j]<=sorted[i]) ||
            (right_side && sorted[i-1]<=sample[j] && sample[j]<sorted[i]))
           {
            out[j] = double(i);
            break;
           }
         else
           {
            if((!right_side && sample[j]<=sorted.Min()) || (right_side && sample[j]<sorted.Min()))
              {
               out[j] = 0;
               break;
              }
            else
              {
               if((!right_side && sample[j]>sorted.Max()) || (right_side && sample[j]>=sorted.Max()))
                 {
                  out[j] = double(sorted.Size());
                  break;
                 }
              }
           }
        }
     }

   return true;
  }
//+------------------------------------------------------------------+
//|matrix and vector addition with broadcasting                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> add(matrix<T>&mat, vector<T>&vec)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size() == mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(mat.Row(i)+vec,i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(mat.Col(i)+vec,i);

   return out;
  }
//+------------------------------------------------------------------+
//|matrix and vector subtraction with broadcasting                   |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> minus(matrix<T>&mat, vector<T>&vec)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size() == mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(mat.Row(i)-vec,i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(mat.Col(i)-vec,i);
   return out;
  }
//+------------------------------------------------------------------+
//|matrix and vector subtraction with broadcasting                   |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> minus(vector<T>&vec, matrix<T>&mat)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size()==mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(vec-mat.Row(i),i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(vec-mat.Col(i),i);
   return out;
  }
//+------------------------------------------------------------------+
//| vector subtraction with broadcasting                             |
//+------------------------------------------------------------------+
template<typename T>
vector<T> minus(vector<T>& va, vector<T>& vb)
  {
   if(va.Size()==vb.Size())
      return va-vb;
   else
     {
      if(va.Size()!=1 && vb.Size()!=1)
        {
         Print(__FUNCTION__, " operands could not be broadcast together ");
         return vector<T>::Zeros(0);
        }
      else
        {
         if(vb.Size()>va.Size())
            return va[0] - vb;
         else
            return va - vb[0];
        }
     }
  }
//+------------------------------------------------------------------+
//| vector addition with broadcasting                                |
//+------------------------------------------------------------------+
template<typename T>
vector<T> add(vector<T>& va, vector<T>& vb)
  {
   if(va.Size()==vb.Size())
      return va+vb;
   else
     {
      if(va.Size()!=1 && vb.Size()!=1)
        {
         Print(__FUNCTION__, " operands could not be broadcast together ");
         return vector<T>::Zeros(0);
        }
      else
        {
         if(vb.Size()>va.Size())
            return va[0] + vb;
         else
            return va + vb[0];
        }
     }
  }
//+------------------------------------------------------------------+
//| vector multiplication with broadcasting                          |
//+------------------------------------------------------------------+
template<typename T>
vector<T> multiply(vector<T>& va, vector<T>& vb)
  {
   if(va.Size()==vb.Size())
      return va*vb;
   else
     {
      if(va.Size()!=1 && vb.Size()!=1)
        {
         Print(__FUNCTION__, " operands could not be broadcast together ");
         return vector<T>::Zeros(0);
        }
      else
        {
         if(vb.Size()>va.Size())
            return va[0]*vb;
         else
            return va*vb[0];
        }
     }
  }
//+------------------------------------------------------------------+
//| vector division with broadcasting                                 |
//+------------------------------------------------------------------+
template<typename T>
vector<T> divide(vector<T>& va, vector<T>& vb)
  {
   if(va.Size()==vb.Size())
      return va/vb;
   else
     {
      if(va.Size()!=1 && vb.Size()!=1)
        {
         Print(__FUNCTION__, " operands could not be broadcast together ");
         return vector<T>::Zeros(0);
        }
      else
        {
         if(vb.Size()>va.Size())
            return va[0]/vb;
         else
            return va/vb[0];
        }
     }
  }
//+------------------------------------------------------------------+
//|matrix and vector mulitiplication with broadcasting               |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> multiply(matrix<T>&mat, vector<T>&vec)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size() == mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(mat.Row(i)*vec,i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(mat.Col(i)*vec,i);
   return out;
  }
//+------------------------------------------------------------------+
//|matrix and vector division with broadcasting                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> divide(matrix<T>&mat, vector<T>&vec)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size() == mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(mat.Row(i)/vec,i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(mat.Col(i)/vec,i);
   return out;
  }
//+------------------------------------------------------------------+
//|matrix and vector division with broadcasting                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> divide(vector<T>&vec, matrix<T>&mat)
  {
   if(vec.Size()!=mat.Cols() && vec.Size()!=mat.Rows())
     {
      Print(__FUNCTION__, " operands could not be broadcast together ", " (",mat.Rows(),":",mat.Cols(), ") and (", vec.Size(),")");
      return matrix<T>::Zeros(0,0);
     }

   matrix<T>out=matrix::Zeros(mat.Rows(),mat.Cols());
   if(vec.Size() == mat.Cols())
      for(ulong i = 0; i<mat.Rows(); ++i)
         out.Row(vec/mat.Row(i),i);
   else
      for(ulong i = 0; i<mat.Cols(); ++i)
         out.Col(vec/mat.Col(i),i);
   return out;
  }

//+------------------------------------------------------------------+
//|  Flattens the array                                              |
//+------------------------------------------------------------------+
template<typename T>
vector<T>ravelMultiIndex(vector<T> &in[], ulong&dims[])
  {
   vector<T> vdims;

   if(!vdims.Assign(dims))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(1);
     }

   vector<T> vec = np::sliceVector(vdims,1);

   if(!np::reverseVector(vec))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(1);
     }

   vec = vec.CumProd();

   if(!np::reverseVector(vec))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(1);
     }

   vector<T> nvec = vector::Ones(vec.Size()+1);

   if(!vectorCopy(nvec,vec,1))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return vector::Zeros(1);
     }

   matrix<T>out(in.Size(),in[0].Size());

   for(ulong i = 0; i<out.Rows(); ++i)
      if(!out.Row(in[i],i))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return vector::Zeros(1);
        }

   return nvec.MatMul(out);

  }

//+------------------------------------------------------------------+
//| get the bin edges                                                |
//+------------------------------------------------------------------+
template<typename T>
vector<T> binEdges(vector<T>&in, ulong numbins=3)
  {
   T first_edge = in.Min();
   T last_edge = in.Max();

   return np::linspace(first_edge,last_edge,numbins);
  }
/*+------------------------------------------------------------------+
//| numpy digitize for scalar values                                 |
//+------------------------------------------------------------------+
template<typename T>
long digitize(T sample, vector<T>&edges, bool right_side=false)
  {
   for(long i = 1; i<long(edges.Size()); ++i)
     {
      if((right_side && edges[i-1]<sample && sample<=edges[i]) ||
         (!right_side && edges[i-1]<=sample && sample<edges[i]))
         return i;
      else
        {
         if((right_side && sample<=edges.Min()) || (!right_side && sample<edges.Min()))
            return i;
         else
           {
            if((right_side && sample>edges.Max()) || (!right_side && sample>=edges.Max()))
               return long(edges.Size());
           }
        }
     }
   return -1;
  }
//+------------------------------------------------------------------+
//| Compute the multidimensional histogram of some data.             |
//+------------------------------------------------------------------+
template<typename T>
bool histogramdd(matrix<T>&in, ulong bins, vector &hist, vector &edges[])
  {
   hist = vector::Zeros(ulong(pow(bins,in.Cols())));

   if(edges.Size())
      ArrayFree(edges);

   if(!bins)
      bins=3;

   if(ArrayResize(edges,int(in.Cols()))<0)
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return false;
     }

   for(ulong i = 0; i<in.Cols(); ++i)
      edges[i] = np::linspace(in.Col(i).Min(), in.Col(i).Max(),bins+1);

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      vector<T>row = in.Row(i);
      long final_index=0;
      for(ulong j =0; j<row.Size(); ++j)
        {
         long index = np::digitize(row[j],edges[j]) - 1;
         if(ulong(index) == (edges[j].Size()-1))
            index -=1;
         final_index+=(index*long(pow(bins,row.Size()-(j+1))));
        }
      if(final_index>=long(hist.Size()) || final_index<0)
        {
         Print(__FUNCTION__, " hist index out of bounds. Index ", final_index, " hist size ", hist.Size());
         return false;
        }
      hist[final_index] +=1.0;
     }

   return true;
  }
//+------------------------------------------------------------------+
//| Compute the multidimensional histogram of some data.             |
//+------------------------------------------------------------------+
template<typename T>
bool histogramdd(matrix<T>&in, ulong &bins[], vector &hist, vector &edges[])
  {
   hist = vector::Zeros(np::internal::prodArray(bins));

   if(edges.Size())
      ArrayFree(edges);

   if(ulong(bins.Size())<in.Cols())
     {
      Print(__FUNCTION__, " invalid parameter, bins array is of insufficient size ");
      return false;
     }

   if(bins[ArrayMinimum(bins)]==0)
     {
      Print(__FUNCTION__, " invalid parameter, bins array contains zero value ");
      return false;
     }

   if(ArrayResize(edges,int(in.Cols()))<0)
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return false;
     }

   for(ulong i = 0; i<in.Cols(); ++i)
      edges[i] = np::linspace(in.Col(i).Min(), in.Col(i).Max(),bins[i]+1);


   for(ulong i = 0; i<in.Rows(); ++i)
     {
      vector<T>row = in.Row(i);
      long final_index=0;
      for(ulong j =0; j<row.Size(); ++j)
        {
         long index = np::digitize(row[j],edges[j]) - 1;
         if(ulong(index) == (edges[j].Size()-1))
            index -=1;
         final_index+=(index*long(np::internal::prodArray(bins,uint(j+1))));
        }

      if(final_index>=long(hist.Size()) || final_index<0)
        {
         Print(__FUNCTION__, " hist index out of bounds. Index ", final_index, " hist size ", hist.Size());
         return false;
        }

      hist[final_index] +=1.0;
     }

   return true;
  }

//+------------------------------------------------------------------+
//| Compute the one dimensional histogram of some data.             |
//+------------------------------------------------------------------+
template<typename T>
bool histogram(vector<T>&in, ulong bins, vector &hist, vector &edges)
  {
   hist = vector::Zeros(in.Size());

   if(!bins)
      bins=3;

   if(in.Size()/2 < bins)
     {
      Print(__FUNCTION__, " invalid parameters, too few observations relative to the number of bins ");
      return false;
     }

   edges = np::linspace(in.Min(), in.Max(),bins+1);

   for(ulong j =0; j<in.Size(); ++j)
     {
      long index = np::digitize(in[j],edges) - 1;
      if(ulong(index) == (edges.Size()-1))
         index -=1;
      hist[index] +=1.0;
     }

   return true;
  }
//+------------------------------------------------------------------+
//|Computes empirical quantiles for a matrix.                        |
//+------------------------------------------------------------------+
matrix quantiles(matrix &in,vector &probs, double alphap=0.4,double betap=0.4, ulong axis = 0, double lower=-DBL_MIN, double upper=DBL_MAX)
  {
   matrix d = in;

   if(!d.Clip(lower,upper))
     {
      Print(__FUNCTION__, " Clip() failure ", GetLastError());
      return matrix::Zeros(probs.Size(),in.Cols());
     }

   vector m = alphap + probs*(1.0-alphap-betap);

   matrix out(axis?in.Rows():probs.Size(),axis?probs.Size():in.Cols());

   ulong step = (axis)?in.Rows():in.Cols();

   for(ulong i = 0; i<step; ++i)
     {
      vector vec = (axis)?d.Row(i):d.Col(i);
      vector outvec = np::internal::quantile(vec,m,probs);
      if((axis && !out.Row(outvec,i)) || (!axis && !out.Col(outvec,i)))
        {
         Print(__FUNCTION__, " vector insertion failure ", GetLastError());
         return matrix::Zeros(out.Rows(),out.Cols());
        }
     }

   return out;
  }*/
//+------------------------------------------------------------------+
//|BitwiseNot                                                        |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseNot(vector<T>&left)
  {
   int aleft[];
   int res[];

   if(left.Size()<1)
     {
      Print(__FUNCTION__, " invalid function input ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size());
     }

   if(!MathBitwiseNot(aleft,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseAnd                                                        |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseAnd(vector<T>&left,vector<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(left.Size()!=right.Size() || right.Size()<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size());
     }

   if(!MathBitwiseAnd(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseOr                                                         |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseOr(vector<T>&left,vector<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(left.Size()!=right.Size() || right.Size()<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size())
     }

   if(!MathBitwiseOr(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseXor                                                        |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseXor(vector<T>&left,vector<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(left.Size()!=right.Size() || right.Size()<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft) || !vecAsArray(right,aright))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size())
     }

   if(!MathBitwiseXor(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseShiftL                                                     |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseShiftL(vector<T>&left,int shift)
  {
   int aleft[];
   int res[];

   if(!left.Size() || shift<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size())
     }

   if(!MathBitwiseShiftL(aleft,shift,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseShiftR                                                     |
//+------------------------------------------------------------------+
template<typename T>
vector<T>bitwiseShiftR(vector<T>&left,int shift)
  {
   int aleft[];
   int res[];

   if(!left.Size() || shift<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return vector::Zeros(left.Size());
     }

   if(!vecAsArray(left,aleft))
     {
      Print(__FUNCTION__, " vector conversion failed ");
      return vector::Zeros(left.Size())
     }

   if(!MathBitwiseShiftR(aleft,shift,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   vector<T> vec;

   if(!vec.Assign(res))
     {
      Print(__FUNCTION__ " Array assignment to vector failed ", GetLastError());
      return vector::Zeros(left.Size());
     }

   return vec;
  }
//+------------------------------------------------------------------+
//|BitwiseNot                                                        |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseNot(matrix<T>&left)
  {
   int aleft[];
   int res[];

   if(!left.Rows() || !left.Cols())
     {
      Print(__FUNCTION__, " invalid function input ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseNot(aleft,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }
//+------------------------------------------------------------------+
//|BitwiseAnd                                                        |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseAnd(matrix<T>&left,matrix<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft) || !matAsArray(right,aright))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseAnd(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }
//+------------------------------------------------------------------+
//|BitwiseOr                                                         |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseOr(matrix<T>&left,matrix<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft) || !matAsArray(right,aright))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseOr(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }
//+------------------------------------------------------------------+
//|BitwiseXor                                                        |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseXor(matrix<T>&left,matrix<T>&right)
  {
   int aleft[];
   int aright[];
   int res[];

   if(!left.Rows() || !left.Cols() || !right.Cols() || !right.Rows())
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft) || !matAsArray(right,aright))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseXor(aleft,aright,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }
//+------------------------------------------------------------------+
//|BitwiseShiftL                                                     |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseShiftL(matrix<T>&left,int shift)
  {
   int aleft[];
   int res[];

   if(!left.Rows() || !left.Cols() || shift<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseShiftL(aleft,shift,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }
//+------------------------------------------------------------------+
//|BitwiseShiftR                                                     |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>bitwiseShiftR(matrix<T>&left,int shift)
  {
   int aleft[];
   int res[];

   if(!left.Rows() || !left.Cols() || shift<1)
     {
      Print(__FUNCTION__, " invalid function inputs ");
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   if(!matAsArray(left,aleft))
     {
      Print(__FUNCTION__, " matrix conversion failed ");
      return matrix::Zeros(left.Rows(),left.Cols())
     }

   if(!MathBitwiseShiftR(aleft,shift,res))
     {
      Print(__FUNCTION__ " MathBitwiseNot failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   matrix<T> mat;

   if(!mat.Assign(res) || !mat.Reshape(left.Rows(),left.Cols()))
     {
      Print(__FUNCTION__ " Array assignment to matrix failed ", GetLastError());
      return matrix::Zeros(left.Rows(),left.Cols());
     }

   return mat;
  }

//+------------------------------------------------------------------+
//| write matrix to csv file                                         |
//+------------------------------------------------------------------+
template <typename T>
bool matrix2csv(string csv_name, matrix<T> &matrix_, string header_string="", bool common=false, int digits=8)
  {
   FileDelete(csv_name);
   int handle = FileOpen(csv_name,FILE_WRITE|FILE_CSV|FILE_ANSI|(common?FILE_COMMON:FILE_ANSI),",",CP_UTF8);

   if(header_string == "")
      for(ulong i=0; i<matrix_.Cols(); ++i)
         header_string += IntegerToString(long(i)) + (i==matrix_.Cols()-1?"":",");

   if(handle == INVALID_HANDLE)
     {
      printf("Invalid %s handle Error %d ",csv_name,GetLastError());
      return (false);
     }

   string concstring;
   vector<T> row = {};

   datetime time_start = GetTickCount(), current_time;

   string header[];

   ushort u_sep;
   u_sep = StringGetCharacter(",",0);
   StringSplit(header_string,u_sep, header);

   vector<T> colsinrows = matrix_.Row(0);

   if(ArraySize(header) != (int)colsinrows.Size())
     {
      printf("headers=%d and columns=%d from the matrix vary is size ",ArraySize(header),colsinrows.Size());
      return false;
     }

//---

   string header_str = "";
   for(int i=0; i<ArraySize(header); ++i)
      header_str += header[i] + (i+1 == colsinrows.Size() ? "" : ",");

   FileWrite(handle,header_str);

   FileSeek(handle,0,SEEK_SET);

   for(ulong i=0; i<matrix_.Rows() && !IsStopped(); ++i)
     {
      ZeroMemory(concstring);

      row = matrix_.Row(i);
      for(ulong j=0, cols =1; j<row.Size() && !IsStopped(); ++j, ++cols)
        {
         current_time = GetTickCount();

         concstring += (string)NormalizeDouble(row[j],digits) + (cols == matrix_.Cols() ? "" : ",");
        }

      FileSeek(handle,0,SEEK_END);
      FileWrite(handle,concstring);
     }

   FileClose(handle);

   return (true);
  }
//+------------------------------------------------------------------+
//| read csv file into matrix                                        |
//+------------------------------------------------------------------+
matrix readcsv(string file_name,bool common=false,string delimiter=",",bool handle_date_column=false,int date_column_index=-1,bool skip_header=true)
  {
   string Arr[],csv_header[];
   int all_size = 0;

   int cols_total=0;

   int handle = FileOpen(file_name,FILE_SHARE_READ|FILE_CSV|FILE_ANSI|(common?FILE_COMMON:FILE_ANSI),delimiter);

   if(handle == INVALID_HANDLE)
     {
      printf("Invalid %s handle Error %d ",file_name,GetLastError());
      Print(GetLastError()==0?" TIP | File Might be in use Somewhere else or in another Directory":"");
     }
   else
     {
      int column = 0;
      int rows = skip_header?0:1;

      while(!FileIsEnding(handle) && !IsStopped())
        {
         string data = FileReadString(handle);

         //---

         if((rows == 0 && skip_header) || (rows == 1 && !skip_header))
           {
            ArrayResize(csv_header,column+1);
            csv_header[column] = skip_header?data:string(column);
           }

         ++column;

         if(rows>0)  //Avoid the first column which contains the column's header
           {
            ++all_size;
            ArrayResize(Arr,all_size);

            Arr[all_size-1] = data;

           }
         //---

         if(FileIsLineEnding(handle))
           {
            cols_total=column;
            ++rows;
            column = 0;
           }
        }

      FileClose(handle);
     }

   int rows =all_size/cols_total;

   Comment("");

   matrix mat(rows, cols_total);
   string Col[];
   vector col_vector;

   for(int i=0; i<cols_total; ++i)
     {
      internal::col(Arr, Col, i+1, cols_total);
      mat.Col(internal::formatCol(Col,handle_date_column,date_column_index,i), i);
     }

   return(mat);
  }
//+------------------------------------------------------------------+
//| read csv file into matrix                                        |
//+------------------------------------------------------------------+
matrix readcsv_from_string(string csv_string_data,bool common=false,string delimiter=",",bool handle_date_column=false,int date_column_index=-1,bool skip_header=true)
  {
   string Arr[],csv_header[],lines[];
   int all_size = 0;
   int rows = 0;
   int cols_total=0;
   
   int handle = StringSplit(csv_string_data,StringGetCharacter("\n",0),lines);

   if(handle <= 0)
     {
      Print(__FUNCTION__, " no lines found in the string data ");
      return matrix::Zeros(0,0);
     }
   else
     {
      int column = 0;
      int shift  = 0;
      cols_total = StringSplit(lines[shift],StringGetCharacter(delimiter,0),csv_header);
      if(StringLen(csv_header[cols_total-1]) == 0)
       --cols_total;
      all_size = handle - 1;
      if(StringLen(lines[handle-1]) == 0)
        --all_size;
      rows = all_size;
      all_size*=cols_total;
      ArrayResize(Arr,all_size);
      for(int i = 0; i<(all_size); i+=cols_total)
       {
        StringSplit(lines[skip_header?((i+cols_total)/cols_total):(i/cols_total)],StringGetCharacter(delimiter,0),csv_header);
        ArrayCopy(Arr,csv_header,i,0,cols_total);
       }
     }
     
   matrix mat(rows, cols_total);
   string Col[];
   vector col_vector;

   for(int i=0; i<cols_total; ++i)
     {
      internal::col(Arr, Col, i+1, cols_total);
      mat.Col(internal::formatCol(Col,handle_date_column,date_column_index,i), i);
     }

   return(mat);
  }
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the       |
//| elements of the vector.                                          |
//|                                                                  |
//| Arguments:                                                       |
//| in     : vector                                                  |
//| count  : Number of elements to choose                            |
//| result : Output vector                                           |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool sampleVector(vector<T> &in,const ulong count,vector<T> &result)
  {
   if(!in.Size())
      return false;
//--- prepare target array and calculate values
   if(result.Size()<count)
      if(!result.Resize(count))
         return(false);
   for(ulong i=0; i<count; ++i)
     {
      int ind=(int)(in.Size()*MathRand()/32768);
      result[i]=in[ind];
     }
//---
   return true;
  }
//+------------------------------------------------------------------+
//| Generates a random sample from a given 1-D array a               |
//+------------------------------------------------------------------+
template<typename T>
vector<T> choice(vector<T>& a, vector& p, ulong size = 0, bool replace = true)
  {
   ulong pop_size = a.Size();
   vector idx;
   if(!pop_size && size)
     {
      Print(__FUNCTION__, " a cannot be empty unless no samples are taken ");
      return vector<T>::Zeros(0);
     }
   if(p.Size())
     {
      ulong d = p.Size();
      double atol = DBL_EPSILON;
      if(d!=pop_size)
        {
         Print(__FUNCTION__, " a and p must have same size ");
         return vector<T>::Zeros(0);
        }
      double psum = p.Sum();
      if(p.Min()<0.0)
        {
         Print(__FUNCTION__, " probabilities cannot be negative ");
         return vector<T>::Zeros(0);
        }
      if(psum-1.>atol)
        {
         Print(__FUNCTION__, " probabilities must sum to 1 ");
         return vector<T>::Zeros(0);
        }
     }
   CHighQualityRandStateShell rstate;
   CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
   bool is_scalar = (size==0);
   ulong shape = 0;
   vector cdf;
   if(!is_scalar)
      shape = size;
   else
      size = 1;
   if(replace)
     {
      if(p.Size())
        {
         cdf = p.CumSum();
         cdf/=cdf[cdf.Size()-1];
         vector uniform_samples = vector::Zeros(shape);
         for(ulong i = 0; i<shape; uniform_samples[i] = CHighQualityRand::HQRndUniformR(rstate.GetInnerObj()), ++i);
         if(!searchsorted(cdf,uniform_samples,true,idx))
           {
            Print(__FUNCTION__, " searchsorted failed ", __LINE__);
            return vector<T>::Zeros(0);
           }
        }
      else
        {
         idx = vector::Zeros(size);
         for(ulong i = 0; i<shape; idx[i] = (double)CHighQualityRand::HQRndUniformI(rstate.GetInnerObj(),int(pop_size)), ++i);
        }
     }
   else
     {
      if(size>pop_size)
        {
         Print(__FUNCTION__, " cannot take a larger sample thatn population when replace = false  ");
         return vector<T>::Zeros(0);
        }
      if(p.Size())
        {
         ulong n_uniq = 0;
         vector pp = p;
         vector found = vector::Zeros(shape);
         vector cdf,_new,temp;
         long uindices[],counted[];
         while(n_uniq<size)
           {
            vector x = vector::Zeros(size - n_uniq);
            for(ulong i = 0; i<x.Size(); x[i] = CHighQualityRand::HQRndUniformR(rstate.GetInnerObj()), ++i);
            if(n_uniq>0)
              {
               for(ulong i = 0; i<n_uniq; ++i)
                  pp[ulong(found[i])] = 0.0;
              }
            cdf = pp.CumSum();
            cdf /= cdf[cdf.Size()-1];
            if(!searchsorted(cdf,x,true,_new))
              {
               Print(__FUNCTION__, " searchsorted failed ", __LINE__);
               return vector::Zeros(0);
              }
            if(!unique(_new,temp,uindices,counted) || !ArraySort(uindices))
              {
               Print(__FUNCTION__, " unique failed ", __LINE__);
               return vector::Zeros(0);
              }
            _new = select(_new,uindices);
            for(ulong i = n_uniq, j = 0; i<(n_uniq+_new.Size()); ++i)
               found[i] = _new[j++];
            n_uniq+=_new.Size();
           }
         idx = found;
        }
      else
        {
         idx = permutation(pop_size);
         idx = sliceVector(idx,0,long(size));
         idx.Resize(shape);
        }
     }
   return select(a,idx);
  }
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the       |
//| elements of the vector with replacement.                        |
//|                                                                  |
//| Arguments:                                                       |
//| in     : vector with values                           |
//| count       : Number of elements to choose                       |
//| replace     : If true, Sampling with Replacement will be used    |
//| result    : Output vector                                       |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool sampleVector(vector<T> &in,const ulong count,const bool replace,vector<T> &result)
  {
//--- Sampling with Replacement
   if(replace)
      return sampleVector(in,count,result);

   if(!in.Size())
      return(false);
//--- prepare target array
   if(result.Size()<count)
      if(!result.Resize(count))
         return(false);
//--- unique values needed, prepare indices
   ulong indices[];

   if(!arange(indices,int(in.Size())))
      return(false);

   for(ulong i=0; i<count; ++i)
     {
      //--- select random index and swap items [j] and [i]
      ulong j=i+ulong((ulong)(in.Size()-i)*MathRand()/32768);

      if(j!=i)
        {
         ulong t=indices[i];
         indices[i]=indices[j];
         indices[j]=t;
        }
     }
//--- select data according to indices
   for(ulong i=0; i<count; ++i)
      result[i]=in[indices[i]];
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the       |
//| elements of the vector.                                          |
//|                                                                  |
//| Arguments:                                                       |
//| in     : vector                                                  |
//| count  : Number of elements to choose                            |
//| result : Output vector                                           |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool sampleArray(T &array[],const ulong count,T &result[])
  {
   if(!array.Size())
      return false;
//--- prepare target array and calculate values
   if(result.Size()<uint(count))
      if(ArrayResize(result,int(count))!=int(count))
         return(false);
   for(ulong i=0; i<count; ++i)
     {
      int ind=(int)(array.Size()*MathRand()/32768);
      result[i]=array[ind];
     }
//---
   return true;
  }
//+------------------------------------------------------------------+
//| The function takes a sample of the specified size from the       |
//| elements of the vector with replacement.                         |
//|                                                                  |
//| Arguments:                                                       |
//| array     : array with values                                    |
//| count       : Number of elements to choose                       |
//| replace     : If true, Sampling with Replacement will be used    |
//| result    : Output vector                                        |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool sampleArray(T &array[],const ulong count,const bool replace,T &result[])
  {
//--- Sampling with Replacement
   if(replace)
      return sampleArray(array,count,result);

   if(!array.Size())
      return(false);
//--- prepare target array
   if(result.Size()<uint(count))
      if(ArrayResize(result,int(count))!=int(count))
         return(false);
//--- unique values needed, prepare indices
   ulong indices[];

   if(!arange(indices,int(array.Size())))
      return(false);

   for(ulong i=0; i<count; ++i)
     {
      //--- select random index and swap items [j] and [i]
      ulong j=i+ulong((ulong)(array.Size()-i)*MathRand()/32768);

      if(j!=i)
        {
         ulong t=indices[i];
         indices[i]=indices[j];
         indices[j]=t;
        }
     }
//--- select data according to indices
   for(ulong i=0; i<count; ++i)
      result[i]=array[indices[i]];
//---
   return(true);
  }

//+------------------------------------------------------------------+
//| EmpiricalProbabilityDensityEmpirical                             |
//+------------------------------------------------------------------+
template<typename T>
bool epdf(const vector<T> &vec,const ulong count,vector<T> &x,vector<T> &pdf)
  {
   if(count<=1)
      return(false);
   ulong size=vec.Size();
   if(size==0)
      return(false);
//--- check NaN values
   for(ulong i=0; i<size; ++i)
     {
      if(!MathIsValidNumber(vec[i]))
         return(false);
     }
//--- prepare output arrays
   if(x.Size()<count)
      if(!x.Resize(count))
         return(false);
   if(pdf.Size()<count)
      if(pdf.Resize(count))
         return(false);
//--- search for min,max and range
   double minv=vec.Min();
   double maxv=vec.Max();
   double range=maxv-minv;
   if(range==0)
      return(false);
//--- calculate probability density of the empirical distribution
   for(ulong i=0; i<count; ++i)
     {
      x[i]=minv+i*range/(count-1);
      pdf[i]=0;
     }
   for(ulong i=0; i<size; ++i)
     {
      double v=(vec[i]-minv)/range;
      ulong ind=ulong((v*(count-1)));
      ++pdf[ind];
     }
//--- normalize values
   double dx=range/count;
   double sum=0;
   for(ulong i=0; i<count; ++i)
      sum+=pdf[i]*dx;
   if(sum==0)
      return(false);
   double coef=1.0/sum;
   for(ulong i=0; i<count; ++i)
      pdf[i]*=coef;
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| CumulativeDistributionEmpirical                                  |
//+------------------------------------------------------------------+
template<typename T>
bool ecdf(const vector <T> &vec,const ulong count,vector<T> &x,vector<T> &cdf)
  {
   if(count<1)
      return(false);
   ulong size=vec.Size();
   if(size==0)
      return(false);
//--- check NaN values
   if(vec.HasNan())
     {
      Print(__FUNCTION__," input data contains invalid number(s)");
      return false;
     }
//--- prepare output arrays
   if(x.Size()<count)
      if(!x.Resize(count))
         return(false);
   if(cdf.Size()<count)
      if(!cdf.Resize(count))
         return(false);
//--- search for min,max and range
   double minv=vec.Min();
   double maxv=vec.Max();
   double range=maxv-minv;
   if(range==0)
      return(false);
//--- calculate probability density of the empirical distribution
   vector<T> pdf;
   if(!pdf.Resize(count))
      return(false);
   for(ulong i=0; i<count; ++i)
     {
      x[i]=minv+i*range/(count-1);
      pdf[i]=0;
     }
   for(ulong i=0; i<size; ++i)
     {
      double v=(vec[i]-minv)/range;
      ulong ind=ulong((v*(count-1)));
      ++pdf[ind];
     }
//--- normalize values
   double dx=range/count;
   double sum=0;
   for(ulong i=0; i<count; ++i)
      sum+=pdf[i]*dx;
   if(sum==0)
      return(false);
   double coef=1.0/sum;
   for(ulong i=0; i<count; ++i)
      pdf[i]*=coef;
//--- calculate cumulative distribution function
   sum=0.0;
   for(ulong i=0; i<count; ++i)
     {
      sum+=pdf[i]*dx;
      cdf[i]=sum;
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| MathFactorial                                                    |
//+------------------------------------------------------------------+
double factorial(const uint n)
  {
   double ft[21]=
     {
      1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600,
      6227020800,87178291200,1307674368000,20922789888000,355687428096000,
      6402373705728000,121645100408832000,2432902008176640000
     };
   if(n<=20)
      //--- use values from the factorials table
      return(ft[n]);
   else
     {
      //--- calculate product starting from 20th element of factorials table
      double val=ft[20];
      for(uint i=21; i<=n; ++i)
         val*=i;
      //--- return result
      return(val);
     }
  }
//+------------------------------------------------------------------+
//| MathCorrelationSpearman's RHO                                    |
//+------------------------------------------------------------------+
template <typename T>
double spearmanRho(vector<T> &vec1,vector<T> &vec2)
  {
   ulong size=vec1.Size();
   if(size<1 || vec2.Size()!=size)
     {
      Print(__FUNCTION__, " vector lengths donot match ");
      return(EMPTY_VALUE);
     }
//--- calculate ranks
   vector rank_x = rankVector(vec1);
   vector rank_y = rankVector(vec2);
//---
   return rank_x.CorrCoef(rank_y);
  }
//+------------------------------------------------------------------+
//| MathCorrelationKendall                                           |
//+------------------------------------------------------------------+
template <typename T>
double kendalTau(const vector<T> &vec1,const vector<T> &vec2)
  {
   double tau = double("nan");
   ulong size=vec1.Size();
   if(size==0 || vec2.Size()!=size)
     {
      Print(__FUNCTION__, " size of input vectors donot match ");
      return(tau);
     }
//---
   long cnt1=0,cnt2=0,cnt=0;
//---
   for(long i=0; i<long(size); ++i)
     {
      for(long j=i+1; j<long(size); ++j)
        {
         T delta1=vec1[i]-vec1[j];
         T delta2=vec2[i]-vec2[j];
         T delta=delta1*delta2;
         if(delta==0)
           {
            if(delta1!=0)
               ++cnt1;
            if(delta2!=0)
               ++cnt2;
           }
         else
           {
            ++cnt1;
            ++cnt2;
            if(delta>0.0)
               ++cnt;
            else
               cnt--;
           }
        }
     }
//--- calculate Kendall tau
   long den=cnt1*cnt2;
   if(den==0)
     {
      Print(__FUNCTION__, " failed zero check at line 3139 ");
      return tau;
     }
   tau=double(cnt)/MathSqrt(den);
//---
   return(tau);
  }
//+------------------------------------------------------------------+
//| MathRank                                                          |
//+-------------------------------------------------------------------+
template<typename T>
vector<T> rankVector(vector<T> &in_vec)
  {
   vector<T> rank = vector<T>::Ones(in_vec.Size());
   ulong size=in_vec.Size();
   if(size<1)
      return rank;
   if(size==1)
     {
      rank[0]=1;
      return rank;
     }
//--- prepare arrays
   vector<T>values = in_vec;
   ulong indices[];
//---
   if(!arange(indices,int(size)))
      return rank;
//---
   ulong i,j,k,t,tmpi;
   double tmp;
//--- sort
   if(size!=1)
     {
      i=2;
      do
        {
         t=i;
         while(t!=1)
           {
            k=t/2;
            if(values[k-1]>=values[t-1])
               t=1;
            else
              {
               //--- swap
               tmp=values[k-1];
               values[k-1]=values[t-1];
               values[t-1]=tmp;
               tmpi=indices[k-1];
               indices[k-1]=indices[t-1];
               indices[t-1]=tmpi;
               t=k;
              }
           }
         i=i+1;
        }
      while(i<=size);
      i=size-1;
      do
        {
         //--- swap
         tmp=values[i];
         values[i]=values[0];
         values[0]=tmp;
         tmpi=indices[i];
         indices[i]=indices[0];
         indices[0]=tmpi;
         t=1;
         while(t!=0)
           {
            k=2*t;
            if(k>i)
               t=0;
            else
              {
               if(k<i)
                  if(values[k]>values[k-1])
                     ++k;
               if(values[t-1]>=values[k-1])
                  t=0;
               else
                 {
                  //--- swap
                  tmp=values[k-1];
                  values[k-1]=values[t-1];
                  values[t-1]=tmp;
                  tmpi=indices[k-1];
                  indices[k-1]=indices[t-1];
                  indices[t-1]=tmpi;
                  t=k;
                 }
              }
           }
         i=i-1;
        }
      while(i>=1);
     }
//--- compute tied ranks
   i=0;
   while(i<size)
     {
      j=i+1;
      while(j<size)
        {
         //--- check
         if(values[j]!=values[i])
            break;
         j=j+1;
        }
      for(k=i; k<j; ++k)
         values[k]=1+(i+j-1)*0.5;
      i=j;
     }
//--- set output values
   for(i=0; i<size; ++i)
      rank[indices[i]]=values[i];
//---

   return rank;
  }
//+------------------------------------------------------------------+
//| masking vector based on < condition                              |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsLt(vector<T>&in, const T value)
  {
   vector<T>out = vector<T>::Zeros(in.Size());

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(in[i]<value)
         out[i] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking vector based on <= condition                              |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsLte(vector<T>&in, const T value)
  {
   vector<T>out = vector<T>::Zeros(in.Size());

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(in[i]<=value)
         out[i] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking vector based on > condition                              |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsGt(vector<T>&in, const T value)
  {
   vector<T>out = vector<T>::Zeros(in.Size());

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(in[i]>value)
         out[i] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking vector based on >= condition                              |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsGte(vector<T>&in, const T value)
  {
   vector<T>out = vector<T>::Zeros(in.Size());

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(in[i]>=value)
         out[i] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking matrix based on < condition                              |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsLt(matrix<T>&in, const T value)
  {
   matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(in[i][j]<value)
            out[i][j] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking matrix based on <= condition                              |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsLte(matrix<T>&in, const T value)
  {
   matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(in[i][j]<=value)
            out[i][j] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking matrix based on > condition                              |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsGt(matrix<T>&in, const T value)
  {
   matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(in[i][j]>value)
            out[i][j] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking matrix based on >= condition                              |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsGte(matrix<T>&in, const T value)
  {
   matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(in[i][j]>=value)
            out[i][j] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking vector based on == condition by digits                   |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereVectorIsEq(vector<T>&in, const T value, int _digits = 8)
  {
   vector<T>out = vector<T>::Zeros(in.Size());

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(NormalizeDouble(in[i],_digits)==value)
         out[i] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| masking matrix based on == condition by digits                   |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> whereMatrixIsEq(matrix<T>&in, const T value, int _digits = 8)
  {
   matrix<T>out = matrix<T>::Zeros(in.Rows(),in.Cols());

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(NormalizeDouble(in[i][j],_digits)>=value)
            out[i][j] = 1.0;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| replace vector members that are not numbers                      |
//+------------------------------------------------------------------+
template<typename T>
vector<T> replace_whereVectorIsNaN(vector<T>&in,const T fill = 0.0)
  {
   vector<T>out = in;

   if(in.HasNan()==in.Size())
     {
      out.Fill(fill);
      return out;
     }

   for(ulong i = 0; i<in.Size(); ++i)
     {
      if(MathClassify(in[i])!=FP_NORMAL)
         out[i] = fill;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| replace matrix members that are not numbers                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> replace_whereMatrixIsNaN(matrix<T>&in,const T fill = 0.0)
  {
   matrix<T>out = in;

   if(in.HasNan()==(in.Rows()*in.Cols()))
     {
      out.Fill(fill);
      return out;
     }

   for(ulong i = 0; i<in.Rows(); ++i)
     {
      for(ulong j = 0; j<in.Cols(); ++j)
         if(MathClassify(in[i][j])==FP_NAN)
            out[i][j] = fill;
     }

   return out;
  }
//+------------------------------------------------------------------+
//| Copy to select columns                                           |
//+------------------------------------------------------------------+
template<typename T>
bool copyToSelectCols(matrix<T> &copyto, matrix<T> &from, vector &v_mask)
  {
   if(from.Cols()!=copyto.Cols() || v_mask.Size()!=copyto.Cols())
     {
      Print(__FUNCTION__, " invalid parameters ");
      return false;
     }

   for(ulong i = 0; i<copyto.Cols(); ++i)
     {
      if(v_mask[i])
        {
         if(!copyto.Col(from.Col(i),i))
           {
            Print(__FUNCTION__, " failed column assignment ", GetLastError());
            return false;
           }
        }
     }
   return true;
  }
//+------------------------------------------------------------------+
//| Copy to select rows                                              |
//+------------------------------------------------------------------+
template<typename T>
bool copyToSelectRows(matrix<T> &copyto, matrix<T> from, vector &v_mask)
  {
   if(from.Rows()!=copyto.Rows() || v_mask.Size()!=copyto.Rows())
     {
      Print(__FUNCTION__, " invalid parameters ");
      return false;
     }

   for(ulong i = 0; i<copyto.Rows(); ++i)
     {
      if(v_mask[i])
        {
         if(!copyto.Row(from.Row(i),i))
           {
            Print(__FUNCTION__, " failed column assignment ", GetLastError());
            return false;
           }
        }
     }
   return true;
  }
//+------------------------------------------------------------------+
//|Copy select elements                                              |
//+------------------------------------------------------------------+
template<typename T>
bool copySelectElements(vector<T> &copyto, vector<T> &from, vector &v_mask)
  {
   if(from.Size()!=copyto.Size() || v_mask.Size()!=copyto.Size())
     {
      Print(__FUNCTION__, " invalid parameters ");
      return false;
     }

   for(ulong i = 0; i<copyto.Size(); ++i)
     {
      if(v_mask[i])
         copyto[i] = from[i];
     }
   return true;
  }
//+------------------------------------------------------------------+
//|   sign                                                           |
//+------------------------------------------------------------------+
template<typename T>
T sign(const T a)
  {
   if(a<T(0))
      return T(-1);
   else
      if(a>T(0))
         return T(1);
      else
         return T(0);
  }
//+------------------------------------------------------------------+
//| calculation of raw moment about specified center                 |
//+------------------------------------------------------------------+
double moment(vector& a, ulong order, double mean=NULL)
  {
   if(!a.Size())
      return a.Mean();
   if(order == 0 || (order == 1 && !mean))
      return order==0?1.0:0.;
   ulong list[];
   ArrayResize(list,list.Size()+1,100);
   list[list.Size()-1]=order;
   ulong current = order;
   while(current>2)
     {
      if(MathMod(current,2))
         current=(current-1)/2;
      else
         current/=2;

      ArrayResize(list,list.Size()+1,100);
      list[list.Size()-1]= current;
     }

   double mn = (mean == 0.0)?a.Mean():mean;

   vector azeromean = a - mn;
   vector s;
   if(list[list.Size()-1] == 1)
      s = azeromean;
   else
      s = pow(azeromean,2.);

   for(int i = int(list.Size())-2; i>=0; --i)
     {
      s = pow(s,2.0);
      if(MathMod(i,2))
         s*=azeromean;
     }
   return s.Mean();
  }
//+------------------------------------------------------------------+
//|  sin(pi*x)                                                       |
//+------------------------------------------------------------------+
template <typename T>
T sinpi(T x)
  {
   T s = 1.0;

   if(x < 0.0)
     {
      x = -x;
      s = -1.0;
     }

   T r = fmod(x, 2.0);
   if(r < 0.5)
     {
      return s * sin(M_PI * r);
     }
   else
      if(r > 1.5)
        {
         return s * sin(M_PI * (r - 2.0));
        }
      else
        {
         return -s * sin(M_PI * (r - 1.0));
        }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
double lgam_large_x(double x)
  {
   const double LS2PI = 0.91893853320467274178;
   double q = (x - 0.5) * log(x) - x + LS2PI;
   if(x > 1.0e8)
     {
      return (q);
     }
   double p = 1.0 / (x * x);
   p = ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x;
   return q + p;
  }
//+------------------------------------------------------------------+
//| helper                                                           |
//+------------------------------------------------------------------+
double polevl(double x, double &coef[], int N)
  {
   double ans;
   int i;
   int p;

   p = 0;
   ans = coef[p++];
   i = N;

   do
     {
      ans = ans * x + coef[p++];
     }
   while(--i != 0);

   return (ans);
  }
//+------------------------------------------------------------------+
//|   helper                                                         |
//+------------------------------------------------------------------+
double p1evl(double x, double &coef[], int N)
  {
   double ans;
   uint p;
   int i;

   p = 0;
   ans = x + coef[p++];
   i = N - 1;

   do
      ans = ans * x + coef[p++];
   while(--i != 0);

   return (ans);
  }
//+------------------------------------------------------------------+
//|  sign of log gamma                                               |
//+------------------------------------------------------------------+
double lgam_sgn(double x, int &sign)
  {
   double p, q, u, w, z;
   int i;

   const double LS2PI = 0.91893853320467274178;
   const double MAXLGM = 2.556348e305;

   double gamma_A[] =
     {
      8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4,
      -2.77777777730099687205E-3, 8.33333333333331927722E-2
     };

   double gamma_B[] =
     {
      -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5,
         -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5
        };

   double gamma_C[] =
     {
      -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5,
         -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6
        };

   sign = 1;

   if(MathClassify(x)!=FP_NORMAL)
     {
      return x;
     }

   if(x < -34.0)
     {
      q = -x;
      w = lgam_sgn(q, sign);
      p = floor(q);
      if(p == q)
        {
         Print(__FUNCTION__, " singular error ");
         return (double("inf"));
        }
      i = (int)p;
      if((i & 1) == 0)
        {
         sign = -1;
        }
      else
        {
         sign = 1;
        }
      z = q - p;
      if(z > 0.5)
        {
         p += 1.0;
         z = p - q;
        }
      z = q * sinpi(z);
      if(z == 0.0)
        {
         Print(__FUNCTION__, " singular error ");
         return (double("inf"));
        }
      z = log(M_PI) - log(z) - w;
      return (z);
     }

   if(x < 13.0)
     {
      z = 1.0;
      p = 0.0;
      u = x;
      while(u >= 3.0)
        {
         p -= 1.0;
         u = x + p;
         z *= u;
        }
      while(u < 2.0)
        {
         if(u == 0.0)
           {
            Print(__FUNCTION__, " singular error ");
            return (double("inf"));
           }
         z /= u;
         p += 1.0;
         u = x + p;
        }
      if(z < 0.0)
        {
         sign = -1;
         z = -z;
        }
      else
        {
         sign = 1;
        }
      if(u == 2.0)
        {
         return (log(z));
        }
      p -= 2.0;
      x = x + p;
      p = x * polevl(x, gamma_B, 5) / p1evl(x, gamma_C, 6);
      return (log(z) + p);
     }
   if(x>MAXLGM)
      return sign*double("inf");
   if(x>=1000.)
      return lgam_large_x(x);
   q = (x-0.5)*log(x) - x +LS2PI;
   p = 1./(x*x);
   return q+polevl(p,gamma_A,4)/x;
  }
//+-------------------------------------------------------------------------+
//| Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1)|
//+-------------------------------------------------------------------------+
double lbeta_asymp(double a, double b, int &sgn)
  {
   double r = lgam_sgn(b, sgn);
   r -= b * log(a);

   r += b * (1 - b) / (2 * a);
   r += b * (1 - b) * (1 - 2 * b) / (12 * a * a);
   r += -b * b * (1 - b) * (1 - b) / (12 * a * a * a);

   return r;
  }
//+------------------------------------------------------------------+
//| Special case for a negative integer argument                     |
//+------------------------------------------------------------------+
double beta_negint(int a, double b)
  {
   int sgn;
   if(b == int(b) && 1 - a - b > 0)
     {
      sgn = (MathMod(b, 2) == 0) ? 1 : -1;
      return sgn * CBetaF::Beta(1. - a - b, b);
     }
   else
     {
      Print(__FUNCTION__, " overflow error ");
      return double("inf");
     }
  }
//+------------------------------------------------------------------+
//| helper                                                           |
//+------------------------------------------------------------------+
double lbeta_negint(int a, double b)
  {
   double r;
   if(b == int(b) && 1 - a - b > 0)
     {
      r = lbeta(1 - a - b, b);
      return r;
     }
   else
     {
      Print(__FUNCTION__, " overflow error ");
      return double("inf");
     }
  }
//+------------------------------------------------------------------+
//|  helper                                                          |
//+------------------------------------------------------------------+
double lbeta(double a, double b)
  {
   double y;
   int sign;
   const double MAXGAM = 171.624376956302725;
   const double beta_ASYMP_FACTOR = 1e6;
   sign = 1;

   if(a <= 0.0)
     {
      if(a == floor(a))
        {
         if(a == int(a))
           {
            return lbeta_negint(int(a), b);
           }
         else
           {
            Print(__FUNCTION__, " overflow error ");
            return (sign * double("inf"));
           }
        }
     }

   if(b <= 0.0)
     {
      if(b == floor(b))
        {
         if(b == int(b))
           {
            return lbeta_negint(int(b), a);
           }
         else
           {
            Print(__FUNCTION__, " overflow error ");
            return (sign * double("inf"));
           }
        }
     }

   if(MathAbs(a) < MathAbs(b))
     {
      y = a;
      a = b;
      b = y;
     }

   if(MathAbs(a) > beta_ASYMP_FACTOR * MathAbs(b) && a > beta_ASYMP_FACTOR)
     {
      /* Avoid loss of precision in lgam(a + b) - lgam(a) */
      y = lbeta_asymp(a, b, sign);
      return y;
     }

   y = a + b;
   if(MathAbs(y) > MAXGAM || MathAbs(a) > MAXGAM || MathAbs(b) > MAXGAM)
     {
      int sgngam;
      y = lgam_sgn(y, sgngam);
      sign *= sgngam; /* keep track of the sign */
      y = lgam_sgn(b, sgngam) - y;
      sign *= sgngam;
      y = lgam_sgn(a, sgngam) + y;
      sign *= sgngam;
      return (y);
     }

   y = CGammaFunc::GammaFunc(y);
   a = CGammaFunc::GammaFunc(a);
   b = CGammaFunc::GammaFunc(b);
   if(y == 0.0)
     {
      Print(__FUNCTION__, " overflow error ");
      return (sign * double("inf"));
     }

   if(MathAbs(MathAbs(a) - MathAbs(y)) > MathAbs(MathAbs(b) - MathAbs(y)))
     {
      y = b / y;
      y *= a;
     }
   else
     {
      y = a / y;
      y *= b;
     }

   if(y < 0)
     {
      y = -y;
     }

   return (log(y));
  }
//+------------------------------------------------------------------+
//| Binomial coefficient                                             |
//+------------------------------------------------------------------+
double binom(double n, double k)
  {
   double kx,nx,num, den, dk, sgn;
   if(n < 0)
     {
      nx = floor(n);
      if(n == nx)
        {
         // Undefined
         return double("nan");
        }
     }

   kx = floor(k);
   if(k == kx && (MathAbs(n) > 1e-8 || n == 0))
     {
      nx = floor(n);
      if(nx == n && kx > nx / 2 && nx > 0)
        {
         // Reduce kx by symmetry
         kx = nx - kx;
        }

      if(kx >= 0 && kx < 20)
        {
         num = 1.0;
         den = 1.0;
         for(int i = 1; i < 1 + int(kx); i++)
           {
            num *= i + n - kx;
            den *= i;
            if(MathAbs(num) > 1e50)
              {
               num /= den;
               den = 1.0;
              }
           }
         return num / den;
        }
     }

// general case
   if(n >= 1e10 * k && k > 0)
     {
      // avoid under/overflows intermediate results
      return exp(-lbeta(1 + n - k, 1 + k) - log(n + 1));
     }
   if(k > 1e8 * MathAbs(n))
     {
      // avoid loss of precision
      num = CGammaFunc::GammaFunc(1 + n) / MathAbs(k) + CGammaFunc::GammaFunc(1 + n) * n / (2 * k * k); // + ...
      num /= M_PI * pow(MathAbs(k), n);
      if(k > 0)
        {
         kx = floor(k);
         if(int(kx) == kx)
           {
            dk = k - kx;
            sgn = (MathMod(int(kx), 2) == 0) ? 1 : -1;
           }
         else
           {
            dk = k;
            sgn = 1;
           }
         return num * sin((dk - n) * M_PI) * sgn;
        }
      kx = floor(k);
      if(int(kx) == kx)
        {
         return 0;
        }
      return num * sin(k * M_PI);
     }
   return 1 / (n + 1) / CBetaF::Beta(1 + n - k, 1 + k);
  }
//+------------------------------------------------------------------+
//|  N choose K                                                      |
//+------------------------------------------------------------------+

template<typename T>
T comb(ulong N, T k, bool exact=false, bool repetition = false)
  {
   if(repetition)
      return comb(N+ulong(k)-1,(T)k,exact);
   if(exact)
      return comb(N,k);
   else
     {
      bool cond = false;
      cond = ((k<=T(N)) && (N>=0) && (k>=0));
      T val = (T)binom(double(N),double(k));
      if(!cond)
         val = 0;
      return val;
     }
  }

//+------------------------------------------------------------------+
//| Compute the kurtosis (Fisher or Pearson) of a dataset.           |
//+------------------------------------------------------------------+
double kurtosis(vector& in,bool fisher=true, bool bias=true)
  {
   ulong n = in.Size();
   double out;
   double mean = in.Mean();
   double m2 = moment(in,2,mean);
   double m4 = moment(in,4,mean);

   bool zero = (m2<= pow(mean*2.e-16,2.));
   out = (!zero)?m4/pow(m2,2.):double("nan");

   if(!bias)
     {
      bool cancorrect = (~zero&(n>3));
      double nval = 1.0/(n-2)/(n-3) * ((pow(n,2)-1.0)*m4/pow(m2,2.0) - 3*pow((n-1),2.0));
      if(cancorrect)
         out = nval + 3.;
      else
         out = double("nan");
     }

   if(fisher)
      out = out - 3.;
   return out;
  }
//+------------------------------------------------------------------+
//| get element wise maximum                                         |
//+------------------------------------------------------------------+
vector maximum(vector &x, vector &y)
  {
   if(x.Size()!=y.Size())
     {
      Print(__FUNCTION__, " vector size mismatch ");
      return vector::Zeros(0);
     }

   matrix xy(x.Size(),2);
   if(!xy.Col(x,0) || !xy.Col(y,1))
     {
      Print(__FUNCTION__, " column insertion error ", GetLastError());
      return vector::Zeros(0);
     }

   return xy.Max(1);
  }
//+------------------------------------------------------------------+
//| get element wise minimum                                         |
//+------------------------------------------------------------------+
vector minimum(vector &x, vector &y)
  {
   if(x.Size()!=y.Size())
     {
      Print(__FUNCTION__, " vector size mismatch ");
      return vector::Zeros(0);
     }

   matrix xy(x.Size(),2);
   if(!xy.Col(x,0) || !xy.Col(y,1))
     {
      Print(__FUNCTION__, " column insertion error ", GetLastError());
      return vector::Zeros(0);
     }

   return xy.Min(1);
  }
//+------------------------------------------------------------------+
//| get element wise maximum                                         |
//+------------------------------------------------------------------+
vector maximum(vector &x, double v)
  {
   vector y(x.Size());
   y.Fill(v);

   matrix xy(x.Size(),2);
   if(!xy.Col(x,0) || !xy.Col(y,1))
     {
      Print(__FUNCTION__, " column insertion error ", GetLastError());
      return vector::Zeros(0);
     }

   return xy.Max(1);
  }
//+------------------------------------------------------------------+
//| get element wise minimum                                         |
//+------------------------------------------------------------------+
vector minimum(vector &x, double v)
  {
   vector y(x.Size());
   y.Fill(v);

   matrix xy(x.Size(),2);
   if(!xy.Col(x,0) || !xy.Col(y,1))
     {
      Print(__FUNCTION__, " column insertion error ", GetLastError());
      return vector::Zeros(0);
     }

   return xy.Min(1);
  }
//+------------------------------------------------------------------+
//| return vector with elements that meet certain condition          |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereIsNan(vector<T>& in)
  {
   ulong any = in.HasNan();
   if(any)
     {
      vector out = vector::Zeros(any)
                   ulong count = 0;
      for(ulong i = 0; i<in.Size(); ++i)
         if(MathClassify(in[i]) == FP_NAN)
            out[count++] = double(i);
      return out;
     }
   else
      return vector::Zeros(0);
  }
//+------------------------------------------------------------------+
//| return vector with elements that meet certain condition          |
//+------------------------------------------------------------------+
template<typename T>
vector<T> whereIsNotNan(vector<T>& in)
  {
   ulong any = in.HasNan();
   if(any)
     {
      vector out = vector::Zeros(in.Size() - any);
      ulong count = 0;
      for(ulong i = 0; i<in.Size(); ++i)
         if(MathClassify(in[i]) != FP_NAN)
            out[count++] = in[i];
      return out;
     }
   else
      return in;
  }
//+------------------------------------------------------------------+
//| matrix rank                                                      |
//+------------------------------------------------------------------+
ulong matrix_rank(matrix& in)
  {
   matrix A = in;
   matrix u,v;
   vector s;
   if(!A.SVD(u,v,s))
     {
      Print(__FUNCTION__, " svd failure ", GetLastError());
      return 0;
     }

   double rtol = MathMax(A.Rows(),A.Cols())*DBL_EPSILON;
   double tol = s.Max()*rtol;

   vector count = whereVectorIsGt(s,tol);
   return ulong(count.Sum());
  }
//+------------------------------------------------------------------+
//| toeplitz matrix construction                                     |
//+------------------------------------------------------------------+
template<typename T>
matrix<T> toeplitz(vector<T>& x)
  {
   ulong n = x.Size();
   matrix<T> out(n,n);
   for(ulong i = 0; i<n; ++i)
      for(ulong j = 0; j<n; ++j)
         if(i<=j)
            out[i,j] = x[j-i];
         else
            out[i,j] = x[i-j];
   return out;
  }
//+------------------------------------------------------------------+
//| Compute the (Moore-Penrose) pseudo-inverse of a matrix.          |
//+------------------------------------------------------------------+
matrix pinv(matrix& a)
  {
   double rcond = 1.e-15;

   vector vrcond(a.Cols());
   vrcond.Fill(rcond);

   if(!a.Rows())
      return a;

   matrix aa = a.Conjugate();

   matrix u,v;
   vector s;
   if(!aa.SVD(u,v,s))
     {
      Print(__FUNCTION__," svd failed ",GetLastError());
      return matrix::Zeros(0,0);
     }

   CMatrixDouble m = aa;
   double singular_values[];
   CMatrixDouble U,Vt;

   if(!CSingValueDecompose::RMatrixSVD(m,m.Rows(),m.Cols(),2,2,2,singular_values,U,Vt))
     {
      Print(__FUNCTION__," error ", GetLastError());
      return matrix::Zeros(0,0);;
     }
   else
     {
      s.Assign(singular_values);
      u = U.ToMatrix();
      u = sliceMatrixCols(u,0,long(s.Size()));
      v = Vt.ToMatrix();
     }

   vector cuttoff = vrcond*s.Max();
   for(ulong i = 0; i<s.Size(); ++i)
      if(s[i]>cuttoff[i])
         s[i] = 1./s[i];
      else
         s[i] = 0.0;

   matrix mats = matrix::Zeros(s.Size(),1);
   mats.Col(s,0);
   matrix vm = multiply(u.Transpose(),s);
   return v.Transpose().MatMul(vm);
  }
//+------------------------------------------------------------------+
//|  shuffle elements of vector                                      |
//+------------------------------------------------------------------+
template<typename T>
vector  shuffleVector(vector<T> &in,CHighQualityRandStateShell* rstate=NULL)
  {
   vector<T>out(in.Size());
   bool destroy = false;
//---
   if(rstate == NULL)
     {
      destroy = true;
      rstate = new CHighQualityRandStateShell();
      CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
     }
//---
   ulong k=0;
   ulong size=in.Size();
   out=in;
   int error = 0;
   T tmp;
//---
   for(ulong i=0; i<size; ++i)
     {
      k=(int)(CHighQualityRand::HQRndUniformR(rstate.GetInnerObj())*size);
      //---
      if(k>=size)
         k=size-1;
      //---
      tmp = out[i];
      out[i]=out[k];
      out[k]=tmp;
     }
//---
   if(destroy)
      delete rstate;
//---
   return out;
  }
//+------------------------------------------------------------------+
//|  shuffle elements of vector                                      |
//+------------------------------------------------------------------+
template<typename T>
bool  shuffleArray(T &in[],CHighQualityRandStateShell* rstate=NULL)
  {
   bool destroy = false;
//---
   if(rstate == NULL)
     {
      destroy = true;
      rstate = new CHighQualityRandStateShell();
      CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
     }
//---
   ulong k=0;
   uint size=in.Size();
   int error = 0;
   T tmp;
//---
   for(ulong i=0; i<size; ++i)
     {
      k=(int)(CHighQualityRand::HQRndUniformR(rstate)*size);
      //---
      if(k>=size)
         k=size-1;
      //---
      tmp = in[i];
      in[i]=in[k];
      in[k]=tmp;
     }
//---
   if(destroy)
      delete rstate;
//---
   return true;
  }
//+------------------------------------------------------------------+
//|  shuffle elements of vector                                      |
//+------------------------------------------------------------------+
template<typename T>
matrix<T>  shuffleMatrix(matrix<T> &in,bool by_rows=true,CHighQualityRandStateShell* rstate=NULL)
  {
   matrix<T>out(in.Rows(), in.Cols());
   bool destroy = false;
//---
   if(rstate == NULL)
     {
      destroy = true;
      rstate = new CHighQualityRandStateShell();
      CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());
     }
//---
   ulong size=by_rows?in.Rows():in.Cols();
   vector<T> vec;
   out=in;
//---
   for(ulong i=0; i<size; ++i)
     {
      vec = by_rows?in.Row(i):in.Col(i);
      vec = shuffleVector(vec,rstate);
      if((by_rows && !out.Row(vec,i)) ||
         (!by_rows && !out.Col(vec,i)))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return matrix::Zeros(size,size);
        }
     }
//---
   if(destroy)
      delete rstate;
//---
   return out;
  }
//+------------------------------------------------------------------+
//| folds structure                                                  |
//+------------------------------------------------------------------+

struct Folds
  {
   vector            test_indices;
   vector            train_indices;

                     Folds(void)
     {
      test_indices = train_indices = vector::Zeros(0);
     }

                     Folds(Folds& other)
     {
      test_indices = other.test_indices;
      train_indices = other.train_indices;
     }
   void              operator=(Folds& other)
     {
      test_indices = other.test_indices;
      train_indices = other.train_indices;
     }
  };

//+------------------------------------------------------------------+
//|KFold cross validation: returns the indices for the test/train    |
//|split.                                                            |
//+------------------------------------------------------------------+
bool kfold(ulong data_length, Folds& folds[],ulong n_splits = 2, bool shuffle = true, int random_seed = 0)
  {
   if(n_splits==0 || data_length<=n_splits)
     {
      Print(__FUNCTION__, " invalid inputs ");
      return false;
     }

   vector indices = arange(data_length);

   if(!indices.Size())
      return false;

   if(shuffle)
     {
      CHighQualityRandStateShell rstate;
      if(random_seed)
         CHighQualityRand::HQRndSeed(random_seed,random_seed+1,rstate.GetInnerObj());
      else
         CHighQualityRand::HQRndRandomize(rstate.GetInnerObj());

      indices = shuffleVector(indices,(CHighQualityRandStateShell*)GetPointer(rstate));
     }

   long start,stop,current,size;
   vector fold_sizes = vector::Zeros(n_splits);
   ArrayResize(folds,(int)n_splits);
   fold_sizes.Fill(floor(data_length/n_splits));
   ulong to = ulong(MathMod(data_length,n_splits));
   vector temp;
   for(ulong i = 0; i<to; ++i)
      fold_sizes[i]+=1.;

   current = start = stop = size = 0;

   for(ulong i = 0; i<fold_sizes.Size(); ++i)
     {
      start = current;
      size = long(fold_sizes[i]);
      stop = current + size;
      folds[i].test_indices = sliceVector(indices,start,stop);
      size = start+(long(data_length)-stop);
      if(size)
        {
         folds[i].train_indices = vector::Zeros(size);
         if(start)
           {
            temp = sliceVector(indices,0,start);
            if(!vectorCopy(folds[i].train_indices,temp,0,start))
               return false;
           }
         if(stop<long(data_length))
           {
            temp = sliceVector(indices,stop);
            if(!vectorCopy(folds[i].train_indices,temp,start))
               return false;
           }
        }
      current = stop;
     }

   return true;
  }
//+------------------------------------------------------------------+
//| expanding windows                                                |
//+------------------------------------------------------------------+
bool exp_windows(ulong data_length, Folds& folds[],ulong n_splits = 2)
  {
   vector indices = arange(data_length);

   long test_size = long(data_length/(n_splits+1));

   ArrayResize(folds,int(n_splits));

   long test_end, test_start;
   for(ulong i = 0; i<n_splits; ++i)
     {
      test_end = long(data_length)-long(i*test_size);
      test_start = test_end - test_size;

      folds[n_splits-i-1].train_indices=sliceVector(indices,0,test_start);
      folds[n_splits-i-1].test_indices=sliceVector(indices,test_start,test_end);
     }

   return true;

  }
//+------------------------------------------------------------------+
//| The Rosenbrock function.                                          |
//+-------------------------------------------------------------------+
double rosen(vector &x)
  {
   vector x1 = sliceVector(x,1);
   vector x2 = sliceVector(x,0,long(x.Size()-1));
   vector r = 100.0*pow(x1-pow(x2,2.0),2.0)+pow(1.-x2,2.0);
   return r.Sum();
  }
//+------------------------------------------------------------------+
//|Gradient of the Rosenbrock function                               |
//+------------------------------------------------------------------+
vector rosen_gradient(vector &x)
  {
   vector grad = vector::Zeros(x.Size());
   if(x.Size()>2)
     {
      vector xm = sliceVector(x,1,long(x.Size()-1));
      vector xm_m1 = sliceVector(x,0,long(x.Size()-2));
      vector xm_p1 = sliceVector(x,2);
      vector temp = (200.0 * (xm - pow(xm_m1,2.0)) - 400.0 * (xm_p1 - pow(xm,2.0)) * xm - 2.0 * (1.0 - xm));
      if(!vectorCopy(grad,temp,1,long(grad.Size()-1)))
         return vector::Zeros(0);
     }
   grad[0] = -400.0 * x[0] * (x[1] - pow(x[0],2.0)) - 2.0 * (1.0 - x[0]);
   grad[grad.Size()-1] = 200.0 * (x[x.Size()-1] - pow(x[x.Size()-2],2.0));
   return grad;
  }
//+------------------------------------------------------------------+
//| Hessian matrix of the Rosenbrock function.                       |
//+------------------------------------------------------------------+
matrix rosen_hess(vector &x)
  {
   vector x1 = sliceVector(x,0,long(x.Size()-1));
   vector x2 = sliceVector(x,1,long(x.Size()-1));
   vector x3 = vector::Zeros(x.Size());
   matrix h1,h2,h3,H;
   if(!h1.Diag(-400.0*x1,1) || !h2.Diag(400.0*x1,-1))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return matrix::Zeros(0,0);
     }

   H = h1 - h2;

   x3[0] = 1200.0 * pow(x[0],2.0) - 400.0 * x[1] + 2.0;
   vector temp = 202.0 + 1200.0 * pow(x2,2) - 400.0 * sliceVector(x,2);
   x3[x3.Size()-1] = 200.0;

   if(!vectorCopy(x3,temp,1,long(x3.Size()-1)))
      return matrix::Zeros(0,0);

   if(!h3.Diag(x3))
     {
      Print(__FUNCTION__, " error ", GetLastError());
      return matrix::Zeros(0,0);
     }

   return H+h3;
  }
}

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