#ifndef SIM_HPP
#define SIM_HPP

#include <gsl/gsl_rng.h>
#include <complex>
#include <math.h>
#include <sys/time.h>

#include "latlib/neigh.h"

#define EPSILONU 1
#define EPSILONPHI 0.5

class sim : public o815::sim {
public:
  sim(o815 *_O815);
  unsigned int lsize4;
  neigh *nb;
  complex<double> *U, *phi;
  double kappa[2], lambda[2], beta;
  void _newParas();
  
private:
  void _resetConfig();
  void _makeSweep();
  gsl_rng* rangsl;
  double rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi);
  double rhoU(const int& x0, const int& nu0, const complex<double>& candU);
  int updatePhi(const int& iphi, const int& x0);
  int updateU(const int& x0, const int& nu0);
};

sim::sim(o815 *_O815) : o815::sim( _O815, 
				   sizeof(complex<double>)*
				   _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2+4) ) {

  struct timeval tv;

  lsize4 = _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1];

  nb = new neigh(4, _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[1]);

  phi = (complex<double>*)confMem;
  U = (complex<double>*)(confMem + sizeof(complex<double>)*lsize4*2);

  gettimeofday(&tv, NULL);
  rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
  gsl_rng_set(rangsl, 1000000 * tv.tv_sec + tv.tv_usec);
}

void sim::_makeSweep() {  
  for( int ix=0; ix<lsize4; ix++ ) {
    for( int inu=0; inu<4; inu++) updateU(ix, inu);
    for( int iphi=0; iphi<2; iphi++) updatePhi(iphi, ix);
  }
}

void sim::_newParas() {
  kappa[0] = (*O815->paraQ)["kappaone"];
  kappa[1] = (*O815->paraQ)["kappatwo"];
  lambda[0] = (*O815->paraQ)["lambdaone"];
  lambda[1] = (*O815->paraQ)["lambdatwo"];
  beta = (*O815->paraQ)["beta"];
}

void sim::_resetConfig() {
  for(int ix=0; ix<lsize4; ix++) {
    for(int i=0; i<2; i++) phi[ i*lsize4 + ix ] = 0;
    for(int nu=0; nu<4; nu++) U[ ix*4 + nu ] = 1;
  }
}

int sim::updateU(const int& x0, const int& nu0)
{
  complex<double> candU = U[x0*4+nu0] * polar(1.0, 2*EPSILONU*( 0.5 - gsl_rng_uniform(rangsl) ));

  if ( gsl_rng_uniform(rangsl) <= rhoU(x0, nu0, candU) ) {
    U[x0*4 + nu0] = candU;
    return 1;
  }

  return 0;
}

int sim::updatePhi(const int& iphi, const int& x0)
{
  complex<double> candPhi = phi[ iphi*lsize4 + x0 ] + 
    complex<double>(2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ), 
		    2*EPSILONPHI*( 0.5 - gsl_rng_uniform(rangsl) ));

  if ( gsl_rng_uniform(rangsl) <= rhoPhi(iphi, x0, candPhi) ) {
    phi[ iphi*lsize4 + x0 ] = candPhi;
    return 1;
  }

  return 0;
}

double sim::rhoPhi(const int& iphi, const int& x0, const complex<double>& candPhi)
{
  double deltaS=0;

  for( int mu=0; mu<4; mu++) {
    if( iphi == 0 ) {
      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
      deltaS -= 2 * real( conj(candPhi) * U[ x0*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
      deltaS -= 2 * real( conj(candPhi) * conj(U[ (*nb)[x0*8+mu+4]*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
    }
    else if( iphi == 1 ) {
      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
      deltaS += 2 * real( conj(phi[ iphi*lsize4 + x0 ]) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
      deltaS -= 2 * real( conj(candPhi) * conj(U[ x0*4 + mu ]) * phi[ iphi*lsize4 + (*nb)[x0*8+mu] ] );
      deltaS -= 2 * real( conj(candPhi) * U[ (*nb)[x0*8+mu+4]*4 + mu ] * phi[ iphi*lsize4 + (*nb)[x0*8+mu+4] ] );
    }
  }

  deltaS -= kappa[iphi] * norm(phi[ iphi*lsize4 + x0 ]);
  deltaS += kappa[iphi] * norm(candPhi);

  deltaS -= lambda[iphi] * pow(norm(phi[ iphi*lsize4 + x0 ]),2);
  deltaS += lambda[iphi] * pow(norm(candPhi),2);

  return exp(-deltaS);
}

double sim::rhoU(const int& x0, const int& nu0, const complex<double>& candU)
{
  double deltaS=0;

  for( int nu=0; nu<4; nu++ ) {
    if( nu == nu0 ) continue;
    deltaS += beta * real( U[x0*4+nu0] * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
    deltaS += beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(U[ x0*4 + nu0 ]) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
    deltaS -= beta * real( candU * U[ (*nb)[x0*8+nu0]*4 + nu ] * conj(U[ (*nb)[x0*8+nu]*4 + nu0 ]) * conj(U[ x0*4 + nu ]) );
    deltaS -= beta * real( U[ (*nb)[x0*8+nu+4]*4 + nu0 ] * U[ (*nb)[ (*nb)[x0*8+nu+4]*8+nu0 ]*4 + nu ] * conj(candU) * conj(U[ (*nb)[x0*8+nu+4]*4 + nu ]) );
  }

  deltaS += 2 * real( conj(phi[ 0*lsize4 + x0 ]) * U[ x0*4 + nu0 ] * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
  deltaS -= 2 * real( conj(phi[ 0*lsize4 + x0 ]) * candU * phi[ 0*lsize4 + (*nb)[x0*8+nu0] ]  );
  
  deltaS += 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(U[ x0*4 + nu0 ]) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );
  deltaS -= 2 * real( conj(phi[ 1*lsize4 + x0 ]) * conj(candU) * phi[ 1*lsize4 + (*nb)[x0*8+nu0] ]  );

  return exp(-deltaS);
}

#endif