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


#define BAYES1 1
#define GLOPT 2
#define EXKOR 3
#define NLP 4

#define NGLOBIT MAX(GLOPT_MAX_IT,BAYES1_MAX_IT)
#define NLOCIT MAX(NLP_MAX_IT,EXKOR_MAX_IT)
#define NIT MAX(NGLOBIT,NLOCIT)
#define NPARMS (Q+S*R)

int Po,Qo,epsprn,ig,il; //VP
int P,Q,t0, cols[MAXCOLS]; //VP
double fm[NIT], *xn[NIT] ;
double parms[DM], min[DM], max[DM],rnd[SI] ;
double *a,*b, *ao,*bo,*c[S],delta,delta1,opt,opt1;

extern double *z, *z0, *z1, *z2, *z3, *zp; //VP
extern double *zInp; //VP
extern double (*A)[DM], *B; //VP

char inputfile[100], outputfile[100];//VP

int initialize(char* ininame);
void shutdown();

void set_bounds(), get_glob_parms ( double *, double *),
     get_local_parms ( double *, double *),
     init_arrays(),get_z(),set_t(int),out_ab(int,int,int)   ;
double fi0(int , int , int ), fm0(int), frw(int), gauss(), multid_step (int , int , int ), progn_err (int , int ),multi_step (int,int,int);
double* get_a (double *b, double *c[S]); 
main ()
{
int P1,Q1,t01,Po,Qo,to,ts,tmax;

initialize(INIFILE);
printf("Number of Factors M %i \n",M);
printf("Number of of columns of matrix c R %i \n",R);

/*
printf("End of Learning Set T0 %i \n",T0);
printf("End of Testing Set %i \n",T1);
printf("End of Complete Set %i \n",T);
*/
printf("INIFILE %s \n",INIFILE);
printf("Multi-step prognose indicator V %i \n",V);
printf("Multi-step optimization indicator W %i \n",W);
printf("Number of prognose repetitions K %i \n",K);
printf("indicator of simulated annealing SA %i \n",SA);
printf ("temperature of simulated annealing ST %e \n",ST);
printf ("number of simulated annealing iterations  SI %i  \n",SI );
printf ("indicator of residual printing EPS %i  \n",EPS);
printf ("indicator of input control INP %i  \n",INP);
printf ("the first line in datafile for regression, TR  %i  \n",TR);
printf ("the  last line in datafile for regression, TE  %i  \n",TE);
printf ("starting number of AR parameters Ps %i  \n",Ps );
printf ("starting number of MA parameters Qs %i  \n",Qs );
printf ("minimal number of AR parameters Pmin %i  \n",Pmin );
printf ("minimal number of MA parameters Qmin %i  \n",Qmin );
printf ("maximal number of AR parameters Pmax %i  \n",Pmax );
printf ("maximal number of MA parameters Qmax %i  \n",Qmax );
printf ("number of data entries T %i  \n",T );
printf ("number of  entries for a and b optimization T0 %i  \n",T0 );
printf ("number of  entries for P and Q optimization T1 %i  \n",T1 );
printf ("indicator of variance F %i  \n",F );

/*
if (M==1)
printf("Datafile0 %s \n",DATAFILE0);
if (M==2)
{
printf("Datafile0 %s \n",DATAFILE0);
printf("Datafile1 %s \n",DATAFILE1);
}
if (M==3)
{
printf("Datafile0 %s \n",DATAFILE0);
printf("Datafile2 %s \n",DATAFILE1);
printf("Datafile3 %s \n",DATAFILE2);
}
if (M==4)
{
printf("Datafile0 %s \n",DATAFILE0);
printf("Datafile1 %s \n",DATAFILE1);
printf("Datafile2 %s \n",DATAFILE2);
printf("Datafile3 %s \n",DATAFILE3);
}
*/
 ig=0;
 il=0;
int dP=M;
int dQ=1;
int dt=1;
ts=0;
tmax=ts;
printf("Global Max Iterat %i \n",BAYES1_MAX_IT);
printf("Local Max Iterat %i \n",EXKOR_MAX_IT);
printf("b-Bounds %e \n",MAX_B_BOUND);
printf("Diskret Params Init (Ps-AR params, Qs-MA params,ts-starting time) :  Ps Qs ts  = %i %i %i  \n",Ps,Qs,ts);
printf("Diskret  Params Step :  dP dQ dt = %i %i %i \n",dP,dQ,dt);
printf("Diskret  Params Max :  Pmax Qmax tmax = %i %i %i \n",Pmax,Qmax,tmax);




P=Ps;
Q=Qs;
t0=ts;
Po=Ps;
Qo=Qs;
to=ts;
 
opt=fi0(P,Q,t0);
bo=b;
ao=a;
for (int l=0;l<Ps;l++)
printf("Start AR:  l=0_Ps a[l] =  %e \n",a[l]);
for (int k=0;k<Qs;k++)
{
printf("Start MA: k=0_Qs b[k]  =  %e  \n",b[k]);
/*
printf("Start parms: k=0_Qs parms[k]  = %e \n",parms[k]);
*/
}
printf("Start Results: opt delta_perc delta_mean= %e  %e  %e  \n",opt,delta,delta1);
printf("Start Discret: Po Qo to   opt   delta_perc  delta_mean = %i %i %i %e %e %e  \n",Po,Qo,to,opt,delta,delta1);
for (int k=0;k<Qo;k++)
{
printf("Start MA: k=0_Qo b[k] = %e  \n",b[k]);
}
for (int l=0;l<Po;l++)
{
printf("Start  AR:  l=0_Po  a[l] =  %e  \n",a[l]);
}
if (Ps>Pmax)
    printf("Error:Ps>Pmax\n");
if (Ps<Pmin)
    printf("Error:Ps<Pmin\n");
if (Qs>Qmax)
    printf("Error:Qs>Qmax\n");
if (Qs<Qmin)
    printf("Error:Qs<Qmin\n");
if (P>DM/M)
    printf("Warning:P>DM/M\n");
if (T1 != 2*T0)
printf("Warning:T1 != 2*T0\n");
if (P>T0/M)
    printf("Warning:P>T0/M\n");

if (SA==1)
{
double eks;
P=Ps;
Q=Qs;
//for (int j=0;j<SI;j++)
for (int i=0;i<SI;i++)
{

 rnd[i]=RAND ();
if (rnd[i]>= 0.5 && P+M-1 < Pmax)
P=P+M;
if (rnd[i] < 0.5 && P-M+1 > Pmin)
P=P-M;

 rnd[i]=RAND ();
if (rnd[i] >= 0.5 && Q < Qmax)
Q=Q+1;
 rnd[i]=RAND ();
if (rnd[i] < 0.5 && Q > Qmin)
Q=Q-1;
/*
printf (" before opt1= fi0 i P Q t0 rnd[i] opt1 opt   %i   %i   %i   %i     %e    %e    %e \n",i,P,Q,t0,rnd[i],opt1, opt);
for (int i=0;i<P;i++)
printf (" before opt1= fi0 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" before opt1= fi0 i  b[i] %i   %e \n",i,b[i]);
*/

opt1=fi0(P,Q,t0);

/*
for (int i=0;i<P;i++)
printf (" after opt1= fi0 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt1= fi0 i  b[i] %i   %e \n",i,b[i]);

printf (" after opt1=f0 i P Q t0 rnd[i] opt1 opt   %i   %i   %i   %i     %e    %e    %e \n",i,P,Q,t0,rnd[i],opt1, opt); 
*/

if (opt1 < opt)
{opt=opt1;
 Po=P;
 Qo=Q;
to=t0;
 bo=b; ao=a;
/*
printf (" after opt= opt1 i P Q t0 rnd[i] opt1 opt   %i   %i   %i   %i     %e    %e    %e \n",i,P,Q,t0,rnd[i],opt1, opt);
for (int i=0;i<P;i++)
printf (" after opt1= opt1 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt1= opt1 i  b[i] %i   %e \n",i,b[i]);
*/

}

/*
printf (" before opt1>= opt i P Q t0 rnd[i] opt1 opt   %i   %i   %i   %i     %e    %e    %e \n",i,P,Q,t0,rnd[i],opt1, opt); 
*/


if (opt1 >= opt)
{
rnd[i]=RAND();
eks=ldexp(1.,(opt-opt1)/(ST*0.6931));
if (rnd[i] > eks)
{
 Po=P;
 Qo=Q;
to=t0;


/*
printf (" after opt1 >= opt i P Q t0 rnd[i] opt1 opt   %i   %i   %i   %i     %e    %e    %e \n",i,P,Q,t0,rnd[i],opt1, opt); 
for (int i=0;i<P;i++)
printf (" after opt1>= opt i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt1>= opt i  b[i] %i   %e \n",i,b[i]);
*/

}
}


}
}


if (SA==0)
{
for ( P=Pmin; P<Pmax+1; P=P+dP)
for ( Q=Qmin; Q<Qmax+1; Q=Q+dQ)
for ( t0=ts; t0<tmax+1; t0=t0+dt)
{

/*
printf (" before opt1= fi0  P Po Q Qo t0 opt1 opt      %i   %i   %i   %i   %i 
    %e    %e \n",P,Po,Q,Qo,t0,opt1, opt);
for (int i=0;i<P;i++)
printf (" before opt1= fi0 i a[i]  %i    %e  \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" before opt1= fi0 i  b[i] %i   %e  \n",i,b[i]);
*/

opt1=fi0(P,Q,t0);

/*
for (int i=0;i<P;i++)
printf (" after opt1= fi0 i a[i]  %i      %e\n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt1= fi0 i  b[i]  %i    %e \n",i,b[i]);

printf (" after opt1=f0 i P Po Q Qo t0 opt1 opt   %i   %i   %i   %i    %i   %i    %e    %e \n",i,P,Po,Q,Qo,t0,opt1, opt); 
*/

if (opt1 < opt)
{opt=opt1;

/*
printf (" after opt= opt1 i P Q t0 opt1 opt   %i   %i   %i   %i      %e    %e \n",i,P,Q,t0,opt1, opt);
for (int i=0;i<P;i++)
printf (" after opt1= opt1 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt1= opt1 i  b[i] %i   %e \n",i,b[i]);
*/


 Po=P;
 Qo=Q;
to=t0;


/*
printf("before Po=P ... P Po Q Qo
%i  %i  %i   %i  \n",P,Po,Q,Qo);
for (int l=0;l<P;l++)
{
printf("before Po=P ...  a[l] = %e  \n",a[l]);
}
*/

}
}
}

P=Po;
Q=Qo;
t0=to;

/*
printf("before F==1 P Po Q  Qo
%i  %i  %i   %i  \n",P,Po,Q,Qo);

for (int l=0;l<P;l++)
{
printf("before F==1 l=0_P  a[l] = %e  \n",a[l]);
}

printf (" before opt= fi0  P Q t0 opt opt      %i   %i   %i   
    %e    %e \n",P,Q,t0,opt, opt);
for (int i=0;i<P;i++)
printf (" before opt1= fi0 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" before opt1= fi0 i  b[i] %i   %e \n",i,b[i]);
*/
epsprn=1;
opt=fi0 (P,Q,t0);

//printf (" epsprn %i \n",epsprn);
/*
for (int i=0;i<P;i++)
printf (" after opt= fi0 i a[i]  %i   %e \n",i,a[i]);
for (int i=0;i<Q;i++)
printf (" after opt= fi0 i  b[i] %i   %e \n",i,b[i]);

printf (" after opt=f0  P Q t0  opt1 opt     %i   %i   %i  
   %e    %e     \n",P,Q,t0,opt1, opt); 
*/

printf("Finish Discret: Po Qo toP Q t0 opt delta_perc delta_mean= %i  %i  %i       %i  %i  %i      %e  %e  %e   \n",Po,Qo,to,P,Q,t0,opt,delta,delta1);  
for (int k=0;k<Q;k++)
{
printf("Finish MAo: k=0_Q b[k] =  %e \n",b[k]);
} 
for (int l=0;l<P;l++)
{
printf("Finish  ARo:  l=0_P  a[l] = %e  \n",a[l]);
}
/*
double reg=0.;
for (int i=1;i<P;i++)
{
reg=reg+a[i]*z[T/M-i];
}
printf("simple_regression v=0 %e \n",reg);
*/
if (PL==1)
{
double D1=(B1-A1)/DN;
double xplot[100];
 FILE *fp = fopen ( "plot.out", "w") ;
 for (  xplot[PLOT1] = A1 ; xplot[PLOT1] <=B1  ; xplot[PLOT1]+= D1 )
 {
for (int i=0;i<Q;i++)   
{
if (i==PLOT1)
b[i]=xplot[i];
}
 double fiplot=log (fi(b,Q));
/*
printf ( "b[0] b[1]  fiplot  %e  %e  %e\n", b[0],b[1],fiplot ) ;
*/
fprintf (fp, " %e  \n",b[PLOT1]);
 fprintf (fp, " %e  \n",fiplot);
fprintf (fp, "%s   \n"," ");
 }
fclose (fp);
}

if (PL==2)
{
double D1=(B1-A1)/DN;
double D2=(B2-A2)/DN;
double xplot[100];
 FILE *fp = fopen ( "plot.out", "w") ;
 for (  xplot[PLOT1] = A1 ; xplot[PLOT1] <=B1  ; xplot[PLOT1]+= D1 )
 {
for (  xplot[PLOT2] = A2 ; xplot[PLOT2] <=B2  ; xplot[PLOT2]+= D2 )
{
for (int i=0;i<Q;i++)   
{
if (i==PLOT1)
b[i]=xplot[i];

if (i==PLOT2)
b[i]=xplot[i];


}
double fiplot= log (fi(b,Q));
/*
fprintf (fp, " %e  \n",b[PLOT1]);
fprintf (fp, " %e  \n",b[PLOT2]);
 fprintf (fp, " %e  \n",fiplot);
*/

 fprintf (fp, " %e %e %e \n",b[PLOT1],b[PLOT2],fiplot);
//fprintf (fp, "%s   \n"," ");
}
//fprintf (fp, "%s   \n"," ");
}
fclose (fp);
}
if (V==1)
{
multi_step (t0,T1/M,T/M);
double progn;
progn= progn_err (T1/M, T/M);
printf("Finish progn_err %e \n",progn);
/*
double reg=0.;
for (int i=1;i<P;i++)
{
reg=reg+a[i]*z[T/M-i];
}
printf("simple_regression v=1 %e \n",reg);
*/
}
//multi_step (t0,T1/M,T/M);
}

