#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "def.h"

#define MAXPATHLEN 100
#define MAXITER 1000

int main(int argc, char *argv[])
{
  int read_data(char*file, Data*ddata);    
  float fuij(float fpar[NFPAR],Pool pl, float ui);
  void Metroplis(float fpar[NFPAR], float rpar,Data*data, int nburn, int nskip, Reff*rf, Strf*u0, long*idum);
  void init_rf(Reff*rf, Strf*u0, int nsite);
  void score_betaij(float fpar[NFPAR],Pool*pl, float ui, float dfij[NFPAR],float ddfij[NFPAR][NFPAR]);
  void parm_update(Reff*rf,
                 Data*data,
                 float oldfpar[NFPAR],
                 float delta_fpar[NFPAR],
                 float *newrpar);
  void vec_sub(float*mat1,float*mat2,float *ret,int n);
  void vec_add(float*mat1,float*mat2,float*ret,int n);
  float vec_max(float*mat,int n);
  void free_rf(Reff*rf);
  float square(float x);
  void writeline(char*file,float fpar[NFPAR], float rpar);
  int check_cvg(float fpar[NFPAR],float rpar, 
              float del_fpar[NFPAR],float del_rpar,
              float tol);
  
   char  path[MAXPATHLEN];
   char  out[MAXPATHLEN];
   float fpar[NFPAR];
   float rpar;
   int nburn;
   int nskip; 
   int size[5];
   float tol;

   printf("datafile:\n");
   scanf("%s", path);
   printf("outputfile:\n");
   scanf("%s", out);
   printf("nburn:\n");
   scanf("%d",&nburn);
   printf("nskip:\n");
   scanf("%d",&nskip);
   printf("par(2 fixed 1 random):\n");
   scanf("%f %f %f" , &fpar[0],&fpar[1],&rpar);
   printf("monte-carlo size:\n");
   scanf("%d %d %d %d %d", &size[0],&size[1],&size[2],&size[3],&size[4]);
   printf("tol:\n");   
   scanf("%f",&tol);
  
  Data ddata;
  
  //allocate memory and read data from file
  read_data(path,&ddata);
  
  int nsite = ddata.nst;
  float delta_fpar[NFPAR],delta_rpar,rnewpar;
  int  M,i;

  int iter = 1;
  int converge=0;
  long idum=-1090;    

  //starting values
  Strf u0={nsite, malloc(sizeof(float)*nsite)};
  
  Reff rndeff;  
  while((iter < MAXITER)&& (converge==0)){
      if(iter<30)M=size[0];
      else if(iter<50)M=size[1];
      else if(iter<70)M=size[2];
      else if(iter<100)M=size[3];
      else M=size[4];


      //each iteration has to re=allocate 
      //random effect mem
      rndeff.M = M;
      rndeff.u = malloc(sizeof(Strf)*M);

      //allocate memories for random effect,
      //including resetting starting values 
      //u0 to (0,0)
      init_rf(&rndeff, &u0, nsite);
      
      //change seed ?
      idum += -12;
                 
      Metroplis(fpar,rpar,&ddata,nburn,nskip,&rndeff,&u0,&idum);
      parm_update(&rndeff,&ddata,fpar,delta_fpar,&rnewpar);
      delta_rpar=rpar-rnewpar;
      iter++;
                                
      vec_add(fpar,delta_fpar,fpar,NFPAR);
      rpar=rnewpar;
      
      // check for convergence, save up to 5 iteration
      if(check_cvg(fpar,rpar, 
               delta_fpar,delta_rpar,
               tol) == 1){
              converge = 1;
              printf("converges");
            }
           
      writeline(out,fpar,rpar);
      //printf("fixed effect:%f %f \n",fpar[0],fpar[1]);
      //printf("random effect:%f %f %f\n",rnewpar[0],rnewpar[1],rnewpar[2]);
      
      //release mem for rf
      free_rf(&rndeff);
  }

  return 0;
}