X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/blobdiff_plain/2a5f627da5468eb738f1b35aba77652729c1d7a1..7a5913c428ac75aa2de2330fa926e8deb3457112:/u1casc-flux/sim.hpp

diff --git a/u1casc-flux/sim.hpp b/u1casc-flux/sim.hpp
new file mode 100644
index 0000000..9eca3d7
--- /dev/null
+++ b/u1casc-flux/sim.hpp
@@ -0,0 +1,478 @@
+#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;
+
+private:
+  void _makeSweep();
+  void _newParas();
+  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 iPlaq(int i, int j);
+  int lpUpdate(int x, int mu, int i);
+  double I(int p);
+  int wArg(int x, int i);
+  double WF(int nom, int denom, int iflav);
+  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"];
+
+  memset(confMem, 0, sizeof(int)*lsize4*(2*4+2*4+6));
+}
+
+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