/*
 ***************************************************************************
 * This package contains code for the paper "Mobile Money in Tanzania" by  *
 *                                                                         *
 * Nicholas Economides                                                     *
 * Stern School of Business, NYU; NET Institute; economides@stern.nyu.edu  *
 *                                                                         *
 * and                                                                     *
 *                                                                         *
 * Przemyslaw Jeziorski                                                    *
 * Haas School of Business, UC Berkeley; przemekj@haas.berkeley.edu        *
 *                                                                         *
 * Copyright by Przemyslaw Jeziorski, 2016                                 *
 *                                                                         *
 * This program is free software: you can redistribute it and/or modify    *
 * it under the terms of the GNU General Public License as published by    *
 * the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                     *
 *                                                                         *
 * This program is distributed in the hope that it will be useful,         *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
 * GNU General Public License for more details.                            *
 *                                                                         *
 * You should have received a copy of the GNU General Public License       *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 ***************************************************************************
*/


#include<probit.h>
#include<time.h>

mxArray *readMatlab(MATFile *mfp, const char *var_name);

#define COUT_DEBUG

inline long fastint(long *g);

int main(int argc, char *argv[]) {
  time_t start,end;
  double dif;
  MATFile *mfp;
  FILE *myfile;
  mxArray *transfersMx, *distanceMx, *parametersMx, *balanceMx;
  double *parameters;
  int size;
  long *mySeed = (long *) mxMalloc(4*sizeof(long));
  pthread_t thread_id[SMP];
  pthread_attr_t attr;
  ThreadStructure *input;
  int jump;
  int j,i,i1,k,p,P;
  double *momentsOut, *momentsGradientOut;
  double *momentsDataOut, *varcovarOut, *welfareOut;
  mxArray *momentsDataOutMx, *momentsOutMx, *momentsGradientOutMx, *varcovarOutMx, *welfareOutMx, *urbanMx, *price_counterMx;
  long genSeed;
  int jac;
  

  myfile = fopen ("output.txt","a");

  if ((mfp = matOpen(argv[1], "r")) == NULL) {
    fprintf(stderr, "Cannot open MAT-file %s\n", argv[1]);
    #ifdef DEBUG
      fprintf(myfile, "Cannot open MAT-file %s\n", argv[1]);
    #endif
    #ifdef COUT_DEBUG
      printf("Cannot open MAT-file %s\n", argv[1]);
    #endif
    exit(1);
  } else {
    #ifdef DEBUG
      fprintf(myfile, "Opening MAT-file %s\n", argv[1]);
    #endif
    #ifdef COUT_DEBUG
      printf("Opening MAT-file %s\n", argv[1]);
    #endif
  }

  transfersMx = readMatlab(mfp, "transfers");
  transfersData = mxGetPr(transfersMx);
  distanceMx = readMatlab(mfp, "distance");
  distanceData = mxGetPr(distanceMx);
  balanceMx = readMatlab(mfp, "balance");
  balanceData = mxGetPr(balanceMx);
  urbanMx = readMatlab(mfp, "urban");
  urbanData = mxGetPr(urbanMx);


  size = mxGetN(transfersMx);
  jump=size/SMP+1;

  if (matClose(mfp) != 0) {
    printf("Error closing file %s\n", argv[1]);
    exit(EXIT_FAILURE);
  }

  if ((mfp = matOpen(argv[2], "r")) == NULL) {
    fprintf(stderr, "Cannot open MAT-file %s\n", argv[2]);
    #ifdef DEBUG
      fprintf(myfile, "Cannot open MAT-file %s\n", argv[2]);
    #endif
    #ifdef COUT_DEBUG
      printf("Cannot open MAT-file %s\n", argv[2]);
    #endif
    exit(1);
  } else {
    #ifdef DEBUG
      fprintf(myfile, "Opening MAT-file %s\n", argv[2]);
    #endif
    #ifdef COUT_DEBUG
      printf("Opening MAT-file %s\n", argv[2]);
    #endif
  }

  parametersMx = readMatlab(mfp, "parameters");
  parameters = mxGetPr(parametersMx);

  price_before = mxGetPr(readMatlab(mfp, "price_before"));
  price_after = mxGetPr(readMatlab(mfp, "price_after"));
  amounts = mxGetPr(readMatlab(mfp, "amounts"));
  jac = mxGetPr(readMatlab(mfp, "jac"))[0];

  price_counterMx = readMatlab(mfp, "price_counter");
  if(price_counterMx==NULL) {
    price_counter=NULL;
    counterfactual=0;
  } else {
    price_counter=mxGetPr(price_counterMx);
    counterfactual=mxGetPr(readMatlab(mfp, "counterfactual"))[0];
  }

  if (matClose(mfp) != 0) {
    printf("Error closing file %s\n", argv[2]);
    exit(EXIT_FAILURE);
  }

  if ((mfp = matOpen(argv[4], "r")) == NULL) {
    fprintf(stderr, "Cannot open MAT-file %s\n", argv[4]);
    #ifdef DEBUG
      fprintf(myfile, "Cannot open MAT-file %s\n", argv[4]);
    #endif
    #ifdef COUT_DEBUG
      printf("Cannot open MAT-file %s\n", argv[4]);
    #endif
    exit(1);
  } else {
    #ifdef DEBUG
      fprintf(myfile, "Opening MAT-file %s\n", argv[4]);
    #endif
    #ifdef COUT_DEBUG
      printf("Opening MAT-file %s\n", argv[4]);
    #endif
  }

  ext_before = readDoubleMatrix(readMatlab(mfp, "ext_before"));
  ext_after = readDoubleMatrix(readMatlab(mfp, "ext_after"));

  if (matClose(mfp) != 0) {
    printf("Error closing file %s\n", argv[4]);
    exit(EXIT_FAILURE);
  }

  if(jac==0) {
    P=1;
  } else {
    P=PARAMETERS+1;
  }
  time(&start);

  momentsOutMx = mxCreateDoubleMatrix(MOMENTS,1,mxREAL);
  momentsOut=mxGetPr(momentsOutMx);
  momentsDataOutMx = mxCreateDoubleMatrix(MOMENTS,1,mxREAL);
  momentsDataOut=mxGetPr(momentsDataOutMx);

  momentsGradientOutMx = mxCreateDoubleMatrix(MOMENTS,PARAMETERS,mxREAL);
  momentsGradientOut = mxGetPr(momentsGradientOutMx);

  varcovarOutMx = mxCreateDoubleMatrix(MOMENTS,MOMENTS,mxREAL);
  varcovarOut = mxGetPr(varcovarOutMx);

  welfareOutMx = mxCreateDoubleMatrix(1,1,mxREAL);
  welfareOut = mxGetPr(welfareOutMx);

  for(p=0;p<P;++p) {
    printf("%d\n",p);
    if(p>0) {
      parameters[p-1]+=0.001;
    }

    genSeed = 687162341;
    for(i=0;i<SMP;++i) {
      dsfmt_init_gen_rand(dsfmt+i, fastint(&genSeed));
    }

    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
  
    input = (ThreadStructure *) malloc(SMP*sizeof(ThreadStructure));
    for (j=0;j<SMP;j++) {
      input[j].thread=j;
      input[j].start=jump*j;
      input[j].stop=MIN(jump*(j+1),size);
      for (i=0;i<PARAMETERS;++i) {
        input[j].parameters[i]=parameters[i];
      }
      input[j].jac=p;
      pthread_create(&thread_id[j], &attr, moments, (void *) &input[j] );
    }
    pthread_attr_destroy(&attr);
    for (j=0;j<SMP;j++) {
      pthread_join(thread_id[j], NULL);
    }
    if(p>0) {
      parameters[p-1]-=0.001;
      for (i=0;i<MOMENTS;++i) {
        momentsGradientOut[i+(p-1)*MOMENTS]=0;   
      }
      for (j=0;j<SMP;j++) {
        for (i=0;i<MOMENTS;++i) {
          momentsGradientOut[i+(p-1)*MOMENTS]+=input[j].moments[i]/size;
        }
      } 
    } else {
      for (i=0;i<MOMENTS;++i) {
        momentsOut[i]=0;   
        momentsDataOut[i]=0;
      }
      for (j=0;j<SMP;j++) {
        for (i=0;i<MOMENTS;++i) {
          momentsOut[i]+=input[j].moments[i]/size;
          momentsDataOut[i]+=input[j].moments_data[i]/size;
          for (i1=0;i1<=i;++i1) {
            varcovarOut[i1*MOMENTS+i]+=input[j].varcovar[i][i1]/size;
          }
          for (i1=i+1;i1<MOMENTS;++i1) {
            varcovarOut[i1*MOMENTS+i]+=input[j].varcovar[i1][i]/size;
          }
        }
        welfareOut[0]+=input[j].welfare/size;
      }
    }
  }

  time(&end);
  dif = difftime(end,start);
  printf ("\nCalculations took %.2lf seconds to run.\n", dif );

  if ((mfp = matOpen(argv[3], "w")) == NULL) {
    fprintf(stderr, "Cannot open MAT-file %s\n", argv[3]);
    #ifdef DEBUG
      fprintf(myfile, "Cannot open MAT-file %s\n", argv[3]);
    #endif
    #ifdef COUT_DEBUG
      printf("Cannot open MAT-file %s\n", argv[3]);
    #endif
    exit(1);
  } else {
    #ifdef DEBUG
      fprintf(myfile, "Opening MAT-file %s\n", argv[3]);
    #endif
    #ifdef COUT_DEBUG
      printf("Opening MAT-file %s\n", argv[3]);
    #endif
  }

  matPutVariable(mfp, "moments", momentsOutMx);
  if(jac>0) matPutVariable(mfp, "momentsGradient", momentsGradientOutMx);
  matPutVariable(mfp, "momentsData", momentsDataOutMx);
  matPutVariable(mfp, "varcovar", varcovarOutMx);
  matPutVariable(mfp, "welfare", welfareOutMx);

  if (matClose(mfp) != 0) {
    printf("Error closing file %s\n", argv[3]);
    exit(EXIT_FAILURE);
  }


  fclose(myfile);
}

