#ifndef SIM_HPP
#define SIM_HPP

#define BESSELCACHE 10000
#define WCACHE 1000

#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <math.h>
#include <sys/time.h>
#include <boost/math/special_functions/bessel.hpp>

#include "latlib/neigh.h"

class sim : public o815::sim {
public:
  sim(o815 *_O815);
  unsigned int lsize4;
  neigh *nb;
  int *l, *lp, *p;
  double kappa[2], lambda[2], Mu[2], beta;
  int wArg(int x, int i);
  double WF(int nom, int denom, int iflav);
  int iPlaq(int i, int j);
  double I(int p);
  void _newParas();
  
private:
  void _resetConfig();
  void _makeSweep();
  gsl_rng* rangsl;
  int lineSweep();
  int lineUpdate(int x, int vardir);
  int cubeUpdate(int x, int orient);
  void getCube(int cube[3], int orient);
  int plaqUpdate(int x, int mu, int nu, int i);
  int lpUpdate(int x, int mu, int i);
  static double wkern(double t, void *params);
  double wGsl(int n, int iflav);
  struct bcache { 
    double val; 
    bool set; 
  };
  bcache besselCache[BESSELCACHE]; 
  struct kernpara { 
    int n; 
    double kappa; 
    double lambda; 
  };
  struct wcache { 
    double val; 
    double error; 
    bool set; 
  };
  gsl_integration_cquad_workspace *w;
  double rhoLine(int xstart, int vardir, int modifier);
  double rhoPlaq_Fac(int x, int mu, int i, int modifier);
  double rhoPlaq(int x, int mu, int nu, int i, int modifier);
  double rhoLp(int x, int mu, int i, int modifier);
  bool neverComputedB, neverComputedW[2];
  double rhoCube(int x, int *cube, int modifier);
  double wLastLambda[2], wLastKappa[2];
  double lastBesselBeta;
  int ipow(int base, int expo);
  wcache wCache[2][WCACHE*WCACHE];
  double WFl0(int nom, int denom, int iflav);
};

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

  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]);

  l = (int*)confMem;
  lp = (int*)(confMem + sizeof(int)*2*lsize4*4);
  p = (int*)(confMem + sizeof(int)*2*lsize4*4*2);

  neverComputedB = true;
  neverComputedW[0] = true;
  neverComputedW[1] = true;

  w = gsl_integration_cquad_workspace_alloc(10000);

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

void sim::_makeSweep() 
{ 
  lineSweep();
  
  for( int x=0; x<lsize4; x++) {
    for(int orient=0; orient<4; orient++) 
      cubeUpdate(x, orient);
    
    for(int mu=0; mu<3; mu++) 
      for(int nu=mu+1; nu<4; nu++) 
	for(int i=0; i<2; i++) 
	  plaqUpdate(x, mu, nu, i);
    
    for(int mu=0; mu<4; mu++) 
      for(int i=0; i<2; i++) 
	lpUpdate(x, mu, i);
  }
}

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"];
  Mu[0] = (*O815->paraQ)["muone"];
  Mu[1] = (*O815->paraQ)["mutwo"];
}

void sim::_resetConfig() {
  memset(confMem, 0, confmemSize);
}

int sim::lineSweep()
{
  for (int vardir = 0; vardir < 4; vardir++)
    for (int x0 = 0; x0 < O815->comargs.lsize[0]; x0++) {
      for (int x1 = 0; x1 < O815->comargs.lsize[0]; x1++) {
	for (int x2 = 0; x2 < O815->comargs.lsize[0]; x2++) {
	  for (int x3 = 0; x3 < O815->comargs.lsize[1]; x3++) {
	    lineUpdate(x0 
		       + x1 * O815->comargs.lsize[0] 
		       + x2 * O815->comargs.lsize[0] * O815->comargs.lsize[0] 
		       + x3 * O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0], vardir);  
	    if(vardir==3) break;
	  }
	  if(vardir==2) break;
	}
	if(vardir==1) break;
      }
      if(vardir==0) break;
    }
}

int sim::lineUpdate(int x, int vardir)
{
  int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);

  if ( gsl_rng_uniform(rangsl) < rhoLine(x, vardir, modifier) )
    {
      const int stepsize = ipow(O815->comargs.lsize[0], vardir);

      for (int xvar = 0; xvar < O815->comargs.lsize[vardir/3]; xvar++)
	for (int i = 0; i < 2; i++)
	  l[ i*lsize4*4 + (x+xvar*stepsize)*4 + vardir ] += modifier;
      
      return 1;
    }
  
  return 0;
}

