/* File: FI.C */

#include "fi.h"
#include <math.h>
#include <stdio.h>
#include "fitimeb.h"
#include <stdlib.h>
#include "time.h"


extern int P,Q,t0,cols[];
//double z[T],z0[T],z1[T],z2[T],z3[T];
double *z, *z0, *z1, *z2, *z3, *zp; //VP
double *zInp = 0;
static double a[DM], *b, *c[S] ;
double (*A)[DM], *B; //VP
static double  ma[DM][DM], mb[DM],ymin[T*M],ymax[T*M],
yav[T*M],yb[T*M],eps[K][T*M],y[K][T*M] ;
static int Tst ;
extern char inputfile[], outputfile[]; //VP

int number_of_variables = Q + S*R  ;

double gauss (), multid_step (int , int , int ),multi_step (int,int,int),progn_err (int , int );
double* get_a (double *b, double *c[S]);


void    set_t (int t),get_z (),out_ab(int,int,int),init_arrays (),eqgauss ();

/*
void get_z0 () 
{
    FILE* fp = fopen (DATAFILE0,"r") ;
    for ( int t = 0 ; t < T/M ; t++ ) 
	fscanf ( fp, "%lg", &z0[t] ) ;


fclose (fp) ;
}

void get_z1 () 
{
    FILE* fp = fopen (DATAFILE1,"r") ;
    for ( int t = 0 ; t < T/M ; t++ ) 
	fscanf ( fp, "%lg", &z1[t] ) ;

 fclose (fp) ;
}

void get_z2 () 
{
    FILE* fp = fopen (DATAFILE2,"r") ;
    for ( int t = 0 ; t < T/M ; t++ ) 
	fscanf ( fp, "%lg", &z2[t] ) ;

fclose (fp) ;
}

void get_z3 () 
{
    FILE* fp = fopen (DATAFILE3,"r") ;
    for ( int t = 0 ; t < T/M ; t++ ) 

  fclose (fp) ;
}

*/
void get_z () 
{
   static bool zInpFilled = false;
   if (!zInpFilled)
   {
      zInp = new double[T];
      zInpFilled = true;
   }
   else goto get_z_copyz;

   {
   FILE* fp = fopen (inputfile, "r") ;
   char str[1255];

   for ( int t = 0, nCols = 0 ; t < T; t += M, nCols = 0) 
   {
      if (fgets(str, 1255, fp) == NULL)
      {
         printf("ERROR-1: Unexpected end of file !\n");
         return;
      }

      char* strCol;
      strCol = (char*)strtok(str, " ");
      for (int i = 0; strCol != NULL && i < MAXCOLS; i++)
      {
         if (cols[i] == 1)
         {
            zInp[t + nCols] = atof(strCol);
//            printf("%d %6.2f; ", t + nCols, zInp[t + nCols]);
 //T E S T I N G   D A T A   I N P U T:
if (INP==1)
{
 printf("  %i    %6.2f\n ", t + nCols, zInp[t + nCols]);
}  
           nCols++;
         }
         strCol = (char*)strtok(NULL, " ");
      }

//      printf("\n");
      if (nCols != M)
{   
 printf("nCols  M  %d  %i \n",nCols,M);    
 printf("ERROR-2: Unexpected end of line!\n");
   }
}  
 fclose (fp) ;
   }

get_z_copyz:
//   if (memcmp(z, zInp, T * sizeof(double)) != 0)
//      printf("MEMCMP z zImp - differs");
   memcpy(z, zInp, T * sizeof(double));
}


