//+------------------------------------------------------------------+
//|                                                          SSA.mq4 |
//|                      Copyright © 2007, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright © 2007, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
#property library
 
//+------------------------------------------------------------------+
//| My function                                                      |
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
int  MaxCompCol=20; 
int  MaxCatLag=200;
int  MaxLen=5000;
 
//---
/*class Matrix{
public:
  int i; 
  int j; //выап
  double **m;
  Matrix(int i, int j);
  ~Matrix();
};
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
Matrix::Matrix(int ni, int nj)
{
 m=new double*[ni];
 for(int k=0;k<ni;k++) m[k]=new double[nj];
 i=ni; j=nj;
}
Matrix::~Matrix()
{
  for(int k=0;k<i;k++) delete[] m[k];
  delete[] m;
}*/
void Equal(double &x1[][], double x2[][], int t=0) //(Matrix *x1, Matrix *x2, int t=0)
{
 int N=ArrayRange(x1,0);
 int M=ArrayRange(x1,1);
 for(int i=0;i<N;i++) //(int i=0;i<x1->i;i++)
   for(int j=0;j<M;j++){ //(int j=0;j<x1->j;j++)
     if(t!=0) x1[i][j]=x2[j][i]; //x1->m[i][j]=x2->m[j][i];
     else x1[i][j]=x2[i][j]; //x1->m[i][j]=x2->m[i][j];
   } 
}
//---------------------------------------------------------------------------
/*Вычисление функции Sn, нужна для вычисления собственных значений
Там считатются отрицательные определители*/
double gaussSn(double A[][],double l,int n) //(Matrix *A,double l,int n)
{
  int count,cp,i1,i,j,k;
  double B[200][200];//Matrix *B = new Matrix(MaxCatLen,MaxCatLen);
  double c;
  double w[]; ArrayResize(w,n);//double *w=new double[n]; //double w[MaxLen-MaxCatLen];
  double s1,s2;
  //------
  Equal(B,A);
  //------
  for (i=0;i<=n-1;i++)
  {
   B[i][i] = B[i][i]-l; //B->m[i][i] = B->m[i][i]-l;
  }
  cp=0;
  for(k=0;k<=n-2;k++)
  {
    for(i=k+1;i<=n-1;i++)
    {
      if(B[k,k]==0) // (B->m[k,k]==0)
      {
        for(i1=0;i1<=n-1;i1++)
        {
          w[i1] = B[i1][k]; //w[i1] = B->m[i1][k];
          B[i1][k] = B[i1][k+1]; //B->m[i1][k] = B->m[i1][k+1];
          B[i1][k+1] = w[i1]; //B->m[i1][k] = B->m[i1][k+1];
        }
        cp=cp+1;
      }
      c=B[i][k]/B[k][k]; //c=1.0*B->m[i][k]/B->m[k][k];
      for (j=0;j<=n-1;j++)
      {
        B[i][j] = B[i][j]-B[k][j]*c;//B->m[i][j] = B->m[i][j]-B->m[k][j]*c;
      }
    }
  }
   count=0;
   s1=1;
   for( i=0;i<=n-1;i++)
   {
      s2 = B[i][i]; //s2 = B->m[i][i];
      if (s2<0) count = count+1;
   }
   //delete B;
   //delete[] w;
   return (1.0*count);
}
//---------------------------------------------------------------------------
/*Вычисление собственных значений методом бисекции}
  Хорош тем что сколько надо собственных чисел столько и посчитает, 
  сильно экономит ресурс*/
