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