u1casc-flux/sim.hpp

   1 #ifndef SIM_HPP
   2 #define SIM_HPP
   3
   4 #define BESSELCACHE 10000
   5 #define WCACHE 1000
   6
   7 #include <gsl/gsl_rng.h>
   8 #include <gsl/gsl_errno.h>
   9 #include <gsl/gsl_integration.h>
  10 #include <math.h>
  11 #include <sys/time.h>
  12 #include <boost/math/special_functions/bessel.hpp>
  13
  14 #include "latlib/neigh.h"
  15
  16 class sim : public o815::sim {
  17 public:
  18   sim(o815 *_O815);
  19   unsigned int lsize4;
  20   neigh *nb;
  21   int *l, *lp, *p;
  22   double kappa[2], lambda[2], Mu[2], beta;
  23   int wArg(int x, int i);
  24   double WF(int nom, int denom, int iflav);
  25   int iPlaq(int i, int j);
  26   double I(int p);
  27   void _newParas();
  28
  29 private:
  30   void _resetConfig();
  31   void _makeSweep();
  32   gsl_rng* rangsl;
  33   int lineSweep();
  34   int lineUpdate(int x, int vardir);
  35   int cubeUpdate(int x, int orient);
  36   void getCube(int cube[3], int orient);
  37   int plaqUpdate(int x, int mu, int nu, int i);
  38   int lpUpdate(int x, int mu, int i);
  39   static double wkern(double t, void *params);
  40   double wGsl(int n, int iflav);
  41   struct bcache {
  42     double val;
  43     bool set;
  44   };
  45   bcache besselCache[BESSELCACHE];
  46   struct kernpara {
  47     int n;
  48     double kappa;
  49     double lambda;
  50   };
  51   struct wcache {
  52     double val;
  53     double error;
  54     bool set;
  55   };
  56   gsl_integration_cquad_workspace *w;
  57   double rhoLine(int xstart, int vardir, int modifier);
  58   double rhoPlaq_Fac(int x, int mu, int i, int modifier);
  59   double rhoPlaq(int x, int mu, int nu, int i, int modifier);
  60   double rhoLp(int x, int mu, int i, int modifier);
  61   bool neverComputedB, neverComputedW[2];
  62   double rhoCube(int x, int *cube, int modifier);
  63   double wLastLambda[2], wLastKappa[2];
  64   double lastBesselBeta;
  65   int ipow(int base, int expo);
  66   wcache wCache[2][WCACHE*WCACHE];
  67   double WFl0(int nom, int denom, int iflav);
  68 };
  69
  70 sim::sim(o815 *_O815) : o815::sim( _O815,
  71                                    sizeof(int)*
  72                                    _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2*4+2*4+6) ) {
  73
  74   struct timeval tv;
  75
  76   lsize4 = _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1];
  77
  78   nb = new neigh(4, _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
  79
  80   l = (int*)confMem;
  81   lp = (int*)(confMem + sizeof(int)*2*lsize4*4);
  82   p = (int*)(confMem + sizeof(int)*2*lsize4*4*2);
  83
  84   neverComputedB = true;
  85   neverComputedW[0] = true;
  86   neverComputedW[1] = true;
  87
  88   w = gsl_integration_cquad_workspace_alloc(10000);
  89
  90   gettimeofday(&tv, NULL);
  91   rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
  92   gsl_rng_set(rangsl, 1000000 * tv.tv_sec + tv.tv_usec);
  93 }
  94
  95 void sim::_makeSweep()
  96 {
  97   lineSweep();
  98
  99   for( int x=0; x<lsize4; x++) {
 100     for(int orient=0; orient<4; orient++)
 101       cubeUpdate(x, orient);
 102
 103     for(int mu=0; mu<3; mu++)
 104       for(int nu=mu+1; nu<4; nu++)
 105         for(int i=0; i<2; i++)
 106           plaqUpdate(x, mu, nu, i);
 107
 108     for(int mu=0; mu<4; mu++)
 109       for(int i=0; i<2; i++)
 110         lpUpdate(x, mu, i);
 111   }
 112 }
 113
 114 void sim::_newParas() {
 115   kappa[0] = (*O815->paraQ)["kappaone"];
 116   kappa[1] = (*O815->paraQ)["kappatwo"];
 117   lambda[0] = (*O815->paraQ)["lambdaone"];
 118   lambda[1] = (*O815->paraQ)["lambdatwo"];
 119   beta = (*O815->paraQ)["beta"];
 120   Mu[0] = (*O815->paraQ)["muone"];
 121   Mu[1] = (*O815->paraQ)["mutwo"];
 122 }
 123
 124 void sim::_resetConfig() {
 125   memset(confMem, 0, confmemSize);
 126 }
 127
 128 int sim::lineSweep()
 129 {
 130   for (int vardir = 0; vardir < 4; vardir++)
 131     for (int x0 = 0; x0 < O815->comargs.lsize[0]; x0++) {
 132       for (int x1 = 0; x1 < O815->comargs.lsize[0]; x1++) {
 133         for (int x2 = 0; x2 < O815->comargs.lsize[0]; x2++) {
 134           for (int x3 = 0; x3 < O815->comargs.lsize[1]; x3++) {
 135             lineUpdate(x0
 136                        + x1 * O815->comargs.lsize[0]
 137                        + x2 * O815->comargs.lsize[0] * O815->comargs.lsize[0]
 138                        + x3 * O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0], vardir);
 139             if(vardir==3) break;
 140           }
 141           if(vardir==2) break;
 142         }
 143         if(vardir==1) break;
 144       }
 145       if(vardir==0) break;
 146     }
 147 }
 148
 149 int sim::lineUpdate(int x, int vardir)
 150 {
 151   int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
 152
 153   if ( gsl_rng_uniform(rangsl) < rhoLine(x, vardir, modifier) )
 154     {
 155       const int stepsize = ipow(O815->comargs.lsize[0], vardir);
 156
 157       for (int xvar = 0; xvar < O815->comargs.lsize[vardir/3]; xvar++)
 158         for (int i = 0; i < 2; i++)
 159           l[ i*lsize4*4 + (x+xvar*stepsize)*4 + vardir ] += modifier;
 160
 161       return 1;
 162     }
 163
 164   return 0;
 165 }
 166
 167 int sim::cubeUpdate(int x, int orient)
 168 {
 169   int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
 170   int cube[3];
 171   getCube(cube, orient);
 172
 173   if ( gsl_rng_uniform(rangsl) < rhoCube(x, cube, modifier) )
 174     {
 175       p[ x*6 + iPlaq( cube[0],cube[1] ) ] += modifier;
 176       p[ x*6 + iPlaq( cube[0],cube[2] ) ] -= modifier;
 177       p[ x*6 + iPlaq( cube[1],cube[2] ) ] += modifier;
 178       p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] )  ] -= modifier;
 179       p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] )  ] += modifier;
 180       p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] )  ] -= modifier;
 181
 182       return 1;
 183     }
 184
 185   return 0;
 186 }
 187
 188 void sim::getCube(int cube[3], int orient)
 189 {
 190   int skipped=0;
 191   for (int i=0; i<4; i++) {
 192     if ( i == orient ) {
 193       skipped++;
 194       continue;
 195     }
 196     cube[i-skipped] = i;
 197   }
 198 }
 199
 200 int sim::plaqUpdate(int x, int mu, int nu, int i)
 201 {
 202   int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
 203   int lmodsign = 1 - 2*i;
 204
 205   if ( gsl_rng_uniform(rangsl) < rhoPlaq(x, mu, nu, i, modifier) )
 206     {
 207       p[ x*6 + iPlaq(mu,nu) ] = p[ x*6 + iPlaq(mu,nu) ] + modifier;
 208       l[ i*lsize4*4 + x*4 + mu ] = l[ i*lsize4*4 + x*4 + mu ] - lmodsign*modifier;
 209       l[ i*lsize4*4 + x*4 + nu ] = l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier;
 210       l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] = l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier;
 211       l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] = l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier;
 212       return 1;
 213     }
 214
 215   return 0;
 216 }
 217
 218 int sim::iPlaq(int i, int j)
 219 {
 220   if ( i == 0 ) return j-1;
 221   else if ( i == 1 ) return j+1;
 222   else return j+2;
 223 }
 224
 225 int sim::lpUpdate(int x, int mu, int i)
 226 {
 227   int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
 228
 229   if ( (lp[ i*lsize4*4 + x*4 + mu ] + modifier) < 0 )
 230     return 0;
 231
 232   if ( gsl_rng_uniform(rangsl) < rhoLp(x, mu, i, modifier) ) {
 233     lp[ i*lsize4*4 + x*4 + mu ] += modifier;
 234     return 1;
 235   }
 236
 237   return 0;
 238 }
 239
 240 double sim::I(int p)
 241 {
 242   /* http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html */
 243
 244   p = abs(p);
 245
 246   if ( beta != lastBesselBeta || neverComputedB )
 247     {
 248       memset(besselCache, 0, sizeof(bcache)*BESSELCACHE);
 249       lastBesselBeta = beta;
 250       neverComputedB = false;
 251     }
 252
 253   if ( p < BESSELCACHE )
 254     {
 255       if ( besselCache[p].set ) return besselCache[p].val;
 256       else
 257         {
 258           besselCache[p].val = boost::math::cyl_bessel_i(p, beta);
 259           besselCache[p].set = true;
 260           return besselCache[p].val;
 261         }
 262     }
 263   else
 264     {
 265       return boost::math::cyl_bessel_i(p, beta);
 266     }
 267 }
 268
 269 int sim::wArg(int x, int i)
 270 {
 271   int warg = 0;
 272
 273   for (int nu=0; nu<4; nu++)
 274     warg += abs( l[ i*lsize4*4 + x*4 + nu ] )
 275       + 2 * ( lp[ i*lsize4*4 + x*4 + nu ] + lp[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] )
 276       + abs( l[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] );
 277
 278   return warg;
 279 }
 280
 281 double sim::WF(int nom, int denom, int iflav)
 282 {
 283   if ( nom == denom ) return 1;
 284
 285   if ( wLastLambda[iflav] != lambda[iflav] || wLastKappa[iflav] != kappa[iflav] || neverComputedW[iflav] )
 286     {
 287       memset(wCache[iflav], 0, sizeof(wcache)*WCACHE*WCACHE);
 288       wLastLambda[iflav] = lambda[iflav];
 289       wLastKappa[iflav] = kappa[iflav];
 290       neverComputedW[iflav] = false;
 291     }
 292
 293   if ( nom/2 < WCACHE && denom/2 < WCACHE )
 294     {
 295       if ( wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set )
 296         return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
 297       else
 298         {
 299           /* integral can be solved analytically for lambda==0 */
 300           if( lambda == 0 )
 301             wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = WFl0(nom, denom, iflav);
 302           else
 303             wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = wGsl(nom, iflav) / wGsl(denom, iflav);
 304
 305           wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set = true;
 306
 307           return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
 308         }
 309     }
 310   else
 311     {
 312       if( lambda == 0 )
 313         return WFl0(nom, denom, iflav);
 314       else
 315         return wGsl(nom, iflav) / wGsl(denom, iflav);
 316     }
 317 }
 318
 319 double sim::wkern(double t, void *params)
 320 {
 321   kernpara *para = (kernpara*)params;
 322   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);
 323   return wkern;
 324 }
 325
 326 double sim::wGsl(int n, int iflav)
 327 {
 328   double result;
 329   gsl_function gslF;
 330   kernpara para;
 331   int status;
 332
 333   para.n = n;
 334   para.kappa = kappa[iflav];
 335   para.lambda = lambda[iflav];
 336
 337   gslF.function = &wkern;
 338   gslF.params = &para;
 339
 340   status = gsl_integration_cquad(&gslF, 0, 1, 0, 1e-7, w, &result, NULL, NULL);
 341
 342   return result;
 343 }
 344
 345 double sim::rhoLine(int xstart, int vardir, int modifier)
 346 {
 347   double rho=1;
 348   int tmpwarg;
 349
 350   const int stepsize = ipow(O815->comargs.lsize[0],vardir);
 351   for(int xvar=0; xvar<O815->comargs.lsize[vardir/3]; xvar++)
 352     for(int i=0; i<2; i++)
 353       {
 354         const int x = xstart + xvar*stepsize;
 355
 356         if( abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + vardir ] ) )
 357           rho /= abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) + lp[ i*lsize4*4 + x*4 + vardir ];
 358         else
 359           rho *= abs( l[ i*lsize4*4 + x*4 + vardir ] ) + lp[ i*lsize4*4 + x*4 + vardir ];
 360
 361         tmpwarg = wArg(x,i);
 362         rho *= WF( tmpwarg
 363                    - abs( l[ i*lsize4*4 + x*4 + vardir ] )
 364                    + abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier )
 365                    - abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] )
 366                    + abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] + modifier ),
 367                    tmpwarg, i );
 368
 369         if(vardir==3)
 370           rho *= exp(modifier*Mu[i]);
 371       }
 372
 373   return rho;
 374 }
 375
 376 double sim::rhoPlaq(int x, int mu, int nu, int i, int modifier)
 377 {
 378   double rho = 1;
 379   int lmodsign = 1 - 2*i;
 380   int tmpwarg;
 381
 382   rho *= I( p[ x*6 + iPlaq(mu,nu) ] + modifier ) / I( p[ x*6 + iPlaq(mu,nu) ] );
 383
 384   tmpwarg = wArg(x,i);
 385   rho *= WF( tmpwarg
 386              - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
 387              + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
 388              tmpwarg, i );
 389
 390   tmpwarg = wArg((*nb)[x*8+mu],i);
 391   rho *= WF( tmpwarg
 392              - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] )
 393              + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier ),
 394              tmpwarg, i );
 395
 396   tmpwarg = wArg((*nb)[(*nb)[x*8+mu]*8+nu],i);
 397   rho *= WF( tmpwarg
 398              - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] )
 399              + 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 ),
 400              tmpwarg, i );
 401
 402   tmpwarg = wArg((*nb)[x*8+nu],i);
 403   rho *= WF( tmpwarg
 404              - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
 405              + abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
 406              tmpwarg, i );
 407
 408   rho *= rhoPlaq_Fac(x, mu, i, -lmodsign*modifier);
 409   rho *= rhoPlaq_Fac(x, nu, i, +lmodsign*modifier);
 410   rho *= rhoPlaq_Fac((*nb)[x*8+mu], nu, i, -lmodsign*modifier);
 411   rho *= rhoPlaq_Fac((*nb)[x*8+nu], mu, i, +lmodsign*modifier);
 412
 413   return rho;
 414 }
 415
 416 double sim::rhoPlaq_Fac(int x, int mu, int i, int modifier)
 417 {
 418   if( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + mu] ) )
 419     return 1.0 / ( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) + lp[ i*lsize4*4 + x*4 + mu ] );
 420   else
 421     return ( abs( l[ i*lsize4*4 + x*4 + mu ] ) + lp[ i*lsize4*4 + x*4 + mu ] );
 422 }
 423
 424 double sim::rhoLp(int x, int mu, int i, int modifier)
 425 {
 426   double rho=1;
 427   int tmpwarg;
 428
 429   if(modifier > 0)
 430       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 );
 431   else
 432       rho *= ( abs(l[ i*lsize4*4 + x*4 + mu ]) + lp[ i*lsize4*4 + x*4 + mu ] ) * lp[ i*lsize4*4 + x*4 + mu ];
 433
 434   tmpwarg = wArg(x,i);
 435   rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
 436   tmpwarg = wArg((*nb)[x*8+mu],i);
 437   rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
 438
 439   return rho;
 440 }
 441
 442 double sim::rhoCube(int x, int *cube, int modifier)
 443 {
 444   double rho=1;
 445   rho *= I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] );
 446   rho *= I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] - modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] );
 447   rho *= I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] );
 448   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] )  ] );
 449   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] )  ] );
 450   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] )  ] );
 451   return rho;
 452 }
 453
 454 int sim::ipow(int base, int expo)
 455 {
 456   int ipow = 1;
 457   while( expo>0 )
 458     {
 459       ipow *= base;
 460       expo--;
 461     }
 462   return ipow;
 463 }
 464
 465 double sim::WFl0(int nom, int denom, int iflav)
 466 {
 467   double wfl0 = 1;
 468
 469   if (nom > denom) {
 470     for (int i=nom/2; i>denom/2; i--)
 471       wfl0 *= i;
 472     return wfl0 * pow(kappa[iflav], 0.5*(denom-nom));
 473   }
 474   else {
 475     for (int i=denom/2; i>nom/2; i--)
 476       wfl0 *= i;
 477     return (1.0 / wfl0) * pow(kappa[iflav], 0.5*(denom-nom));
 478   }
 479 }
 480
 481 #endif