#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "def.h"
#define MAXLNSZ 100
#define MAXFMT  100
float mean(float*ary,int n){
      float ret=0;
      int i=0;
      for(i=0;i<n;i++){
         ret += ary[i];      
      }
      
      return ret/n;
}

float var(float*ary,int n){
      float sumx = 0;
      float sumx2= 0;
      int i = 0;
      for(i=0;i<n;i++){
         sumx += ary[i];            
         sumx2+= ary[i]*ary[i];
      }                     
      
      return (sumx2-sumx*sumx/n)/(n-1);
}

float sumsquare(float*ary,int n){
      float ret = 0;
      int i=0;
      for(i = 0; i<n; i++){
            ret += ary[i]*ary[i];            
      }
      
      return ret;
}

//read one integer per line from text file
int read_int(FILE*f){
    char buf[MAXLNSZ];
    int  nret;
    
    
    if(!f){
           printf("invalid file handle.\n");
           exit(1);                 
    }
    
    if(fgets(buf,MAXLNSZ,f) == NULL){
           printf("Error reading integer.\n");
           exit(1);                               
    }   
    
    if(!sscanf(buf,"%d\n", &nret)){
        printf("Error in sscanf.\n");
        exit(1);
    }

    return nret;
}

//read multiple floating points number per line
//currently only 2 covariates
//seperated by a space
void read_float(FILE*f,float*cov,int n){
    char buf[MAXLNSZ];
    char *substr=NULL;
    int i;
    if(!f){
           printf("invalid file handle.\n");
           exit(1);                 
    }
    
    if(fgets(buf,MAXLNSZ,f) == NULL){
           printf("Error reading float.\n");
           exit(1);                               
    }   
    
    substr = strtok(buf," ");
    
    if(!substr) {
          printf("Error reading float.\n");
          exit(1);
    }
    
    cov[0]=atof(substr);

    for(i=0;i<n-1;i++){
        substr=strtok(NULL," ");
        cov[i+1]= atof(substr);
    }
   
}
//read data from file
int read_data(char*file, Data*ddata){
    FILE *infile;
    char buf[MAXLNSZ];
    
    /* Open the file.  If NULL is returned there was an error */
    if((infile = fopen(file, "r")) == NULL) {
      printf("Error Opening data File.\n");
      exit(1);
    }

    
    /*start reading from file*/
    int i,j,k;
    int nsite=read_int(infile);
    
    //assigns site info
    (*ddata).nst = nsite;
    (*ddata).st = malloc(sizeof(Site)*nsite);
    
    for(i = 0;i < nsite; i++){         
          int npls=read_int(infile);
          
          //pools info in site i
          (*ddata).st[i].npl=npls;
          (*ddata).st[i].pls=malloc(sizeof(Pool)*npls);
          
          for(j=0; j<npls; j++){                   
                    //subjects info in pool ij
                    int nsubs=read_int(infile);                   
                    (*ddata).st[i].pls[j].nsub = nsubs;
                    (*ddata).st[i].pls[j].tij  = read_int(infile);
                    (*ddata).st[i].pls[j].xijk = malloc(sizeof(float)*nsubs);
                    read_float(infile,(*ddata).st[i].pls[j].xijk,nsubs);
          }
    }//end loop i
    
    //print data
   /* 
    printf("%d\n",(*ddata).nst);
    for(i=0;i<(*ddata).nst;i++){
        printf("%d\n",(*ddata).st[i].npl);
        for(j=0;j<(*ddata).st[i].npl;j++){
          printf("%d\n", (*ddata).st[i].pls[j].nsub);
          printf("%d\n",(*ddata).st[i].pls[j].tij);
          for(k=0;k<(*ddata).st[i].pls[j].nsub;k++){
             printf("%f",(*ddata).st[i].pls[j].xijk[k]);             
          }
          printf("\n");
        }
    }
    
    */
    fclose(infile);
    return 0;
}

//initialize random effect 
void init_rf(Reff*rf, Strf*u0, int nsite){
     int i,j,k;
     
     for(i=0;i<(*rf).M;i++){
        (*rf).u[i].nst = nsite;
        (*rf).u[i].uh  = malloc(sizeof(float)*nsite);
        if((*rf).u[i].uh == NULL){printf("memory allocation failed.\n"); exit(1);}
     }
     
     for(j=0;j<nsite;j++){
       (*u0).uh[j]=0;
    }    
}

//free memory for random effect for one run
void free_strf(Strf*u){      
    free((*u).uh);
}

void free_rf(Reff*rf){
     int i,k;
     for(i=0;i<(*rf).M;i++){
        free_strf((*rf).u+i);
     }   
     
     free((*rf).u);
     
     free(rf);
}

float square(float x){
      return x*x;
}

//inversion of 2*2 matrix
//inverted matrix is stored in x
void Invert22(float x[2][2])
{
     float a0=x[0][0],a1=x[0][1],a2=x[1][0],a3=x[1][1];
     x[0][0]=a3/(a0*a3-a1*a2);
     x[1][0]=-(a2/a3)*x[0][0];
     x[1][1]=a0/(a0*a3-a1*a2);
     x[0][1]=-(a1/a0)*x[1][1];   
}

//scalar multiplication of matrix of dim m*n
//the first argument is &mat[0][0]
/*void mult(float **mat,int m,int n,float c){
     int i,j;
     for(i=0;i<m;i++){
          for(j=0;j<n;j++)
             mat[i][j] *= c;
     }
}*/

//one dimensional matrix subtraction
void vec_sub(float*mat1,float*mat2,float *ret,int n){
        int i;
        
        for(i=0;i<n;i++){
           *(ret+i)=*(mat1+i)-*(mat2+i);       
        }
}

//store added value in ret
//ret might be equal to mat1 or mat2;
void vec_add(float*mat1,float*mat2,float*ret,int n){
     int i;
     
     for(i=0;i<n;i++){
        *(ret+i) = *(mat1+i) + *(mat2+i);     
     }

}
//returns the element with maximum abs in one dimensional vector
float vec_max(float*mat,int n){
        int i;
        float ret = fabs(*mat);
        float temp;
        
        for(i=1;i<n;i++){
             temp = fabs(*(mat+i));
             if(temp>ret){
                 ret = temp;
             }
        }

        return ret;
}

//write 5 parameters to a line
void writeline(char*file,float fpar[NFPAR], float rpar){
  FILE*out;
  
  out=fopen(file,"a");
  
  //if not exist then create
  if(!out){
    out = fopen(file,"w");
    if(!out) {
       printf("creating output file failed.\n");
       exit(1);
    }
  }
  
  fprintf(out, "%f %f %f\n", fpar[0],fpar[1],rpar);
  fclose(out);
}

int check_cvg(float fpar[NFPAR],float rpar, 
              float del_fpar[NFPAR],float del_rpar,
              float tol){
       int j;
       float temp1[NFPAR];
       int ret=0;
       
       for(j=0;j<NFPAR;j++){
              temp1[j] = del_fpar[j]/fpar[j];
       }
             
       if((vec_max(temp1,NFPAR)<tol)&&(fabs(del_rpar/rpar)<tol))ret = 1;       
       else ret=0;
       
       return ret;

}