double fi0(int p, int q , int s )
{
//printf("fi0-1 %i % i % i\n",p,q,s);
    double fval ;
//printf("fi0-1.1 %i % i % i\n",p,q,s);
    get_z () ;
//printf("fi0-2 %e % e % e\n",p,q,s);
    set_bounds () ;
    set_t (T0/M) ;
   init_arrays() ;
//printf("fi0-3 %e % e % e\n",p,q,s);
if (NPARMS > 0)
{
get_glob_parms ( parms, &fval ) ;
get_local_parms ( parms, &fval ) ;
//   init_arrays() ;
//printf("fi0-4 %e % e % e\n",p,q,s);
a=get_a(b,c);
//printf("fi0-5 %e % e % e\n",p,q,s);
}
else
     {  
//  init_arrays() ;
//printf("fi0-6 %e % e % e\n",p,q,s);
a=get_a(b,c);
//printf("fi0-7 %e % e % e\n",p,q,s);
}

double delt_p0;
double rw_p0;
set_t(T0/M);
    delt_p0=fm0(P/M);
    rw_p0=frw(P/M);

double delt_01;
    double rw_01;
    set_t(T1/M);
if (W==0)
{
delt_01=fm0(T0/M);
}

if (W==1)
{
delt_01=multid_step(t0,T0/M,T1/M);

}
  rw_01=frw(T0/M);

double delt_1t;
    double rw_1t;
 set_t(T/M);
 delt_1t=fm0(T1/M);
rw_1t=frw(T1/M);
/*
printf ("fi0_1t: Tst delt_1t rw_1t %i %e %e\n",T/M,delt_1t,rw_1t);
*/
delta=100*(delt_1t-rw_1t)/rw_1t;
delta1=delt_1t;
/*
printf ("fi0_delta: Tst delta %i %e \n",T/M,delta);
*/
    double delt_0t;
    double rw_0t;
set_t(T/M);
delt_0t=fm0(T0/M);
rw_0t=frw(T0/M);
/*
printf ("fi0_0t: Tst delt_0t rw_0t %i %e %e\n",T0/M,delt_0t,rw_0t);
*/
   double delt_pt;
    double rw_pt;
set_t(T/M);
 delt_pt=fm0(P/M);
rw_pt=frw(P/M);
/*
printf ("fi0_pt: Tst delt_pt rw_pt %i %e %e\n",T/M,delt_pt,rw_pt);
printf("fi0_end: NPARMS P Q t0 delta_perc  %i  %i %i %i  %e \n",NPARMS,P,Q,t0,delta);
printf("fi0_end:a[0]_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]);
printf("fi0_end: b[0]_b[3] %e %e %e %e \n",b[0],b[1],b[2],b[3]);
*/
return delt_01;
/*
P=p;
Q=q;
t0=s;
*/
}