void get_AB ()
{
    for (int t = 0 ; t < T ; t++ )
	for ( int i = 0 ; i < P ; i++ )
	    A[t][i] = 0 ;
    for ( int t = 0 ; t < P ; t++ ) 
	B[t] = 0 ;
//int t0=P/M;

    for (int t = t0 ; t < Tst ; t=t+1) 
	for ( int i = 0 ; i < P && M*t-i >0  ; i++ ) {
	    A[M*t][i] = - z[M*t-i-1] ;
for ( int j = 0 ; j < Q && M*(t-j) >0  ; j++ )

		A[M*t][i] += b[j] * A[M*(t-j-1)][i] ;
	    for ( int k = 0 ; k < S && M*( t-k) > 0 ; k++ )
		for ( int l = 0 ; l < R && M*(t-l) > 0 ; l++ )
		    A[M*t][i] -= c[k][l]*z[M*(t-k-1)]*A[M*(t-l-1)][i] ;
	}
    for (int t = t0 ; t < Tst ; t++ ) {
	B[M*t] = z[M*t] ; 
for ( int j = 0 ; j < Q && M*(t-j)> 0 ; j++ ) 
	    B[M*t] += b[j]*B[M*(t-j-1)] ;
        for ( int k = 0 ; k < S && M*(t-k) > 0 ; k++ )
	    for ( int l = 0 ; l < R && M*(t-l) > 0 ; l++ )
		B[M*t] -= c[k][l]*z[M*(t-k-1)]*B[M*(t-l-1)];
    }
/*
printf ("AB: B[M*t] A[M*t] %e %e \n",B[M*t],A[M*t]);
printf ("AB: t t0 Tst M  %i %i %i %i\n",t,t0,Tst,M);
printf ("2 z[0]  z[1]  z[2]  z[3]  z[4] %e %e %e %e %e \n", z[0],z[1],z[2],z[3],z[4]);
*/
/*
printf (" get_AB:
A[2][0] A[2][1]  B[2] b[0] b[1] a[0] a[1] %e %e %e %e %e %e %e \n", A[2][0], A[2][1], B[2],b[0],b[1],a[0],a[1]);
printf ("2 get_AB:: A[4][0] A[4][1] B[4] %e %e %e \n", A[4][0], A[4][1], B[4]);
*/
}
void get_ma_mb()
{
    for ( int i=0 ; i < P ; i++ ) {
	mb[i] = 0 ;
	for ( int j=0 ; j < P ; j++ ) {
            ma[i][j] = 0 ;
for ( int t = t0 ; t < Tst ; t++)
		ma[i][j] += A[M*t][j] * A[M*t][i] ;
    	} 
	mb[i] = 0 ;
for ( int t = t0 ; t < Tst ; t++)
	    mb[i] -= B[M*t] * A[M*t][i] ;
    }
/*
printf ("get_ma: ma[0][0], ma[0][1], ma[1][0], ma[1][1], mb[0], mb[1]b[0] b[1] a[0] a[1] %e %e %e %e
%e %e%e %e %e %e\n",ma[0][0],ma[0][1], ma[1][0],ma[1][1], mb[0],
 mb[1],b[0],b[1],a[0],a[1] );
*/
}

#define AMIN 1e-100

void eqgauss()
{


    int i,j,max;
    double t;
    for(i=0;i<P-1;i++) {
	int k;
	for(j=i+1,max=i;j<P;j++) {
	    if( fabs(ma[j][i]) > fabs(ma[max][i]) )
		max = j ;
	}
	if ( max != i ) {
	    for(k=0;k<P;k++) {
		t = ma[max][k];
		ma[max][k] = ma[i][k];
		ma[i][k] = t;
	    }
	    t = mb[max];
	    mb[max] = mb[i];
	    mb[i] = t;
	}

	for ( j=i+1;j<P;j++) {
	    double s;
	    s=ma[j][i]/( fabs(ma[i][i]) > AMIN ? ma[i][i] : AMIN);
       /*   s=ma[j][i]/ma[i][i];   */
	    for ( k=P-1; k>=i; k--)
		ma[j][k] -= ma[i][k]*s ;
	    mb[j] -= mb[i]*s ;
	}
    }
/*
printf ("01 eqgauss:get ma[0][0],ma[0][1],ma[1][0],ma[1][1], mb[0] mb[1]%e %e %e %e %e %e\n",ma[0][0],ma[0][1],ma[1][0],ma[1][1],mb[0],mb[1] );
*/

 for(i = P-1;i>=0;i--) {
	a[i] = mb[i];
/* 
printf ("for i: eqgauss:get a[0] a[1] b[0] b[1] %e %e %e %e\n", a[0],a[1],b[0], b[1]);
*/
	for(j=P-1; j>i;j--)
	    a[i] -= ma[i][j]*a[j] ;
/* 
printf ("0 for j: eqgauss:get a[0] a[1] b[0] b[1] %e %e %e %e\n", a[0],a[1],b[0], b[1]);
*/
a[i] /= (fabs(ma[i][i])>AMIN ? ma[i][i] : AMIN) ;
/* 
printf ("1 for j: eqgauss:get a[0] a[1] b[0] b[1] %e %e %e %e\n", a[0],a[1],b[0], b[1]);
*/
   }
/*
printf ("1 eqgauss:get ma[0][0],ma[0][1],ma[1][0],ma[1][1], mb[0] mb[1]%e %e %e %e %e %e\n",ma[0][0],ma[0][1],ma[1][0],ma[1][1],mb[0],mb[1] );
 
printf ("eqgauss:get a[0] a[1] b[0] b[1] %e %e %e %e\n", a[0],a[1],b[0], b[1]);
*/
//      a[i] /= ma[i][i];                         
/*
printf (" eqgauss:get ma[0][0],ma[0][1],ma[1][0],ma[1][1], mb[0] mb[1]b[0] b[1] a[0] a[1]%e %e %e %e %e%e %e %e %e %e\n",ma[0][0],ma[0][1],ma[1][0],ma[1][1],mb[0],mb[1],b[0],b[1],a[0],a[1] );   
*/
}


