/*
 ***************************************************************************
 * 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>

void vectorNormal(dsfmt_t *g, int n, double *output);
void cholesky(double *A, double *L, int n);

void *moments(void *dataIn) {
  int i,j,t,d,r,m,f,m1,z;
  int choice;
  double umax, u;
  double draws[36], draws2[38], choice_draws[33];
  double u0[33][14];
  double first_half[33], persistent[33];
  double distance;
  double aggr;
  double balance, balance_log;
  int balance_idx;

  double moments_obs[MOMENTS];
  double moments_data_obs[MOMENTS];
  int firstn, secondn;
  double omega[1089], L[1089];
  const double bound = sqrt(1/(4*pow(cos(PI/(33+1)),2))-1e-6);
  double rho1;
  double random_coef[3];
  double sum = 0;
  int j1;
  double alpha;
  double distance_nolog;
  double distance_log;

  double price_rural[33];
  double price_urban[33];

  ThreadStructure *data = (ThreadStructure *) dataIn;

  if(counterfactual==1) {
    for(i=0;i<33;++i) {
      price_rural[i]=price_before[i]+price_counter[0];
      price_urban[i]=price_before[i]+price_counter[1];
    }
  } 


  for(m=0;m<MOMENTS;++m) {
    data->moments[m]=0;
    data->moments_data[m]=0;

    if(data->jac==0) {
      for(m1=0;m1<=m;++m1) {
        data->varcovar[m][m1]=0;
      }
    }
  }
  data->welfare=0;

  for(f=0;f<33;++f) {
    u0[f][0]=data->parameters[PARAM_FIXED+f];
    for(t=1;t<7;++t) {
      u0[f][t]=data->parameters[PARAM_FIXED+f]+data->parameters[PARAM_TIME+t-1];
    }
    for(t=7;t<14;++t) {
      u0[f][t]=data->parameters[PARAM_FIXED+f]+data->parameters[PARAM_TIME+t-1];
    }
  }

  // Prepare omega
  rho1=2*atan(data->parameters[PARAM_EPSILON_CORRELATION])/PI*bound;
  for(i=0;i<1089;++i) omega[i]=0;
  omega[0]=1;omega[1]=rho1;
  for(i=1;i<32;++i) {
    omega[34*i-1]=rho1;
    omega[34*i]=1;
    omega[34*i+1]=rho1;
  }
  omega[1087]=rho1;omega[1088]=1;
  cholesky(omega, L, 33);
  
  j=data->start*14*3;
  j1=data->start*14*3;
  for(i=data->start;i<data->stop;++i) {
    for(m=0;m<MOMENTS;++m) {
      moments_obs[m]=0;
      moments_data_obs[m]=0;
    }
    for(r=0;r<RSIM;++r) {
      firstn=0; secondn=0;
      vectorNormal(dsfmt+data->thread, 38, draws2);
      aggr=draws2[33]*data->parameters[PARAM_AGGR_STDDEV];   
      for(f=0;f<33;++f) {
        persistent[f]=draws2[f]*data->parameters[PARAM_PERSISTENT_STDDEV];
        first_half[f]=0;
      }
      // Random coefficients
      random_coef[0]=data->parameters[PARAM_PRICE_TRANSFER_MEAN]-data->parameters[PARAM_VARCOVAR_L11]*draws2[34];
      random_coef[1]=MIN(data->parameters[PARAM_PRICE_CASHOUT_MEAN]-data->parameters[PARAM_VARCOVAR_L21]*draws2[34]-data->parameters[PARAM_VARCOVAR_L22]*draws2[35],0);
      random_coef[2]=data->parameters[PARAM_VARCOVAR_L31]*draws2[34]+data->parameters[PARAM_VARCOVAR_L32]*draws2[35]+data->parameters[PARAM_VARCOVAR_L33]*draws2[36];
      sum += random_coef[1]*random_coef[1];

      for(t=0;t<7;++t) {
        for(d=0;d<3;++d) {
          vectorNormal(dsfmt+data->thread, 36, draws);
          choice_draws[0]=draws[0];
          z=34;
          for(f=1;f<33;++f) {
            choice_draws[f]=draws[f-1]*L[z-1]+draws[f]*L[z];
            z+=34;
          }
          distance=exp(data->parameters[PARAM_DISTANCE_MEAN]+data->parameters[PARAM_DISTANCE_STDDEV]*draws[34]);
          balance=exp(data->parameters[PARAM_BALANCE_BEFORE_MEAN]+data->parameters[PARAM_BALANCE_BEFORE_STDDEV]*draws[35]);
          balance_log=log(balance+1);
          distance_log=log(distance+1);

          balance_idx=(balance<BALANCE_UPPER) ? (int) balance/BALANCE_STEP : 1000;

          if(counterfactual==1) {
            if(urbanData[j1]>0.5) {
              price_before=price_urban;
            } else {
              price_before=price_rural;
            }
          } else if(counterfactual==2) {
            if((distance<price_counter[0]) || (distance>=price_counter[1])) {
              continue;
            }
          } 

          choice=33;
          umax=data->parameters[PARAM_GAMMA]*distance_log+draws[33];
          alpha=MIN(random_coef[0]+(urbanData[j1]>0.5)*data->parameters[65],0);

          for(f=0;f<33;++f) {
            u=u0[f][t]+amounts[f]*random_coef[2]+persistent[f]+choice_draws[f]+price_before[f]*alpha+random_coef[1]*ext_before[balance_idx][f]+aggr+(urbanData[j1]>0.5)*data->parameters[64];
            if(u>umax) {
              choice=f;
              umax=u;
            }
          }

          data->welfare+=umax;

          if(choice<33) {
            moments_obs[choice]++;
            moments_obs[66+t]++;

            firstn++;

            moments_obs[80]+=distance;
            moments_obs[81]+=distance*distance;

            moments_obs[94]+=balance_log;
            moments_obs[95]+=balance_log*balance_log;

            moments_obs[98]+=distance*choice;
            moments_obs[99]+=balance_log*choice;

            first_half[choice]=1;

            moments_obs[102]+=(urbanData[j1]>0.5);
          }
          j1++;
        }
      }

      for(t=7;t<14;++t) { 
        for(d=0;d<3;++d) {
          vectorNormal(dsfmt+data->thread, 36, draws);
          choice_draws[0]=draws[0];
          z=34;
          for(f=1;f<33;++f) {
            choice_draws[f]=draws[f-1]*L[z-1]+draws[f]*L[z];
            z+=34;
          }

          distance=exp(data->parameters[PARAM_DISTANCE_MEAN]+data->parameters[PARAM_DISTANCE_STDDEV]*draws[34]);
          balance=exp(data->parameters[PARAM_BALANCE_AFTER_MEAN]+data->parameters[PARAM_BALANCE_AFTER_STDDEV]*draws[35]);
          balance_log=log(balance+1);
          distance_log=log(distance+1);

          balance_idx=(balance<BALANCE_UPPER) ? (int) balance/BALANCE_STEP : 1000;

          if(counterfactual==1) {
            if(urbanData[j1]>0.5) {
              price_before=price_urban;
            } else {
              price_before=price_rural;
            } 
          } else if(counterfactual==2) {
            if((distance<price_counter[0]) || (distance>=price_counter[1])) {
              continue;
            }
          } 

          choice=33;
          umax=data->parameters[PARAM_GAMMA]*distance_log+draws[33];
          alpha=MIN(random_coef[0]+(urbanData[j1]>0.5)*data->parameters[65],0);

          for(f=0;f<33;++f) {
            u=u0[f][t]+amounts[f]*random_coef[2]+persistent[f]+choice_draws[f]+price_after[f]*alpha+random_coef[1]*ext_after[balance_idx][f]+aggr+(urbanData[j1]>0.5)*data->parameters[64];
            if(u>umax) {
              choice=f;
              umax=u;
            }
          }
          if(choice<33) {
            moments_obs[33+choice]++;
            moments_obs[66+t]++;

            secondn++;

            moments_obs[82]+=distance;
            moments_obs[83]+=distance*distance;

            moments_obs[96]+=balance_log;
            moments_obs[97]+=balance_log*balance_log;

            moments_obs[100]+=distance*choice;
            moments_obs[101]+=balance_log*choice;

            moments_obs[103]+=(urbanData[j1]>0.5);

            if(first_half[choice]==1) {
              moments_obs[93]++;
            }
          }
          j1++;
        }
      }

      moments_obs[84+MIN(firstn,2)*3+MIN(secondn,2)]++;
    }

    for(f=0;f<33;f++) {
      first_half[f]=0;
    }
    firstn=0;secondn=0;
    for(t=0;t<7;++t) {
      for(d=0;d<3;++d) {
        if(transfersData[j]<33) {
          balance_log=log(balanceData[j]+1);

          moments_data_obs[(int) transfersData[j]]++;
          moments_data_obs[66+t]++;

          firstn++;

          moments_data_obs[80]+=distanceData[j];
          moments_data_obs[81]+=distanceData[j]*distanceData[j];

          moments_data_obs[94]+=balance_log;
          moments_data_obs[95]+=balance_log*balance_log;

          moments_data_obs[98]+=transfersData[j]*distanceData[j];
          moments_data_obs[99]+=transfersData[j]*balance_log;

          moments_data_obs[102]+=(urbanData[j]>0.5);

          first_half[(int) transfersData[j]]=1;
        }
        j++;
      }
    }
    for(t=7;t<14;++t) {  
      for(d=0;d<3;++d) {
        if(transfersData[j]<33) {
          balance_log=log(balanceData[j]+1);

          moments_data_obs[33+(int) transfersData[j]]++;
          moments_data_obs[66+t]++;

          secondn++;

          moments_data_obs[82]+=distanceData[j];
          moments_data_obs[83]+=distanceData[j]*distanceData[j];

          moments_data_obs[96]+=balance_log;
          moments_data_obs[97]+=balance_log*balance_log;

          moments_data_obs[100]+=transfersData[j]*distanceData[j];
          moments_data_obs[101]+=transfersData[j]*balance_log;

          moments_data_obs[103]+=(urbanData[j]>0.5);

          if(first_half[(int) transfersData[j]]==1) {
            moments_data_obs[93]++;
          }
        }
        j++;
      } 
    }

    moments_data_obs[84+MIN(firstn,2)*3+MIN(secondn,2)]++;

    for(m=0;m<MOMENTS;++m) {
      moments_obs[m]/=RSIM;
      data->moments[m]+=moments_obs[m];
      data->moments_data[m]+=moments_data_obs[m];
    }
    data->welfare/=RSIM;

    if(data->jac==0) {
      for(m=0;m<MOMENTS;++m) {
        for(m1=0;m1<=m;++m1) {
          data->varcovar[m][m1]+=(moments_obs[m]-moments_data_obs[m])*(moments_obs[m1]-moments_data_obs[m1]);
        }
      }
    }
  }
}

void vectorNormal(dsfmt_t *g, int n, double *output) {
  int i;
  float x1, x2, w, y1, y2;
  for(i=0;i<n;i+=2) {
    do {
      x1 = 2.0 * dsfmt_genrand_open_open(g) - 1.0;
      x2 = 2.0 * dsfmt_genrand_open_open(g) - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
    w = sqrt( (-2.0 * log( w ) ) / w );
    output[i] = x1 * w;
    output[i+1] = x2 * w;
  }
}

void cholesky(double *A, double *L, int n) {
  int i, j, k;
  double s;
  for(i = 0; i < n*n; ++i) {
    L[i]=0;
  }
  for(i = 0; i < n; i++)
    for(j = 0; j < (i+1); j++) {
      s = 0;
      for (k = 0; k < j; k++)
        s += L[i * n + k] * L[j * n + k];
      L[i * n + j] = (i == j) ?
        sqrt(A[i * n + i] - s) :
        (1.0 / L[j * n + j] * (A[i * n + j] - s));
    }
}