int sim::cubeUpdate(int x, int orient)
{
  int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
  int cube[3];
  getCube(cube, orient);

  if ( gsl_rng_uniform(rangsl) < rhoCube(x, cube, modifier) )
    {
      p[ x*6 + iPlaq( cube[0],cube[1] ) ] += modifier;
      p[ x*6 + iPlaq( cube[0],cube[2] ) ] -= modifier;
      p[ x*6 + iPlaq( cube[1],cube[2] ) ] += modifier;
      p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] )  ] -= modifier;
      p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] )  ] += modifier;
      p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] )  ] -= modifier;
      
      return 1;
    }
  
  return 0;
}

void sim::getCube(int cube[3], int orient)
{
  int skipped=0;
  for (int i=0; i<4; i++) {
    if ( i == orient ) {
      skipped++;
      continue;
    }
    cube[i-skipped] = i;
  }
}

int sim::plaqUpdate(int x, int mu, int nu, int i)
{
  int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
  int lmodsign = 1 - 2*i;

  if ( gsl_rng_uniform(rangsl) < rhoPlaq(x, mu, nu, i, modifier) )
    {
      p[ x*6 + iPlaq(mu,nu) ] = p[ x*6 + iPlaq(mu,nu) ] + modifier;
      l[ i*lsize4*4 + x*4 + mu ] = l[ i*lsize4*4 + x*4 + mu ] - lmodsign*modifier;
      l[ i*lsize4*4 + x*4 + nu ] = l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier;
      l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] = l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier;
      l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] = l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier;
      return 1;
    }
  
  return 0;
}

int sim::iPlaq(int i, int j)
{
  if ( i == 0 ) return j-1;
  else if ( i == 1 ) return j+1;
  else return j+2;
}

int sim::lpUpdate(int x, int mu, int i)
{
  int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);

  if ( (lp[ i*lsize4*4 + x*4 + mu ] + modifier) < 0 ) 
    return 0;

  if ( gsl_rng_uniform(rangsl) < rhoLp(x, mu, i, modifier) ) {
    lp[ i*lsize4*4 + x*4 + mu ] += modifier;
    return 1;
  }

  return 0;
}

double sim::I(int p)
{
  /* http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html */

  p = abs(p);

  if ( beta != lastBesselBeta || neverComputedB ) 
    {
      memset(besselCache, 0, sizeof(bcache)*BESSELCACHE);
      lastBesselBeta = beta;
      neverComputedB = false;
    }

  if ( p < BESSELCACHE )
    {
      if ( besselCache[p].set ) return besselCache[p].val;
      else
	{
	  besselCache[p].val = boost::math::cyl_bessel_i(p, beta);
	  besselCache[p].set = true;
	  return besselCache[p].val;
	}
    }
  else
    {
      return boost::math::cyl_bessel_i(p, beta);
    }
}

int sim::wArg(int x, int i)
{
  int warg = 0;

  for (int nu=0; nu<4; nu++)
    warg += abs( l[ i*lsize4*4 + x*4 + nu ] ) 
      + 2 * ( lp[ i*lsize4*4 + x*4 + nu ] + lp[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] )
      + abs( l[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] ); 

  return warg;
}

double sim::WF(int nom, int denom, int iflav)
{
  if ( nom == denom ) return 1;
  
  if ( wLastLambda[iflav] != lambda[iflav] || wLastKappa[iflav] != kappa[iflav] || neverComputedW[iflav] )
    {
      memset(wCache[iflav], 0, sizeof(wcache)*WCACHE*WCACHE);
      wLastLambda[iflav] = lambda[iflav]; 
      wLastKappa[iflav] = kappa[iflav];
      neverComputedW[iflav] = false;
    }
  
  if ( nom/2 < WCACHE && denom/2 < WCACHE )
    {
      if ( wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set )
	return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
      else
	{
	  /* integral can be solved analytically for lambda==0 */
	  if( lambda == 0 )
	    wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = WFl0(nom, denom, iflav);
	  else
	    wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = wGsl(nom, iflav) / wGsl(denom, iflav);

	  wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set = true;

	  return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
	}
    }
  else
    {
      if( lambda == 0 )
	return WFl0(nom, denom, iflav);
      else
	return wGsl(nom, iflav) / wGsl(denom, iflav);
    }
}

double sim::wkern(double t, void *params)
{
  kernpara *para = (kernpara*)params;
  double wkern = pow((1-t)/t, para->n + 1) * exp( -para->kappa*pow((1-t)/t, 2) -para->lambda*pow((1-t)/t, 4) ) / pow(t,2);
  return wkern;
}

double sim::wGsl(int n, int iflav)
{
  double result;
  gsl_function gslF;
  kernpara para;
  int status;

  para.n = n;
  para.kappa = kappa[iflav];
  para.lambda = lambda[iflav];

  gslF.function = &wkern;
  gslF.params = &para;

  status = gsl_integration_cquad(&gslF, 0, 1, 0, 1e-7, w, &result, NULL, NULL);

  return result;
}