double fm0 (int t0) {
    double sum = 0 ;
for (int t=0;t<T;t++)
yb[t]=0.;
get_AB();
 for ( int t=t0 ; t < Tst ; t++ )
{
	double et = B[M*t] ;
	for ( int j=0 ; j < P ; j++ ) 
	    et += A[M*t][j]*a[j] ;
	sum += et*et ;
/*
   printf ("t  epsilon(t) sumepsilon(t)^2 %i    %e   %e\n",t,et,sum );
*/
/*
	printf ("fm0:sum et B[M*t] A[M*t] %e %e %e %e \n",sum,et,B[M*t],A[M*t]);
printf ("fm0(cont): t t0 Tst M  %i %i %i %i\n",t,t0,Tst,M);
         printf ("sum %e \n", sum );
*/
}
for ( int t=0 ; t < Tst ; t++ )
{
yb[M*t]=B[M*t];
for (int j=0;j<Q && M*(t-j) >0;j++)
{
yb[M*t]=yb[M*t]-B[M*(t-j-1)]*b[j];
}
/*
printf ("fm0:yb[M*t] B[M*t] t  t0  %e   %e   %i   %i \n",yb[M*t],B[M*t],t,t0);
*/
}

sum=sum/(Tst-t0);
return sqrt(sum) ;
}

double frw ( int t0)
{
double sum=0.;
for ( int t = t0 ; t < Tst ; t++ )
{
sum=sum+ (z[M*t] - z[M*(t-1)])*(z[M*t] - z[M*(t-1)]);
}
sum=sum/(Tst-t0);
return sqrt(sum);
	}

double fm () {
    double sum = 0 ;

    for ( int t=t0 ; t < T0/M ; t++ ) {
	double et = B[M*t] ;
	for ( int j=0 ; j < P ; j++ ) 
	    et += A[M*t][j]*a[j] ;
	sum += et*et ;
    }		      
sum=sum/(T0/M-t0);
return sqrt(sum) ;
}

/* Goal function definition */

