+#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 = ¶
+
+ 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
projects / phys / u1casc.git / commitdiff