void set_bounds ()

{
    for (int i = 0 ; i < Q ; i++ ) {
	min[i] = - MAX_B_BOUND ;
	max[i] = MAX_B_BOUND ;
    }
    for ( ; i < NPARMS ; i++ ) {
	min[i] = - MAX_C_BOUND ;
	max[i] = MAX_C_BOUND ;
    }
}


void get_glob_parms ( double *parms, double *fval)
{
#if (GLOBAL_METH==BAYES1)
    bayes1 ( parms, min, max, NPARMS, fval, fi, BAYES1_MAX_IT, BAYES1_LT, fm, xn) ;
if (ig==0)
{
printf ("BAYES,\n");
ig++;
}
#else

#if (GLOBAL_METH==GLOPT)
    glopt ( parms, min, max, NPARMS, fval, fi, GLOPT_MAX_IT, GLOPT_LT, GLOPT_MAXL, fm, xn) ;
if (ig==0)
{
printf ("GLOPT,\n");
ig++;
}
#endif
#endif
} 


void get_local_parms ( double *parms, double *fval)
{
#if (LOCAL_METH==NLP)
    fnlpql ( parms, min, max, NPARMS, fval, fi, constr, NLP_MAX_IT, NLP_M, NLP_ME, fm, xn) ;

if (il==0)
{
printf ("NLP,\n");
il++;
}
#else
#if (LOCAL_METH==EXKOR)
    exkor ( parms, min, max, NPARMS, fval, fi, EXKOR_MAX_IT, EXKOR_INIT_POINTS, fm, xn) ;

if (il==0)
{
printf ("EXKOR,\n");
il++;
}
#else
    fprintf ( stderr, "LOCAL_METH in fitimeb.h should be NLP or EXKOR.\n") ;
    exit(1) ;
#endif
#endif
}