double fi (const double *x, int n)
{
    b = x ;
/*
printf ("fi 0: n Q P t0 Tst M  %i  %i %i %i %i %i \n",n,Q,P,t0,Tst,M );
printf ("fi 0: x[0]_x[3] %e %e %e %e \n",x[0],x[1],x[2],x[3]);
printf ("fi 0: b[0]- b[3] %e %e %e %e \n",b[0],b[1],b[2],b[3]);
printf ("fi 0: a[0]_a[3] %e %e %e %e \n",a[0],a[1],a[2],a[3]);
*/
for ( int i = 0 ; i < S ; i++)
	c[i] = x + Q + R*i ;
//get_a(b,c);

   get_AB () ;
/*
printf ("after AB : b[0]_b[3] %e %e %e %e \n", b[0],b[1],b[2],b[3]);
printf ("after AB: a[0]a_a[7] %e %e %e %e %e %e %e %e \n", a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7]);
*/
 get_ma_mb () ;
/*
printf ("after ma:  b[0]_b[3] %e %e %e %e  \n",
b[0],b[1],b[2],b[3]);
printf ("after ma: a[0]a_a[7] %e %e %e %e %e %e %e %e \n", a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7]);
*/

   eqgauss () ;
/*
printf ("after eqgauss : b[0]_b[3] %e %e %e %e \n", b[0],b[1],b[2],b[3]);
printf ("aftereq gauss: a[0]a_a[7] %e %e %e %e %e %e %e %e \n", a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7]);  
*/
double sum = 0 ;

    for ( int t=t0 ; t < T0/M ; t++ ) {
	double et = B[M*t] ;
	for ( int j=0 ; j < P ; j++ ) 
	    et += A[M*t][j]*a[j] ;
	sum += et*et ;
    }		      
sum=sum/(T0/M-t0);
double sqsum=sqrt(sum);
/*
printf ("fi 1 : n Q P t0 Tst M  %i %i %i %i  %i %i\n",n,Q,P,t0,Tst,M);
         printf ("fi 1 : sqsum b[0]_b[3] %e %e %e %e %e \n", sqsum,b[0],b[1],b[2],b[3]);
printf ("fi1: a[0]a_a[7] %e %e %e %e %e %e %e %e \n", a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7]);
 */ 
return sqsum;

//    return fm0(t0) ;
}

void constr (const double *x, int n, double *g, int ng)
{
/*
	ng > 0 - calculates all constaints
	ng < 0  - calculates only constraint with number (-ng)
*/
   if ( ng > 0 ) 
	for ( int i = 0 ; i < n ; i++ ) g[i] = 0 ;
   else g[-ng] = 0 ;
}

double* get_a ( double *bb, double *cc[S] )
{
/*
for (int l=0;l<P;l++)
a[l]=0.;
   b = bb ;
for (int k=0;k<P;k++)
a[k]=0.;
*/
get_z();
/*
printf("get_a: z[0] z[1] z[2] z[3] %e %e %e  %e %e %e %e  %e\n",z[0],
z[1],z[2],z[3],z[4],z[5],z[6],z[7]);
*/
/*
printf ("get_a: b[0]_b[3] a[0]_a[3] %e %e %e %e %e %e %e %e \n", b[0],b[1],b[3],b[4],a[0],a[1],a[2],a[3]);
*/
for ( int i = 0 ; i < S ; i++ )
	c[i] = cc[i] ;
    get_AB () ;
/*
printf ("get_a after_AB: a[0]-a[3] b[0]_b[3] %e %e %e %e  %e %e %e %e\n", a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3]);
*/
    get_ma_mb () ;
/*
printf ("get_a after_ma: a[0]_a[1] b[0]_b[3] %e %e %e %e %e %e %e %e\n",a[0],a[1],a[2],a[3], b[0],
 b[1],b[2],b[3]);
*/
   eqgauss () ;
/*
printf ("get_a after_eqgauss: a[0]_a[3] b[0]_b[3]  %e %e %e %e %e %e %e %e\n", a[0],a[1],a[2],a[3],b[0],b[1],b[2],b[3]);
*/
    return a ;
}