mxArray *readMatlab(MATFile *mfp, const char *var_name) {
  mxArray *array_ptr;
  int *from, *to;

  if ((array_ptr = matGetVariable(mfp, var_name)) == NULL) {
    #ifdef COUT_DEBUG
    printf("Cannot read matlab variable %s\n", var_name);
      #endif
      #ifdef DEBUG
        fprintf(myfile, "Cannot read matlab variable %s\n", var_name);
      #endif
    return NULL;
  } 

  if ((mxGetNumberOfElements(array_ptr) == 1) && mxIsNumeric(array_ptr)) {
    #ifdef COUT_DEBUG
      printf("Reading matlab variable %s = %f\n", var_name, mxGetScalar(array_ptr));
    #endif
    #ifdef DEBUG
      fprintf(myfile, "Reading matlab variable %s = %f\n", var_name, mxGetScalar(array_ptr));
    #endif
  } else {
    #ifdef COUT_DEBUG
      printf("Reading matlab variable %s (", var_name);
    #endif
    #ifdef DEBUG
      fprintf(myfile, "Reading matlab variable %s (", var_name);
    #endif
    from = (int *) mxGetDimensions(array_ptr);
    to = from+mxGetNumberOfDimensions(array_ptr)-1;
    for (; from<to; from++) {
      #ifdef COUT_DEBUG
        printf("%dx", *from);
      #endif
      #ifdef DEBUG
        fprintf(myfile, "%dx", *from);
      #endif
    }
    #ifdef COUT_DEBUG
      printf("%d)\n", *from);
    #endif
    #ifdef DEBUG
      fprintf(myfile, "%d)\n", *from);
    #endif
  }

 
 return array_ptr;
}

inline long fastint(long *g) {
  *g = (214013*(*g)+2531011);

  return (*g>>32)& 0x7FFFFFFF;
}

double **readDoubleMatrix(mxArray *in) {
  int columns=mxGetN(in);
  int rows=mxGetM(in);
  double *inPr=mxGetPr(in);
  double **out = (double **) mxMalloc(rows*sizeof(double *));
  int i,j,k=0;

  for(i=0;i<rows;++i) {
    out[i]=(double *) mxMalloc(columns*sizeof(double));
  }

  for(j=0;j<columns;++j) {
    for(i=0;i<rows;++i) {
      out[i][j]=inPr[k];
      k++;
    }
  }

  mxDestroyArray(in);

  return out;
}