double sim::rhoLine(int xstart, int vardir, int modifier)
{
  double rho=1;
  int tmpwarg;

  const int stepsize = ipow(O815->comargs.lsize[0],vardir);
  for(int xvar=0; xvar<O815->comargs.lsize[vardir/3]; xvar++)
    for(int i=0; i<2; i++)
      {
	const int x = xstart + xvar*stepsize;
	
	if( abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + vardir ] ) ) 
	  rho /= abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) + lp[ i*lsize4*4 + x*4 + vardir ];
	else
	  rho *= abs( l[ i*lsize4*4 + x*4 + vardir ] ) + lp[ i*lsize4*4 + x*4 + vardir ];

	tmpwarg = wArg(x,i);
	rho *= WF( tmpwarg 
		   - abs( l[ i*lsize4*4 + x*4 + vardir ] ) 
		   + abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier )
		   - abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] ) 
		   + abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] + modifier ),
		   tmpwarg, i );

	if(vardir==3)
	  rho *= exp(modifier*Mu[i]);
      }
  
  return rho;
}

double sim::rhoPlaq(int x, int mu, int nu, int i, int modifier)
{
  double rho = 1;
  int lmodsign = 1 - 2*i;
  int tmpwarg;

  rho *= I( p[ x*6 + iPlaq(mu,nu) ] + modifier ) / I( p[ x*6 + iPlaq(mu,nu) ] );
  
  tmpwarg = wArg(x,i);
  rho *= WF( tmpwarg 
	     - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
	     + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
	     tmpwarg, i );
  
  tmpwarg = wArg((*nb)[x*8+mu],i);
  rho *= WF( tmpwarg 
	     - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] )
	     + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier ),
	     tmpwarg, i );

  tmpwarg = wArg((*nb)[(*nb)[x*8+mu]*8+nu],i);
  rho *= WF( tmpwarg
	     - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] ) 
	     + abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier ),
	     tmpwarg, i );
  
  tmpwarg = wArg((*nb)[x*8+nu],i);
  rho *= WF( tmpwarg
	     - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
	     + abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
	     tmpwarg, i );
  
  rho *= rhoPlaq_Fac(x, mu, i, -lmodsign*modifier);
  rho *= rhoPlaq_Fac(x, nu, i, +lmodsign*modifier);
  rho *= rhoPlaq_Fac((*nb)[x*8+mu], nu, i, -lmodsign*modifier);
  rho *= rhoPlaq_Fac((*nb)[x*8+nu], mu, i, +lmodsign*modifier);

  return rho;
}

double sim::rhoPlaq_Fac(int x, int mu, int i, int modifier)
{
  if( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + mu] ) )
    return 1.0 / ( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) + lp[ i*lsize4*4 + x*4 + mu ] );
  else
    return ( abs( l[ i*lsize4*4 + x*4 + mu ] ) + lp[ i*lsize4*4 + x*4 + mu ] );
}

double sim::rhoLp(int x, int mu, int i, int modifier)
{
  double rho=1;
  int tmpwarg;

  if(modifier > 0) 
      rho /= ( abs(l[ i*lsize4*4 + x*4 + mu ]) + lp[ i*lsize4*4 + x*4 + mu ] + modifier ) * ( lp[ i*lsize4*4 + x*4 + mu ] + modifier );
  else
      rho *= ( abs(l[ i*lsize4*4 + x*4 + mu ]) + lp[ i*lsize4*4 + x*4 + mu ] ) * lp[ i*lsize4*4 + x*4 + mu ];
  
  tmpwarg = wArg(x,i);
  rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
  tmpwarg = wArg((*nb)[x*8+mu],i);
  rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
  
  return rho;
}

double sim::rhoCube(int x, int *cube, int modifier)
{
  double rho=1;
  rho *= I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] );
  rho *= I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] - modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] );
  rho *= I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] );
  rho *= I( p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] )  ] - modifier ) / I( p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] )  ] );
  rho *= I( p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] )  ] + modifier ) / I( p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] )  ] );
  rho *= I( p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] )  ] - modifier ) / I( p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] )  ] );
  return rho;
}

int sim::ipow(int base, int expo)
{
  int ipow = 1;
  while( expo>0 ) 
    {
      ipow *= base;
      expo--;
    }
  return ipow;
}

double sim::WFl0(int nom, int denom, int iflav)
{
  double wfl0 = 1;

  if (nom > denom) {
    for (int i=nom/2; i>denom/2; i--)
      wfl0 *= i;
    return wfl0 * pow(kappa[iflav], 0.5*(denom-nom));
  }
  else {
    for (int i=denom/2; i>nom/2; i--)
      wfl0 *= i;
    return (1.0 / wfl0) * pow(kappa[iflav], 0.5*(denom-nom));
  }
}

#endif