void set_t (int t)
{
    Tst = t ;
}
/*
void out_ab(int Po,int Qo,int to))
{
P=Po;
Q=Qo;
t0=to;
fi0(P,Q,t0)
// FILE *fp = fopen ( "params.out", "wr") ;
printf ( "a %e %e %e %e %e %e %e %e %e %e\n", a[0], a[1], a[2]), a[3], a[4], a[5], a[6], a[7], a[8], a[9];
	//for ( m = 0 ; m < NPARMS ; m++ )
	    printf ( "b %e %e %e %e %e %e %e %e\n", b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7]);
	printf ( "v %e\n", fval) ;
//fclose (fp);
}
*/
double multi_step (int t1, int t2, int t3)
{

double large=10.e+20;
double sigma=0.;
double average=0.;

 FILE *fp = fopen ( "test.out", "w") ;
for (int t=0;t<T*M;t++)
{
ymin[t]=large;
ymax[t]=-large;
yav[t]=0.;
for (int k=0;k<K;k++)
{
eps[k][t]=0.;
y[k][t]=0.;
}
}

for (int k=0;k<K;k++)
for (int t=t1;t<t2;t++)
{
for (int i=0;i<P && M*t-i > 0;i++)
eps[k][M*t]=eps[k][M*t]-a[i]*z[M*t-i-1];
for (int j=0;j<Q && M*(t-j) > 0;j++)
{
eps[k][M*t]=eps[k][M*t]+b[j]*eps[k][M*(t-j-1)];
}
eps[k][M*t]=eps[k][M*t]+z[M*t];

fprintf (fp, "equal0 eps[0][M*t] eps[1][M*t]  t  k  %e   %e    %i  %i   %i   
\n",eps[0][M*t],eps[1][M*t],t,k,i);
}

for (int t=t1;t<t2;t++) 
{
sigma=sigma+eps[0][M*t]*eps[0][M*t];


fprintf (fp, "sigma eps[0][M*t] eps[0][0] b[0] a[0] a[1] t t1 t2 %e  %e  %e  %e  %e  %e  %i %i  %i  \n",sigma,eps[0][M*t],eps[0][0],b[0],a[0],a[1],t,t1,t2);
fprintf (fp, "z[0] z[M*t] t  %e  %e  %i  \n",z[0],z[M*t],t);

average=average+eps[0][M*t];


fprintf (fp, "average[M*t] t %e  %i \n",average,t) ;

}
sigma=sigma/(t2-t1);
average=average/(t2-t1);
average=average*average;

fprintf (fp, "average**2   %e  \n",average) ;

sigma=sigma-average;
sigma=sqrt (sigma);
if (F==0)
{
sigma=0.;
}
printf ("Finish: sigma %e  \n",sigma) ;
fprintf (fp, "sqrt  %e  \n",sigma) ;

for (int k=0;k<K;k++)
for (int t=t1;t<t2*M;t++)
{
y[k][t]=z[t];
/*
fprintf (fp, "det eps[k][M*t]  y[k][M*t]  z[M*t]  t  k  t2 t3 %e   %e   %e    %i  %i   %i  %i
\n",eps[k][M*t],y[k][M*t],z[M*t],t,k,t2,t3);
*/
}

double noise;
for (int k=0;k<K;k++)
for (int t=t2;t<t3;t++)
{
long ltime=time(NULL);
int stime=(unsigned int) ltime/2;
srand(stime);
noise=gauss();
//noise = noise + noise;
//noise =sigma* noise/(t3-t2);

eps[k][M*t]=noise*sigma;

fprintf (fp, "rnd eps[k][M*t] nois  t k  t2  t3  %e  %e   %i  %i  %i  %i
\n",eps[k][M*t],noise,t,k,t2,t3);
}
fprintf (fp, "rnd noise   %e \n",noise);
for (int k=0;k<K;k++)
for (int t=t2;t<t3;t++)
for (int i=1;i<M;i++)
{
y[k][M*t-i]=z[M*t2-i];


fprintf (fp, "detygr y[k][M*t-i]  z[M*t2-i]  t  k  i %e   %e    %i  %i   %i   
\n",y[k][M*t-i],z[M*t2-i],t,k,i);
}

for (int k=0;k<K;k++)
{
for (int t=t2;t<t3;t++)
{
for (int i=0;i<P && M*t-i > 0;i++)
{
y[k][M*t]=y[k][M*t]+a[i]*y[k][M*t-i-1];

fprintf (fp, "y[k][M*t] y[k][M*t-i-1] a[i] i t k t2 t3  %e    %e    %e      %i   %i   %i   %i  %i 
\n",y[k][M*t],y[k][M*t-i-1],a[i],i,t,k,t2,t3);
}
for (int j=0;j<Q && M*(t-j) > 0 ;j++)
{
y[k][M*t]= y[k][M*t]- b[j]*eps[k][M*(t-j-1)];

fprintf (fp, "y[k][M*t] eps[k][M*(t-j-1)] b[j] t k t2 t3  %e    %e   %e     %i   %i   %i   %i  %i 
\n",y[k][M*t],eps[k][M*(t-j-1)],b[j],j,t,k,t2,t3);
}
for (int i=0;i<P && M*t-i >0;i++)
{

fprintf (fp, "first y[k][M*t-i]  eps[k][M*t] t k i  %e  %e   %i  %i  %i
\n",y[k][M*t-i],eps[k][M*t],t,k,i);
}
y[k][M*t]=y[k][M*(t-1)]+eps[k][M*t];
}
}
for (int k=0;k<K;k++)
{
fprintf (fp, "second y[k][M*t2] eps[k][M*t2] t2 k  %e  %e  %i  %i
\n",y[k][M*t2],eps[k][M*t2],t2,k);
}
for (int t=t2;t<t3;t++)
for (int k=0;k<K;k++)
{
yav[M*t]=yav[M*t]+y[k][M*t];

fprintf (fp, "yav[M*t] y[k][M*t] y[0][M*t]  y[1][M*t]t k %e %e %e %e %i %i
\n",yav[M*t],y[k][M*t],y[0][M*t],y[1][M*t],t,k);
if (y[k][M*t] < ymin[M*t])
ymin[M*t] = y[k][M*t];

fprintf (fp, "ymin[M*t]y[k][M*t] t k %e %e %i %i \n",ymin[M*t],y[k][M*t],t,k);
if (y[k][M*t] > ymax[M*t])
ymax[M*t] = y[k][M*t];

fprintf (fp, "ymax[M*t] y[k][M*t] t k %e %e  %i  %i
\n",ymax[M*t],y[k][M*t],t,k);
}
for (int t=t2;t<t3;t++)
{
yav[M*t]=yav[M*t]/K;
fprintf (fp, "gale: yav[M*t] t k %e  %i 
\n",yav[M*t],t);
}
//fclose (fp);
  for (int t=t2;t<t3;t++)
   {
fprintf (fp, "gale: t  ymin[M*t]  yav[M*t]  ymax[M*t]      %i    %e    %e    %e
\n",t,ymin[M*t],yav[M*t],ymax[M*t]);
   }
fprintf (fp, "sigma and average %e   %e  \n",sigma,average);
double delt=0.;

for (int t=t2;t<t3;t++)
{
delt=delt+(yav[M*t]-z[M*t])*(yav[M*t]-z[M*t]);
fprintf (fp, "multi-step for: t t2  t3  ymin[M*t]  yav[M*t]  ymax[M*t]  z[M*t]    %i   %i   %i   %e   %e    %e    %e
\n",t,t2,t3,ymin[M*t],yav[M*t],ymax[M*t],z[M*t]);
}
delt=delt/(t3-t2);
delt=sqrt (delt);
fprintf (fp, "delt t0 T1 T:
%e      %i  %i  %i   \n",delt,t1,t2,t3);
fclose (fp) ;
return delt;
}


double progn_err (int t2, int t3)
{
double delt=0.;
 FILE *fp = fopen ( "progn.out", "w") ;

for (int t=t2;t<t3;t++)
{
delt=delt+(yav[M*t]-z[M*t])*(yav[M*t]-z[M*t]);
}
delt=delt/(t3-t2);
delt=sqrt (delt);
if (TR <0)
{
 for (int t=t2;t<t3;t++)
   {
/*
fprintf (fp, "progn t  ymin[M*t]  yav[M*t]  ymax[M*t]      %i    %e    %e    %e
\n",t,ymin[M*t],yav[M*t],ymax[M*t]);/

fprintf (fp, "    %i    %e  
\n",t,ymin[M*t]);


fprintf (fp, "    %i    %e  
\n",t,yav[M*t]);
*/

fprintf (fp, "    %i    %e  
\n",t,ymax[M*t]);

  }
//fprintf (fp, "progn_err %e   \n",delt);
return delt;
}
//fclose (fp);
//FILE *fp = fopen ( "regr.out", "w") ;
else
{
double reg[T/M];
int tpr=TR/M;
int tend=TE/M;
for (int j=tpr;j<tend;j++)
{
reg[j]=0.;
}

for (int j=tpr;j<tend;j++)
for (int i=0;i<P-1;i++)
{
reg[j]=reg[j]+a[i]*z[M*j-i-1+M];
/*
fprintf (fp, "loop values reg z[TR-i-1+M] a[i] i =  %e  %e %e %i \n",reg,z[M*tpr-i-1+M],a[i],i);
*/
}
double reg_f[T/M];
for (int j=tpr;j<tend;j++)
{
reg_f[j]=0.;
}
for (int j=tpr;j<tend;j++)
for (int i=0;i<P-1;i++)
{
reg_f[j]=reg_f[j]+a[i]*z[M*j-i-1+M];
/*
fprintf (fp, "loop values reg_f z[TR-i-1+M] a[i] i =  %e  %e %e %i \n",reg_f,z[M*tpr-i-1+M],a[i],i);
*/
}
double reg_err[T/M];
double reg_averr=0.;
for (int j=tpr;j<tend;j++)
{
reg_err[j]=0.;
}
for (int j=tpr;j<tend;j++)
{
reg_err[j]=reg[j]-z[M*j-P+1-1+M];
fprintf (fp, "simple regression and true value =  %e    %e   %i   \n",reg[j],z[M*j-P+1-1+M],j);
fprintf (fp, "regression error  =  %e    %i \n",reg_err[j],j);
reg_averr=reg_averr+fabs(reg_err[j]);
/*
fprintf (fp, "complete regression reg_f[j]  j =  %e   %i\n",reg_f[j],j);
*/
}
reg_averr=reg_averr/(tend-tpr);
//reg_averr=sqrt(reg_averr);
fprintf (fp, "in the line  TR/M  TE/M  reg_averr=   %i   %i
%e   \n",tpr,tend,reg_averr);


for (int i=0;i<P-1;i++)
{
fprintf (fp, "regression parameters a[i] i =  %e   %i \n",a[i],i);
}
return delt;
}
fclose (fp);
}



double multid_step (int t1, int t2, int t3)
{

double large=10.e+20;
double sigma=0.;
double average=0.;

 FILE *fp = fopen ( "ms.out", "w") ;
for (int t=0;t<T;t++)
{
ymin[t]=large;
ymax[t]=-large;
yav[t]=0.;
for (int k=0;k<K;k++)
{
eps[k][t]=0.;
y[k][t]=0.;
}
}

for (int k=0;k<K;k++)
for (int t=t1;t<t2;t++)
{
for (int i=0;i<P && M*t-i > 0;i++)
eps[k][M*t]=eps[k][M*t]-a[i]*z[M*t-i-1];
for (int j=0;j<Q && M*(t-j) > 0;j++)
{
eps[k][M*t]=eps[0][M*t]+b[j]*eps[k][M*(t-j-1)];
}
eps[k][M*t]=eps[k][M*t]+z[M*t];
}

for (int t=t1;t<t2;t++) 
{
sigma=sigma+eps[0][M*t]*eps[0][M*t];


fprintf (fp, "sigma eps[0][M*t] eps[0][0] b[0] a[0] a[1] t t1 t2 %e  %e  %e  %e  %e  %e  %i %i  %i  \n",sigma,eps[0][M*t],eps[0][0],b[0],a[0],a[1],t,t1,t2);
fprintf (fp, "z[0] z[M*t] t  %e  %e  %i  \n",z[0],z[M*t],t);

average=average+eps[0][M*t];


fprintf (fp, "average[M*t] t %e  %i \n",average) ;

}
sigma=sigma/(t2-t1);
average=average/(t2-t1);
average=average*average;

fprintf (fp, "average**2   %e  \n",average) ;

sigma=sigma-average;
sigma=sqrt (sigma);

if (F==0)
{
sigma=0.;
}

fprintf (fp, "sqrt  %e  \n",sigma) ;

for (int k=0;k<K;k++)
for (int t=t1;t<t2*M;t++)
{
y[k][t]=z[t];

fprintf (fp, "det eps[k][M*t]  y[k][M*t]  z[M*t]  t  k  t2 t3 %e   %e   %e    %i  %i   %i  %i
\n",eps[k][M*t],y[k][M*t],z[M*t],t,k,t2,t3);
}

for (int k=0;k<K;k++)
for (int t=t2;t<t3;t++)
for (int i=1;i<M;i++)
{
y[k][M*t-i]=z[M*t2-i];
fprintf (fp, "detygr y[k][M*t-i]  z[M*t2-i]  t  k  i %e   %e    %i  %i   %i   
\n",y[k][M*t-i],z[M*t],t,k,i);
}


double noise;
for (int k=0;k<K;k++)
for (int t=t2;t<t3;t++)
{
long ltime=time(NULL);
int stime=(unsigned int) ltime/2;
srand(stime);
noise=gauss();
eps[k][M*t]=sigma*noise;
noise = noise + noise;
fprintf (fp, "rnd eps[k][M*t] nois  t k  t2  t3  %e  %e   %i  %i  %i  %i
\n",eps[k][M*t],noise,t,k,t2,t3);

}
noise = noise/(t3-t2);
fprintf (fp, "rnd noise   %e   
\n",noise);
for (int k=0;k<K;k++)
for (int t=t2;t<t3;t++)
{
for (int i=0;i<P && M*t-i > 0;i++)
{
y[k][M*t]=y[k][M*t]+a[i]*y[k][M*t-i-1];

fprintf (fp, "y[k][M*t]  a[i] t k t2 t3  %e    %e    %i   %i   %i   %i  %i 
\n",y[k][M*t],a[i],i,t,k,t2,t3);

}
for (int j=0;j<Q && M*(t-j) > 0 ;j++)
{
y[k][M*t]= y[k][M*t]- b[j]*eps[k][M*(t-j-1)];
fprintf (fp, "y[k][M*t] eps[k][M*(t-j-1)] b[j] t k t2 t3  %e    %e   %e     %i   %i   %i   %i  %i 
\n",y[k][M*t],eps[k][M*(t-j-1)],b[j],j,t,k,t2,t3);

}
y[k][M*t]=y[k][M*t]+eps[k][M*t];
fprintf (fp, "first y[k][M*t] y[k][M*(t-1)] eps[k][M*t] t k  %e  %e  %e    %i  %i
\n",y[k][M*t],y[k][M*(t-1)],eps[k][M*t],t,k);
}
for (int k=0;k<K;k++)
{
fprintf (fp, "second y[k][M*t2] eps[k][M*t2] t2 k  %e  %e  %i  %i
\n",y[k][M*t2],eps[k][M*t2],t2,k);
}

for (int t=t2;t<t3;t++)
for (int k=0;k<K;k++)
{
yav[M*t]=yav[M*t]+y[k][M*t];

fprintf (fp, "yav[M*t] y[k][M*t] t k %e %e  %i %i
\n",yav[M*t],y[k][M*t],t,k);

if (y[k][M*t] < ymin[M*t])
ymin[M*t] = y[k][M*t];

fprintf (fp, "ymin[M*t]y[k][M*t] t k %e %e %i %i \n",ymin[M*t],y[k][M*t],t,k);

if (y[k][M*t] > ymax[M*t])
ymax[M*t] = y[k][M*t];

fprintf (fp, "ymax[M*t] y[k][M*t] t k %e %e  %i  %i
\n",ymax[M*t],y[k][M*t],t,k);

}

for (int t=t2;t<t3;t++)
{
yav[M*t]=yav[M*t]/K;
fprintf (fp, "yav[M*t] t k %e  %i 
\n",yav[M*t],t);
}
   for (int t=t2;t<t3;t++)
   {
fprintf (fp, "t  ymin[M*t]  yav[M*t]  ymax[M*t]      %i    %e    %e    %e
\n",t,ymin[M*t],yav[M*t],ymax[M*t]);
   }
fprintf (fp, "sigma and average %e   %e  \n",sigma,average);
double delt=0.;

for (int t=t2;t<t3;t++)
{
delt=delt+(yav[M*t]-z[M*t])*(yav[M*t]-z[M*t]);
}
delt=delt/(t3-t2);
delt=sqrt (delt);
fprintf (fp, "delt t0 T1 T:
%e      %i  %i  %i   \n",delt,t1,t2,t3);
fclose (fp) ;
return delt;
}


double gauss()
{
	static int iset=0;
	static float gset;
	float fac,r,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*RAND()-1.0;
			v2=2.0*RAND()-1.0;
			r=v1*v1+v2*v2;
		} while (r >= 1.0);
		fac=sqrt(-2.0*log(r)/r);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}