double gaussbisectionl(double &A[][],int k,int n) //(Matrix *A,int k,int n)
{
    double e1,maxnorm,cn,a1,b1,c;
    int i,j;
    maxnorm = 0;
    for (i=0;i<=n-1;i++)
    {
        cn = 0;
        for(j=0 ;j<=n-1;j++)
            cn = cn+A[i][j]; //cn = cn+A->m[i][j];
        if (maxnorm<cn) maxnorm=cn;
    }
    a1=0;
    b1=10*maxnorm;
    e1=1.0*maxnorm/10000000;
    while (MathAbs(b1-a1)>e1)
    {
        c = 1.0*(a1+b1)/2;
        if (gaussSn(A,c,n)<k)
        a1=c;
        else
        b1=c;
    }
 return ((a1+b1)/2.0);
}
//---------------------------------------------------------------------------
/*Вычисляет собственные вектора для уже вычисленных собственных чисел l*/
void svector(double &A[][],double l,int n,double &V[]) //(Matrix *A,double l,int n,double *V)
{
    int cp,i1,i,j,k;
    double B[200][200]; //Matrix *B = new Matrix(MaxCatLen,MaxCatLen);
    double c;
    double w[]; ArrayResize(w,n); //double *w=new double[n]; //double w[MaxLen-MaxCatLen];
    Equal(B,A);
    for (i=0;i<=n-1;i++)  B[i][i]=B[i][i]-l; //B->m[i][i]=B->m[i][i]-l;
    cp=0;
    for (k=0;k<=n-2;k++)
        for (i=k+1;i<=n-1;i++)
        {
            if (B[k][k]==0) //(B->m[k][k]==0)
            {
                for (i1=0;i<=n-1;i++)
                {
                    w[i1]=B[i1][k]; //w[i1]=B->m[i1][k];
                    B[i1][k]=B[i1][k+1]; //B->m[i1][k]=B->m[i1][k+1];
                    B[i1][k+1]=w[i1]; //B->m[i1][k+1]=w[i1];
                }
                cp=cp+1;
            }
            c=1.0*B[i][k]/B[k][k]; //c=1.0*B->m[i][k]/B->m[k][k];
            for (j=0;j<=n-1;j++)
            {
                B[i][j]=B[i][j]-B[k][j]*c; //B->m[i][j]=B->m[i][j]-B->m[k][j]*c;
            }
        }
    V[n-1]=1;
    c=1;
    for (i=n-2;i>=0;i--)
    {
        V[i]=0;
        for (j=i;j<=n-1;j++)    V[i]= V[i]-B[i][j]*V[j]; //V[i]= 1.0*V[i]-B->m[i][j]*V[j];
        V[i]=V[i]/B[i][i]; //V[i]=V[i]/B->m[i][i];
        c = c+V[i]*V[i]; //c = c+1.0*V[i]*V[i];
    }
    for (i=0;i<=n-1;i++)    V[i]= V[i]/MathSqrt(c);
    //delete B;
    //delete[] w;
}
/*А вот и сама гусеница}
{Х-вектор исходного ряда
n-его длина
l-длина лага
s-число собственных компонент
(там исходный ряд разбивается на компонеты, а потом восстанавливается, тут ты и задаешь сколько компонент тебе надо)
Y - восстановленный ряд (сглаженный гусеницей) там я в конце еще добавил комментарий, где находятся отдельные
компоненты ряда
*/
//---------------------------------------------------------------------------
void fastsingular( double X[],int n,int l,int s,double &Y[])
{
    if( l>MaxCatLag ) l=MaxCatLag;
   if( s>MaxCompCol ) s=MaxCompCol;
   if( s>l ) s=l;
   if( n>MaxLen ) n=MaxLen;
    double A[200][200]; //Matrix *A = new Matrix(MaxCatLag,MaxCatLag);
    double B[5000][200]; //ArrayResize(B,n); //Matrix *B = new Matrix(n,MaxCatLag);
    double Bn[200][5000]; //Matrix *Bn = new Matrix(MaxCatLag,n);
    
    int kb,lb,m,k,i,j,i1;
    double V[200][5000]; //Matrix *V = new Matrix(MaxCompCol,n);
    double Yn[200][5000]; //Matrix *Yn = new Matrix(MaxCompCol,n);
    double ls[200];
    //---
    double Vtemp[200];
    //---
    j=0;
    k=n-l+1;
    /*Формируем матрицу А в методе который я скачал на сайте создателей эта матрица S*/
    for (i=0;i<=l-1;i++)
    {
        for (j=0;j<=l-1;j++)
        {
            A[i][j]=0; //A->m[i][j]=0;
            for (m=0;m<=k-1;m++) // (m=0;m<=k-1;m++)
            {
                A[i][j]=A[i][j]+X[i+m]*X[m+j]; //A->m[i][j]=A->m[i][j]+1.0*X[i+m]*X[m+j];
                B[m][j]=X[m+j]; //B->m[m][j]=X[m+j];
            }
        }
   }
    /*Находим собственные числа и векторы матрицы А*/
    for (i=0;i<=s-1;i++)
    {
      ls[i]=gaussbisectionl(A,l-i,l);
      svector(A,ls[i],l,Vtemp); //svector(A,ls[i],l,V->m[i]);
      for( j=0;j<=l-1; j++) V[i][j]=Vtemp[j];
    }
   //------
    /*Формируется восстановленная матрица*/
    for (i1=0;i1<=s-1;i1++)
    {
        //------
        for (i=0;i<=k-1;i++)
        {
            Yn[i1][i]=0; //Yn->m[i1][i]=0;
            for (j=0;j<=l-1;j++)  Yn[i1][i]=Yn[i1][i]+B[i][j]*V[i1][j];
                                 //Yn->m[i1][i]=Yn->m[i1][i]+1.0*B->m[i][j]*V->m[i1][j];
        }
        //------
        for (i=0;i<=l-1;i++){
            for (j=0;j<=k-1;j++)    Bn[i][j]=V[i1][i]*Yn[i1][j]; 
                                 //Bn->m[i][j]=1.0*V->m[i1][i]*Yn->m[i1][j];
        }
        //-------
        //Диагональное усреднение (восстановление ряда)
        kb=k;
        lb=l;
        for (i=0;i<=n-1;i++)
        {
            Yn[i1][i]=0; //Yn->m[i1][i]=0;
            if (i<lb-1)
            {
                for (j=0;j<=i;j++)
                {
                    if (l<=k) Yn[i1][i] = Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i] = Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k)  Yn[i1][i] = Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i] = Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*(i+1)); //Yn->m[i1][i]=1.0*Yn->m[i1][i]/(i+1);
            }
            if ((lb-1<=i)&&(i<kb-1))
            {
                for (j=0;j<=lb-1;j++)
                {
                    if (l<=k) Yn[i1][i]=Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k) Yn[i1][i]=Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*lb); //Yn->m[i1][i] = 1.0*Yn->m[i1][i]/lb;
            }
            if (kb-1<=i)
            {
                for (j=i-kb+1;j<=n-kb;j++)
                {
                    if (l<=k) Yn[i1][i]=Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k) Yn[i1][i]=Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*(n-i));//(?! что такое N)//Yn->m[i1][i]=1.0*Yn->m[i1][i]/(n-i);
                //Yn[i1][i]=1.0*Yn[i1][i]/kb;/*(?! что такое N)*/
            }
        }
    }
    /* Вот здесь если не суммировать, то будут отдельные компоненты разложения
    процесса по собственным функциям */
    for (i=0;i<=n-1;i++)
    {
        Y[i]=0;
        //Y[i]=Yn[3][i];
        for (i1=0;i1<=s-1;i1++)  Y[i]=Y[i]+Yn[i1][i]; //Y[i]=Y[i]+Yn->m[i1][i];
    }
    //delete A;
    //delete B;
    //delete Bn;
    //delete V;
    //delete Yn;
    //delete[] X;
}
//-----------------------------------------------------------------------------+
/*{Х-вектор исходного ряда
n-его длина
l-длина лага
s-число собственных компонент
(там исходный ряд разбивается на компонеты, а потом восстанавливается, тут ты и задаешь сколько компонент тебе надо)
Y - восстановленный ряд (сглаженный гусеницей) там я в конце еще добавил комментарий, где находятся отдельные
компоненты ряда
*/
//---------------------------------------------------------------------------
void singularPCA( double X[],int n,int l,int s,double &Y[])
{
    if( l>MaxCatLag ) l=MaxCatLag;
   if( s>MaxCompCol ) s=MaxCompCol;
   if( s>l ) s=l;
   if( n>MaxLen ) n=MaxLen;
    double A[200][200]; //Matrix *A = new Matrix(MaxCatLag,MaxCatLag);
    double B[5000][200]; //ArrayResize(B,n); //Matrix *B = new Matrix(n,MaxCatLag);
    double Bn[200][5000]; //Matrix *Bn = new Matrix(MaxCatLag,n);
    
    int kb,lb,m,k,i,j,i1;
    double V[200][5000]; //Matrix *V = new Matrix(MaxCompCol,n);
    double Yn[200][5000]; //Matrix *Yn = new Matrix(MaxCompCol,n);
    double ls[200];
    //---
    double Vtemp[200];
    //---
    j=0;
    k=n-l+1;
    /*Формируем матрицу А в методе который я скачал на сайте создателей эта матрица S*/
    for (i=0;i<=l-1;i++)
    {
        for (j=0;j<=l-1;j++)
        {
            A[i][j]=0; //A->m[i][j]=0;
            for (m=0;m<=k-1;m++) // (m=0;m<=k-1;m++)
            {
                A[i][j]=A[i][j]+X[i+m]*X[m+j]; //A->m[i][j]=A->m[i][j]+1.0*X[i+m]*X[m+j];
                B[m][j]=X[m+j]; //B->m[m][j]=X[m+j];
            }
        }
   }
    /*Находим собственные числа и векторы матрицы А*/
    //for (i=0;i<=s-1;i++)
   // {
      ls[s]=gaussbisectionl(A,l-s,l);
      svector(A,ls[s],l,Vtemp); //svector(A,ls[i],l,V->m[i]);
      for( j=0;j<=l-1; j++) V[s][j]=Vtemp[j];
   // }
   //------
    /*Формируется восстановленная матрица*/
    i1=s;
    //for (i1=0;i1<=s-1;i1++)
   // {
        //------
        for (i=0;i<=k-1;i++)
        {
            Yn[i1][i]=0; //Yn->m[i1][i]=0;
            for (j=0;j<=l-1;j++)  Yn[i1][i]=Yn[i1][i]+B[i][j]*V[i1][j];
                                 //Yn->m[i1][i]=Yn->m[i1][i]+1.0*B->m[i][j]*V->m[i1][j];
        }
        //------
        for (i=0;i<=l-1;i++){
            for (j=0;j<=k-1;j++)    Bn[i][j]=V[i1][i]*Yn[i1][j]; 
                                 //Bn->m[i][j]=1.0*V->m[i1][i]*Yn->m[i1][j];
        }
        //-------
        //Диагональное усреднение (восстановление ряда)
        kb=k;
        lb=l;
        for (i=0;i<=n-1;i++)
        {
            Yn[i1][i]=0; //Yn->m[i1][i]=0;
            if (i<lb-1)
            {
                for (j=0;j<=i;j++)
                {
                    if (l<=k) Yn[i1][i] = Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i] = Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k)  Yn[i1][i] = Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i] = Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*(i+1)); //Yn->m[i1][i]=1.0*Yn->m[i1][i]/(i+1);
            }
            if ((lb-1<=i)&&(i<kb-1))
            {
                for (j=0;j<=lb-1;j++)
                {
                    if (l<=k) Yn[i1][i]=Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k) Yn[i1][i]=Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*lb); //Yn->m[i1][i] = 1.0*Yn->m[i1][i]/lb;
            }
            if (kb-1<=i)
            {
                for (j=i-kb+1;j<=n-kb;j++)
                {
                    if (l<=k) Yn[i1][i]=Yn[i1][i]+Bn[j][i-j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[j][i-j];
                    if (l>k) Yn[i1][i]=Yn[i1][i]+Bn[i-j][j]; //Yn->m[i1][i]=Yn->m[i1][i]+Bn->m[i-j][j];
                }
                Yn[i1][i] = Yn[i1][i]/(1.0*(n-i));//(?! что такое N)//Yn->m[i1][i]=1.0*Yn->m[i1][i]/(n-i);
                //Yn[i1][i]=1.0*Yn[i1][i]/kb;/*(?! что такое N)*/
            }
        }
    //}
    /* Вот здесь если не суммировать, то будут отдельные компоненты разложения
    процесса по собственным функциям */
    for (i=0;i<=n-1;i++)
    {
        Y[i]=0;
        Y[i]=Yn[s][i];
        //for (i1=0;i1<=s-1;i1++)  Y[i]=Y[i]+Yn[i1][i]; //Y[i]=Y[i]+Yn->m[i1][i];
    }
    //delete A;
    //delete B;
    //delete Bn;
    //delete V;
    //delete Yn;
    //delete[] X;
}
//-----------------------------------------------------------------------------+