void init_arrays ()
{
    for ( int i = 0 ; i < NIT ; i++ )
	if (( xn[i] = (double*) malloc ( sizeof(double)*NPARMS)) == NULL ) {
	    fprintf (stderr, "Could'nt allocate memory\n") ;
	    exit(1) ;
	}
    b = parms ;
/*
printf(" after b=parms MA:  P Q t0  b[0] parms[0] =  %i   %i    %i    %e       %e     \n",P,Q,t0,b[0],parms[0]);
*/
for ( i = 0 ; i < S ; i++ )
	c[i] = parms + Q + R*i ;
}

//VP
int initialize(char* ininame)
{
   char str[100];
   int nCol = 0;
   memset(cols, 0, sizeof(int) * MAXCOLS);

   FILE* fp = fopen(ininame, "r");

   int i = 0;
   while (fgets(str, 100, fp) != NULL)
   {
       i++;
      char str1[100], str2[100];
      sscanf(str, "%s %s", str1, str2);
//      printf("XX%sXX\n", str);

      if (strcmp(str1, "INP") == 0)
      {
         strcpy(inputfile, str2);
//	     printf("INP %s\n", inputfile);

      }
      else if (strcmp(str1, "OUTP") == 0)
      {
         strcpy(outputfile, str2);
//	     printf("OUTP %s\n", outputfile);

      }
/*
      else if (strcmp(str1, "T") == 0)
      {
         T = atoi(str2);
//	     printf("T %d\n", T);

      }
      else if (strcmp(str1, "T0") == 0)
      {
         T0 = atoi(str2);
//	     printf("T0 %d\n", T0);

      }
      else if (strcmp(str1, "T1") == 0)
      {
         T1 = atoi(str2);
//	     printf("T1 %d\n", T1);

      }
      else if (strcmp(str1, "QMAX") == 0)
      {
         Qmax = atoi(str2);
//	     printf("QMAX %d\n", Qmax);

      }
      else if (strcmp(str1, "PMAX") == 0)
      {
         Pmax = atoi(str2);
//	     printf("PMAX %d\n", Pmax);

      }
*/
/*
     else if (strcmp(str1, "QO") == 0)
      {
         Qo = atoi(str2);
//	     printf("Qo %d\n", Qo);

      }

      else if (strcmp(str1, "PO") == 0)
      {
         Po = atoi(str2);
//	     printf("Po %d\n", Po);

      }
*/
      else if (strcmp(str1, "COL") == 0)
      {
         int iCol = atoi(str2);
         iCol--;

         if (iCol >= MAXCOLS || iCol < 0)
	 {
            printf("Wrong column: %s\n", str2);
	    continue;
	 }
//	 else
//	     printf("COL %d\n", iCol);

         cols[iCol] = 1;
         nCol++;
      }
      else

        printf("Check initialization strings!
 %s\n", str);

  }

//   printf("Lines in INI file: %d\n", i);
/*
   M = nCol;
   T *= nCol;
   T0 *= nCol;
   T1 *= nCol;
   Pmax *= nCol;
   Po *= nCol;
*/
//   printf("T: %d, T0 %d, T1 %d, M %d\n", T, T0,
//	  T1, M);

   A = new double[T][DM];
   B = new double[T];
   z = new double[T];
//   z0 = new double[T];
//   z1 = new double[T];
//   z2 = new double[T];
//   z3 = new double[T];
//   zp = new double[T];

//   printf("NIT: %d\n", NIT);
/*
   fm = new double[NIT];
   xn = new double*[NIT];
*/
   // VP double fm[NIT], *xn[NIT];

   return 1; //we might return 0 on error
}

void shutdown()
{
   delete [] A;
   delete [] B;
   delete [] z;
//   delete [] z0;
//   delete [] z1;
//   delete [] z2;
//   delete [] z3;
//   delete [] zp;

   delete [] zInp;

   delete [] fm;
   delete [] xn;
}