{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Monte Carlo integration by Vegas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Brute force Monte Carlo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we will implement the simple Monte Carlo method." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# simple class to take care of statistically determined quantities\n", "class Cumulants:\n", " def __init__(self):\n", " self.sum=0.0 # f_0 + f_1 +.... + f_N\n", " self.sqsum=0.0 # f_0^2 + f_1^2 +....+ f_N^2\n", " self.avg = 0.0 # I_best when many iterations, otherwise = 1/N\\sum_i f_i\n", " self.err = 0.0 # sigma of I_best when many iterations, otherwise sqrt( -^2 )/sqrt(N)\n", " self.chisq = 0.0\n", " self.weightsum=0.0 # \\sum_i 1/sigma_i^2\n", " self.avgsum=0.0 # \\sum_i _i/sigma_i^2\n", " self.avg2sum=0.0 # \\sum_i _i^2/sigma_i^2\n", "\n", "def SimpleMC(integrant, ndim, unit, maxeval, cum):\n", " nbatch=1000 # function will be evaluated in bacthes of 1000 evaluations at one time (for efficiency and storage issues)\n", " neval=0\n", " for nsamples in range(maxeval,0,-nbatch): # loop over all_nsample evaluations in batches of nbatch\n", " n = min(nbatch,nsamples) # How many evaluations in this pass?\n", " xr = unit*random.random((n,ndim)) # generates 2-d array of random numbers in the interval [0,1)\n", " wfun = integrant(xr) # n function evaluations required in single call\n", " neval += n # We just added so many fuction evaluations\n", " cum.sum += sum(wfun) # sum_i f_i = * neval\n", " cum.sqsum += sum(wfun*wfun) # sum_i f_i^2 = * neval\n", " \n", " cum.avg = cum.sum/neval\n", " w0 = sqrt(cum.sqsum/neval) # sqrt()\n", " cum.err = sqrt((w0-cum.avg)*(w0+cum.avg)/neval) # sqrt(sqsum-sum**2)\n", " cum.avg *= unit**ndim # adding units if not [0,1] interval\n", " cum.err *= unit**ndim # adding units if not [0,1] interval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We used here a simple trick to avoid overflow, i.e.,\n", "\n", "\\begin{eqnarray}\n", "\\sqrt{\\frac{\\langle f^2\\rangle-\\langle f\\rangle^2}{N}} = \n", "\\sqrt{\\frac{(\\sqrt{\\langle f^2\\rangle}-\\langle f\\rangle)(\\sqrt{\\langle f^2\\rangle}+\\langle f\\rangle)}{N}}\n", "\\end{eqnarray}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_integrant2(x):\n", " \"\"\" For testing, we are integration the function\n", " 1/(1-cos(x)*cos(y)*cos(z))/pi^3\n", " in the interval [0,pi]**3\n", " \"\"\"\n", " #nbatch,ndim = shape(x)\n", " return 1.0/(1.0-cos(x[:,0])*cos(x[:,1])*cos(x[:,2]))/pi**3" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "unit=pi\n", "ndim=3\n", "maxeval=200000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "SimpleMC(my_integrant2, ndim, pi, maxeval, cum)\n", "print(cum.avg, '+-', cum.err, 'exact=', exact)\n", "print('how many sigma away', abs(cum.avg-exact)/cum.err)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vegas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we define the grid points $g(x)$. At the beginning, we just set $g(x)=x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class Grid:\n", " \"\"\"Contains the grid points g_n(x) with x=[0...1], and g=[0...1]\n", " for Vegas integration. There are n-dim g_n functions.\n", " Constraints : g(0)=0 and g(1)=1.\n", " \"\"\"\n", " def __init__(self, ndim, nbins):\n", " self.g = zeros((ndim,nbins+1)) \n", " # a bit dirty trick: We will later use also g[-1] in interpolation, which should be set to zero, hence\n", " # we allocate dimension nbins+1, rather than nbinx\n", " self.ndim=ndim\n", " self.nbins=nbins\n", " # At the beginning we set g(x)=x\n", " # The grid-points are x_0 = 1/N, x_1 = 2/N, ... x_{N-1}=1.0. Note x_N=x_{-1}=0 \n", " # Note that g(0)=0, and we skip this point on the mesh.\n", " for idim in range(ndim):\n", " self.g[idim,:nbins] = arange(1,nbins+1)/float(nbins)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Vegas_step1(integrant, unit, maxeval, nstart, nincrease, grid, cum):\n", " ndim, nbins = grid.ndim,grid.nbins # dimension of the integral, size of the grid for binning in each direction\n", " unit_dim = unit**ndim # converts from unit cube integration to generalized cube with unit length\n", " nbatch=1000 # function will be evaluated in bacthes of 1000 evaluations at one time (for efficiency and storage issues)\n", " neval=0\n", " print (\"\"\"Vegas parameters:\n", " ndim = \"\"\"+str(ndim)+\"\"\"\n", " unit = \"\"\"+str(unit)+\"\"\"\n", " maxeval = \"\"\"+str(maxeval)+\"\"\"\n", " nstart = \"\"\"+str(nstart)+\"\"\"\n", " nincrease = \"\"\"+str(nincrease)+\"\"\"\n", " nbins = \"\"\"+str(nbins)+\"\"\"\n", " nbaths = \"\"\"+str(nbatch)+\"\\n\")\n", "\n", " bins = zeros((nbatch,ndim),dtype=int) # in which sampled bin does this point fall?\n", " \n", " all_nsamples = nstart\n", " for nsamples in range(all_nsamples,0,-nbatch): # loop over all_nsample evaluations in batches of nbatch\n", " n = min(nbatch,nsamples) # How many evaluations in this pass?\n", " # We are integrating f(g_1(x),g_2(y),g_3(z))*dg_1/dx*dg_2/dy*dg_3/dz dx*dy*dz\n", " # This is represented as 1/all_nsamples \\sum_{x_i,y_i,z_i} f(g_1(x_i),g_2(y_i),g_3(z_i))*dg_1/dx*dg_2/dy*dg_3/dz\n", " # where dg_1/dx = diff*NBINS\n", " wgh = zeros(nbatch) # weights for each random point in the batch\n", " xr = random.random((n,ndim)) # generates 2-d array of random numbers in the interval [0,1)\n", " # This part takes quite a lot of time and it would be nice to rewrite with numba!\n", " for i in range(n):\n", " weight = 1.0/all_nsamples\n", " for dim in range(ndim):\n", " # We want to evaluate the function f at point g(x), i.e, f(g_1(x),g_2(y),...)\n", " # Here we transform the points x,y,z -> g_1(x), g_2(y), g_3(z)\n", " # We hence want to evaluate g(x) ~ g(x[i]), where x is the random number and g is the grid function\n", " # The discretized g(t) is defined on the grid :\n", " # t[-1]=0, t[0]=1/N, t[1]=2/N, t[2]=3/N ... t[N-1]=1.\n", " # We know that g(0)=0 and g(1)=1, so that g[-1]=0.0 and g[N-1]=1.0\n", " # To interpolate g at x, we first compute i=int(x*N) and then we use linear interpolation\n", " # g(x) = g[i-1] + (g[i]-g[i-1])*(x*N-i) ; if i>0\n", " # g(x) = 0 + (g[0]-0)*(x*N-0) ; if i=0\n", " #\n", " pos = xr[i,dim]*nbins # which grid would it fit ? (x*N)\n", " ipos = int(pos) # the grid position is ipos : int(x*N)==i \n", " diff = grid.g[dim,ipos] - grid.g[dim,ipos-1] # g[i]-g[i-1]\n", " # linear interpolation for g(x) : \n", " xr[i,dim] = (grid.g[dim,ipos-1] + (pos-ipos)*diff)*unit # g(xr) ~ ( g[i-1]+(g[i]-g[i-1])*(x*N-i) )*[units]\n", " #gx = grid.g[dim,ipos-1] + (pos-ipos)*diff\n", " #xr[i,dim]= gx*(b-a) + a\n", " bins[i,dim]=ipos # remember in which bin this random number falls.\n", " weight *= diff*nbins # weight for this dimension is dg/dx = (g[i]-g[i-1])*N\n", " # because dx = i/N - (i-1)/N = 1/N\n", " wgh[i] = weight # total weight is (df/dx)*(df/dy)*(df/dx).../N_{samples}\n", " # Here we evaluate function f on all randomly generated x points above\n", " fx = integrant(xr) # n function evaluations required in single call\n", " neval += n # We just added so many fuction evaluations\n", " \n", " # Now we compute the integral as weighted average, namely, f(g(x))*dg/dx\n", " wfun = wgh * fx # weight * function ~ f_i*w_i, here w_i has 1/N in its weight, hence actually we are \n", " # evaluating 1/N * sum_i f_i*w_i\n", " cum.sum += sum(wfun) # 1/N sum_i f_i*w_i = \n", " wfun *= wfun # carefull : we need 1/N (f_i w_i)^2, while this gives (f_i w_i/N)^2\n", " cum.sqsum += sum(wfun) # sum_i (f_i*w_i/N)^2 = /all_nsamples\n", " # \n", " w0 = sqrt(cum.sqsum*all_nsamples) # w0 = sqrt()\n", " w1 = (w0 + cum.sum)*(w0 - cum.sum) # w1 = (w0^2 - ^2) = (-^2)\n", " w = (all_nsamples-1)/w1 # w ~ 1/sigma_i^2 = (N-1)/(-^2)\n", " # Note that variance of the MC sampling is Var(monte-f) = (-^2)/N == 1/sigma_i^2\n", " cum.weightsum += w # weightsum ~ \\sum_i 1/sigma_i^2\n", " cum.avgsum += w*cum.sum # avgsum ~ \\sum_i _i / sigma_i^2\n", " \n", " cum.avg = cum.avgsum/cum.weightsum # I_best = (\\sum_i _i/sigma_i^2 )/(\\sum_i 1/sigma_i^2)\n", " cum.err = sqrt(1/cum.weightsum) # err ~ sqrt(best sigma^2) = sqrt(1/(\\sum_i 1/sigma_i^2))\n", " \n", " cum.avg *= unit**ndim\n", " cum.err *= unit**ndim" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "\n", "unit=pi\n", "ndim=3\n", "maxeval=200000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "nbins=128\n", "nstart =10000\n", "nincrease=5000\n", "grid = Grid(ndim,nbins)\n", "random.seed(0)\n", "\n", "#SimpleMC(my_integrant2, ndim, pi, maxeval, cum)\n", "Vegas_step1(my_integrant2, pi, maxeval, nstart, nincrease, grid, cum)\n", "print (cum.avg, '+-', cum.err, 'exact=', exact)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are going to insert the sampling of the projection of the function $f(x)$ to all axis, $f_1(x), f_2(y)...$, from which we will calculate new grids $g_1(x), g_2(y)...$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Vegas_step2(integrant, unit, maxeval, nstart, nincrease, grid, cum):\n", " ndim, nbins = grid.ndim,grid.nbins # dimension of the integral, size of the grid for binning in each direction\n", " unit_dim = unit**ndim # converts from unit cube integration to generalized cube with unit length\n", " nbatch=1000 # function will be evaluated in bacthes of 1000 evaluations at one time (for efficiency and storage issues)\n", " neval=0\n", " print (\"\"\"Vegas parameters:\n", " ndim = \"\"\"+str(ndim)+\"\"\"\n", " unit = \"\"\"+str(unit)+\"\"\"\n", " maxeval = \"\"\"+str(maxeval)+\"\"\"\n", " nstart = \"\"\"+str(nstart)+\"\"\"\n", " nincrease = \"\"\"+str(nincrease)+\"\"\"\n", " nbins = \"\"\"+str(nbins)+\"\"\"\n", " nbaths = \"\"\"+str(nbatch)+\"\\n\")\n", "\n", " bins = zeros((nbatch,ndim),dtype=int) # in which sampled bin does this point fall?\n", " \n", " all_nsamples = nstart\n", " fxbin = zeros((ndim,nbins)) #new2: after each iteration we reset the average function being binned\n", " for nsamples in range(all_nsamples,0,-nbatch): # loop over all_nsample evaluations in batches of nbatch\n", " n = min(nbatch,nsamples) # How many evaluations in this pass?\n", " # We are integrating f(g_1(x),g_2(y),g_3(z))*dg_1/dx*dg_2/dy*dg_3/dz dx*dy*dz\n", " # This is represented as 1/all_nsamples \\sum_{x_i,y_i,z_i} f(g_1(x_i),g_2(y_i),g_3(z_i))*dg_1/dx*dg_2/dy*dg_3/dz\n", " # where dg_1/dx = diff*NBINS\n", " wgh = zeros(nbatch) # weights for each random point in the batch\n", " xr = random.random((n,ndim)) # generates 2-d array of random numbers in the interval [0,1)\n", " for i in range(n):\n", " weight = 1.0/all_nsamples\n", " for dim in range(ndim):\n", " # We want to evaluate the function f at point g(x), i.e, f(g_1(x),g_2(y),...)\n", " # Here we transform the points x,y,z -> g_1(x), g_2(y), g_3(z)\n", " # We hence want to evaluate g(x) ~ g(x[i]), where x is the random number and g is the grid function\n", " # The discretized g(t) is defined on the grid :\n", " # t[-1]=0, t[0]=1/N, t[1]=2/N, t[2]=3/N ... t[N-1]=1.\n", " # We know that g(0)=0 and g(1)=1, so that g[-1]=0.0 and g[N-1]=1.0\n", " # To interpolate g at x, we first compute i=int(x*N) and then we use linear interpolation\n", " # g(x) = g[i-1] + (g[i]-g[i-1])*(x*N-i) ; if i>0\n", " # g(x) = 0 + (g[0]-0)*(x*N-0) ; if i=0\n", " #\n", " pos = xr[i,dim]*nbins # which grid would it fit ? (x*N)\n", " ipos = int(pos) # the grid position is ipos : int(x*N)==i\n", " diff = grid.g[dim,ipos] - grid.g[dim,ipos-1] # g[i]-g[-1]\n", " # linear interpolation for g(x) : \n", " xr[i,dim] = (grid.g[dim,ipos-1] + (pos-ipos)*diff)*unit # g(xr) ~ ( g[i-1]+(g[i]-g[i-1])*(x*N-i) )*[units]\n", " bins[i,dim]=ipos # remember in which bin this random number falls.\n", " weight *= diff*nbins # weight for this dimension is dg/dx = (g[i]-g[i-1])*N\n", " # because dx = i/N - (i-1)/N = 1/N\n", " wgh[i] = weight # total weight is (df/dx)*(df/dy)*(df/dx).../N_{samples}\n", " # Here we evaluate function f on all randomly generated x points above\n", " fx = integrant(xr) # n function evaluations required in single call\n", " neval += n # We just added so many fuction evaluations\n", "\n", " # Now we compute the integral as weighted average, namely, f(g(x))*dg/dx\n", " wfun = wgh * fx # weight * function ~ f_i*w_i\n", " cum.sum += sum(wfun) # sum_i f_i*w_i = \n", " wfun *= wfun # carefull : this is like (f_i * w_i/N)^2 hence 1/N (1/N (f_i*w_i)^2)\n", " cum.sqsum += sum(wfun) # sum_i (f_i*w_i)^2 = /all_nsamples\n", " # \n", " for dim in range(ndim): #new2\n", " # projection of the sampled function^2 to each dimension, which will be used to improve the grid.\n", " for i in range(n): #new2\n", " fxbin[dim, bins[i,dim] ] += wfun[i] #new2: just bin the function f. We saved the bin position before.\n", " \n", " w0 = sqrt(cum.sqsum*all_nsamples) # w0 = sqrt()\n", " w1 = (w0 + cum.sum)*(w0 - cum.sum) # w1 = (w0^2 - ^2) = (-^2)\n", " w = (all_nsamples-1)/w1 # w ~ 1/sigma_i^2 = (N-1)/(-^2)\n", " # Note that variance of the MC sampling is Var(monte-f) = (-^2)/N == 1/sigma_i^2\n", " cum.weightsum += w # weightsum ~ \\sum_i 1/sigma_i^2\n", " cum.avgsum += w*cum.sum # avgsum ~ \\sum_i _i / sigma_i^2\n", " \n", " cum.avg = cum.avgsum/cum.weightsum # I_best = (\\sum_i _i/sigma_i^2 )/(\\sum_i 1/sigma_i^2)\n", " cum.err = sqrt(1/cum.weightsum) # err ~ sqrt(best sigma^2) = sqrt(1/(\\sum_i 1/sigma_i^2))\n", " \n", " cum.avg *= unit**ndim\n", " cum.err *= unit**ndim\n", " \n", " return fxbin" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "\n", "unit=pi\n", "ndim=3\n", "maxeval=200000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "nbins=128\n", "nstart =10000\n", "nincrease=5000\n", "grid = Grid(ndim,nbins)\n", "\n", "#SimpleMC(my_integrant2, ndim, pi, maxeval, cum)\n", "fxbin = Vegas_step2(my_integrant2, pi, maxeval, nstart, nincrease, grid, cum)\n", "print(cum.avg, '+-', cum.err, 'exact=', exact)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pylab import *\n", "%matplotlib inline\n", "\n", "plot(fxbin[0])\n", "plot(fxbin[1])\n", "plot(fxbin[2])\n", "xlim([0,128])\n", "ylim([0,1e-7])\n", "show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def smfun(x):\n", " if (x>0):\n", " return ((x-1.)/log(x))**(1.5)\n", " else:\n", " return 0.\n", "vsmfun = vectorize(smfun)\n", "\n", "def Smoothen(fxbin):\n", " (ndim,nbins) = shape(fxbin)\n", " final = zeros(shape(fxbin))\n", " for idim in range(ndim):\n", " fxb = copy(fxbin[idim,:])\n", " #**** smooth the f^2 value stored for each bin ****\n", " # f[i] <- (f[i+1]+f[i]+f[i-1])/3.\n", " fxb[:nbins-1] += fxbin[idim,1:nbins]\n", " fxb[1:nbins] += fxbin[idim,:nbins-1]\n", " fxb[1:nbins-1] *= 1/3.\n", " fxb[0] *= 1/2.\n", " fxb[nbins-1] *= 1/2.\n", " norm = sum(fxb)\n", " if( norm == 0 ):\n", " print ('ERROR can not refine the grid with zero grid function')\n", " return # can not refine the grid if the function is zero.\n", " fxb *= 1.0/norm # we normalize the function.\n", " # Note that normalization is such that the sum is 1.\n", " final[idim,:] = vsmfun(fxb)\n", " return final" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pylab import *\n", "%matplotlib inline\n", "\n", "imp = Smoothen(fxbin)\n", "\n", "plot(imp[0])\n", "plot(imp[1])\n", "plot(imp[2])\n", "xlim([0,128])\n", "show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class Grid:\n", " \"\"\"Contains the grid points g_n(x) with x=[0...1], and g=[0...1]\n", " for Vegas integration. There are n-dim g_n functions.\n", " Constraints : g(0)=0 and g(1)=1.\n", " \"\"\"\n", " def __init__(self, ndim, nbins):\n", " self.g = zeros((ndim,nbins+1)) \n", " # a bit dirty trick: We will later use also g[-1] in interpolation, which should be set to zero, hence\n", " # we allocate dimension nbins+1, rather than nbinx\n", " self.ndim=ndim\n", " self.nbins=nbins\n", " # At the beginning we set g(x)=x\n", " # The grid-points are x_0 = 1/N, x_1 = 2/N, ... x_{N-1}=1.0. \n", " # Note that g(0)=0, and we skip this point on the mesh.\n", " for idim in range(ndim):\n", " self.g[idim,:nbins] = arange(1,nbins+1)/float(nbins)\n", " \n", " def RefineGrid(self, imp):\n", " (ndim,nbins) = shape(imp)\n", " gnew = zeros((ndim,nbins+1))\n", " for idim in range(ndim):\n", " avgperbin = sum(imp[idim,:])/nbins\n", " #**** redefine the size of each bin ****\n", " newgrid = zeros(nbins)\n", " cur=0.0\n", " newcur=0.0\n", " thisbin = 0.0\n", " ibin = -1\n", " # we are trying to determine\n", " # Int[ f(g) dg, {g, g_{i-1},g_i}] == I/N_g\n", " # where I == avgperbin\n", " for newbin in range(nbins-1): # all but the last bin, which is 1.0\n", " while (thisbin < avgperbin) :\n", " ibin+=1\n", " thisbin += imp[idim,ibin]\n", " prev = cur\n", " cur = self.g[idim,ibin]\n", " # Explanation is in order : \n", " # prev -- g^{old}_{l-1}\n", " # cur -- g^{old}_l\n", " # thisbin -- Sm = f_{l-k}+.... +f_{l-2}+f_{l-1}+f_l\n", " # we know that Sm is just a bit more than we need, i.e., I/N_g, hence we need to compute how much more\n", " # using linear interpolation :\n", " # g^{new} = g_l - (g_l-g_{l-1}) * (f_{l-k}+....+f_{l-2}+f_{l-1}+f_l - I/N_g)/f_l\n", " # clearly\n", " # if I/N_g == f_{l-k}+....+f_{l-2}+f_{l-1}+f_l\n", " # we will get g^{new} = g_l\n", " # and if I/N_g == f_{l-k}+....+f_{l-2}+f_{l-1}\n", " # we will get g^{new} = g_{l-1}\n", " # and if I/N_g is between the two possibilities, we will get linear interpolation between\n", " # g_{l-1} and g_l\n", " # \n", " thisbin -= avgperbin # thisbin <- (f_{l-k}+....+f_{l-2}+f_{l-1}+f_l - I/N_g)\n", " delta = (cur - prev)*thisbin # delta <- (g_l-g_{l-1})*(f_{l-k}+....+f_{l-2}+f_{l-1}+f_l - I/N_g)\n", " # cur is the closest point from the old mesh, while delta/imp is the correction using linear interpolation.\n", " newgrid[newbin] = cur - delta/imp[idim,ibin] \n", "\n", " newgrid[nbins-1]=1.0\n", " gnew[idim,:nbins]= newgrid\n", " self.g = gnew\n", " return gnew" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Update the class Grid, so that it can refine the grid" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "\n", "unit=pi\n", "ndim=3\n", "maxeval=200000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "nbins=128\n", "nstart =10000\n", "nincrease=5000\n", "\n", "grid = Grid(ndim,nbins)\n", "\n", "fxbin = Vegas_step2(my_integrant2, pi, maxeval, nstart, nincrease, grid, cum)\n", "imp = Smoothen(fxbin)\n", "grid.RefineGrid(imp)\n", "print(cum.avg, '+-', cum.err, 'exact=', exact)\n", "\n", "plot(grid.g[0,:nbins])\n", "plot(grid.g[1,:nbins])\n", "plot(grid.g[2,:nbins])\n", "xlim([0,128])\n", "show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def Vegas_step3(integrant, unit, maxeval, nstart, nincrease, grid, cum):\n", " ndim, nbins = grid.ndim,grid.nbins # dimension of the integral, size of the grid for binning in each direction\n", " unit_dim = unit**ndim # converts from unit cube integration to generalized cube with unit length\n", " nbatch=1000 # function will be evaluated in bacthes of 1000 evaluations at one time (for efficiency and storage issues)\n", " neval=0\n", " print (\"\"\"Vegas parameters:\n", " ndim = \"\"\"+str(ndim)+\"\"\"\n", " unit = \"\"\"+str(unit)+\"\"\"\n", " maxeval = \"\"\"+str(maxeval)+\"\"\"\n", " nstart = \"\"\"+str(nstart)+\"\"\"\n", " nincrease = \"\"\"+str(nincrease)+\"\"\"\n", " nbins = \"\"\"+str(nbins)+\"\"\"\n", " nbaths = \"\"\"+str(nbatch)+\"\\n\")\n", "\n", " bins = zeros((nbatch,ndim),dtype=int) # in which sampled bin does this point fall?\n", " \n", " \n", " all_nsamples = nstart\n", " for iter in range(1000): # NEW in step 3\n", " wgh = zeros(nbatch) # weights for each random point in the batch\n", " fxbin = zeros((ndim,nbins)) # after each iteration we reset the average function being binned\n", " for nsamples in range(all_nsamples,0,-nbatch): # loop over all_nsample evaluations in batches of nbatch\n", " n = min(nbatch,nsamples) # How many evaluations in this pass?\n", " # We are integrating f(g_1(x),g_2(y),g_3(z))*dg_1/dx*dg_2/dy*dg_3/dz dx*dy*dz\n", " # This is represented as 1/all_nsamples \\sum_{x_i,y_i,z_i} f(g_1(x_i),g_2(y_i),g_3(z_i))*dg_1/dx*dg_2/dy*dg_3/dz\n", " # where dg_1/dx = diff*NBINS\n", " xr = random.random((n,ndim)) # generates 2-d array of random numbers in the interval [0,1)\n", " for i in range(n):\n", " weight = 1.0/all_nsamples\n", " for dim in range(ndim):\n", " # We want to evaluate the function f at point g(x), i.e, f(g_1(x),g_2(y),...)\n", " # Here we transform the points x,y,z -> g_1(x), g_2(y), g_3(z)\n", " # We hence want to evaluate g(x) ~ g(x[i]), where x is the random number and g is the grid function\n", " # The discretized g(t) is defined on the grid :\n", " # t[-1]=0, t[0]=1/N, t[1]=2/N, t[2]=3/N ... t[N-1]=1.\n", " # We know that g(0)=0 and g(1)=1, so that g[-1]=0.0 and g[N-1]=1.0\n", " # To interpolate g at x, we first compute i=int(x*N) and then we use linear interpolation\n", " # g(x) = g[i-1] + (g[i]-g[i-1])*(x*N-i) ; if i>0\n", " # g(x) = 0 + (g[0]-0)*(x*N-0) ; if i=0\n", " #\n", " pos = xr[i,dim]*nbins # which grid would it fit ? (x*N)\n", " ipos = int(pos) # the grid position is ipos : int(x*N)==i\n", " diff = grid.g[dim,ipos] - grid.g[dim,ipos-1] # g[i]-g[-1]\n", " # linear interpolation for g(x) : \n", " xr[i,dim] = (grid.g[dim,ipos-1] + (pos-ipos)*diff)*unit # g(xr) ~ ( g[i-1]+(g[i]-g[i-1])*(x*N-i) )*[units]\n", " bins[i,dim]=ipos # remember in which bin this random number falls.\n", " weight *= diff*nbins # weight for this dimension is dg/dx = (g[i]-g[i-1])*N\n", " # because dx = i/N - (i-1)/N = 1/N\n", " wgh[i] = weight # total weight is (df/dx)*(df/dy)*(df/dx).../N_{samples}\n", " \n", " # Here we evaluate function f on all randomly generated x points above\n", " fx = integrant(xr) # n function evaluations required in single call\n", " neval += n # We just added so many fuction evaluations\n", " \n", " # Now we compute the integral as weighted average, namely, f(g(x))*dg/dx\n", " wfun = wgh * fx # weight * function ~ f_i*w_i \n", " cum.sum += sum(wfun) # sum_i f_i*w_i = \n", " wfun *= wfun # carefull : this is like (f_i * w_i/N)^2 hence 1/N (1/N (f_i*w_i)^2)\n", " cum.sqsum += sum(wfun) # sum_i (f_i*w_i)^2 = /all_nsamples\n", " # \n", " for dim in range(ndim): #new2\n", " # Here we make a better approximation for the function, which we are integrating.\n", " for i in range(n): #new2\n", " fxbin[dim, bins[i,dim] ] += wfun[i] #new2: just bin the function f. We saved the bin position before.\n", " \n", " w0 = sqrt(cum.sqsum*all_nsamples) # w0 = sqrt()\n", " w1 = (w0 + cum.sum)*(w0 - cum.sum) # w1 = (w0^2 - ^2) = (-^2)\n", " w = (all_nsamples-1)/w1 # w ~ 1/sigma_i^2 = (N-1)/(-^2)\n", " # Note that variance of the MC sampling is Var(monte-f) = (-^2)/N == 1/sigma_i^2\n", " cum.weightsum += w # weightsum ~ \\sum_i 1/sigma_i^2\n", " cum.avgsum += w*cum.sum # avgsum ~ \\sum_i _i / sigma_i^2\n", " cum.avg2sum += w*cum.sum**2 # avg2cum ~ \\sum_i _i^2/sigma_i^2\n", " \n", " cum.avg = cum.avgsum/cum.weightsum # I_best = (\\sum_i _i/sigma_i^2 )/(\\sum_i 1/sigma_i^2)\n", " cum.err = sqrt(1/cum.weightsum) # err ~ sqrt(best sigma^2) = sqrt(1/(\\sum_i 1/sigma_i^2))\n", " \n", " # NEW in this step3\n", " if iter>0:\n", " cum.chisq = (cum.avg2sum - 2*cum.avgsum*cum.avg + cum.weightsum*cum.avg**2)/iter\n", " \n", " print (\"Iteration {:3d}: I= {:10.8f} +- {:10.8f} chisq= {:10.8f} number of evaluations = {:7d} \".format(iter+1, cum.avg*unit_dim, cum.err*unit_dim, cum.chisq, neval))\n", " imp = Smoothen(fxbin)\n", " grid.RefineGrid(imp)\n", " \n", " cum.sum=0 # clear the partial sum for the next step\n", " cum.sqsum=0\n", " all_nsamples += nincrease # for the next time, increase the number of steps a bit\n", " if (neval>=maxeval): break\n", " \n", " cum.avg *= unit**ndim\n", " cum.err *= unit**ndim" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "\n", "unit=pi\n", "ndim=3\n", "maxeval=2000000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "nbins=128\n", "nstart =100000\n", "nincrease=5000\n", "\n", "grid = Grid(ndim,nbins)\n", "\n", "random.seed(0)\n", "\n", "Vegas_step3(my_integrant2, pi, maxeval, nstart, nincrease, grid, cum)\n", "\n", "print (cum.avg, '+-', cum.err, 'exact=', exact, 'real error=', abs(cum.avg-exact)/exact)\n", "\n", "plot(grid.g[0,:nbins])\n", "plot(grid.g[1,:nbins])\n", "plot(grid.g[2,:nbins])\n", "show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we are going to speed up the code by using vectorized numpy capability, and numba capability.\n", "\n", "Here numba will be used to set fxbin array, using wfun array, which contains projection of the function integrated to all axis.\n", "\n", "To interpolate grid at the random points xr, we will use numpy vectorized functionality and fancy indexing to remove the loop over batches. All loops over batch size n is gone in the final version of the code." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from numba import jit\n", "\n", "@jit(nopython=True)\n", "def SetFxbin(fxbin,bins,wfun):\n", " (n,ndim) = bins.shape\n", " for dim in range(ndim):\n", " # Here we make a better approximation for the function, which we are integrating.\n", " for i in range(n):\n", " fxbin[dim, bins[i,dim] ] += abs(wfun[i]) # just bin the function f. We saved the bin position before. \n", "\n", "\n", "def Vegas_step3b(integrant, unit, maxeval, nstart, nincrease, grid, cum):\n", " ndim, nbins = grid.ndim,grid.nbins # dimension of the integral, size of the grid for binning in each direction\n", " unit_dim = unit**ndim # converts from unit cube integration to generalized cube with unit length\n", " nbatch=1000 # function will be evaluated in bacthes of 1000 evaluations at one time (for efficiency and storage issues)\n", " neval=0\n", " print (\"\"\"Vegas parameters:\n", " ndim = \"\"\"+str(ndim)+\"\"\"\n", " unit = \"\"\"+str(unit)+\"\"\"\n", " maxeval = \"\"\"+str(maxeval)+\"\"\"\n", " nstart = \"\"\"+str(nstart)+\"\"\"\n", " nincrease = \"\"\"+str(nincrease)+\"\"\"\n", " nbins = \"\"\"+str(nbins)+\"\"\"\n", " nbaths = \"\"\"+str(nbatch)+\"\\n\")\n", "\n", " bins = zeros((nbatch,ndim),dtype=int) # in which sampled bin does this point fall?\n", " \n", " \n", " all_nsamples = nstart\n", " for iter in range(1000): # NEW in step 3\n", " wgh = zeros(nbatch) # weights for each random point in the batch\n", " fxbin = zeros((ndim,nbins)) # after each iteration we reset the average function being binned\n", " for nsamples in range(all_nsamples,0,-nbatch): # loop over all_nsample evaluations in batches of nbatch\n", " n = min(nbatch,nsamples) # How many evaluations in this pass?\n", " # We are integrating f(g_1(x),g_2(y),g_3(z))*dg_1/dx*dg_2/dy*dg_3/dz dx*dy*dz\n", " # This is represented as 1/all_nsamples \\sum_{x_i,y_i,z_i} f(g_1(x_i),g_2(y_i),g_3(z_i))*dg_1/dx*dg_2/dy*dg_3/dz\n", " # where dg_1/dx = diff*NBINS\n", " xr = random.random((n,ndim)) # generates 2-d array of random numbers in the interval [0,1)\n", " pos = xr*nbins # (x*N)\n", " bins = array(pos,dtype=int) # which grid would it fit ? (x*N)\n", " wgh = ones(nbatch)/all_nsamples\n", " for dim in range(ndim): \n", " # We want to evaluate the function f at point g(x), i.e, f(g_1(x),g_2(y),...)\n", " # Here we transform the points x,y,z -> g_1(x), g_2(y), g_3(z)\n", " # We hence want to evaluate g(x) ~ g(x[i]), where x is the random number and g is the grid function\n", " # The discretized g(t) is defined on the grid :\n", " # t[-1]=0, t[0]=1/N, t[1]=2/N, t[2]=3/N ... t[N-1]=1.\n", " # We know that g(0)=0 and g(1)=1, so that g[-1]=0.0 and g[N-1]=1.0\n", " # To interpolate g at x, we first compute i=int(x*N) and then we use linear interpolation\n", " # g(x) = g[i-1] + (g[i]-g[i-1])*(x*N-int(x*N))\n", " gi = grid.g[dim,bins[:,dim]] # g[i]\n", " gm = grid.g[dim,bins[:,dim]-1] # g[i-1]\n", " diff = gi - gm # g[i]-g[i-1]\n", " gx = gm + (pos[:,dim]-bins[:,dim])*diff # linear interpolation\n", " xr[:,dim] = gx*unit # xr <- g(xr)\n", " wgh *= diff*nbins # wgh = prod_{dim} dg/dx\n", " \n", " # Here we evaluate function f on all randomly generated x points above\n", " fx = integrant(xr) # n function evaluations required in single call\n", " neval += n # We just added so many fuction evaluations\n", " \n", " # Now we compute the integral as weighted average, namely, f(g(x))*dg/dx\n", " wfun = wgh * fx # weight * function ~ f_i*w_i \n", " cum.sum += sum(wfun) # sum_i f_i*w_i = \n", " wfun *= wfun # carefull : this is like (f_i * w_i/N)^2 hence 1/N (1/N (f_i*w_i)^2)\n", " cum.sqsum += sum(wfun) # sum_i (f_i*w_i)^2 = /all_nsamples\n", " # \n", " SetFxbin(fxbin,bins,wfun)\n", " #for dim in range(ndim): #new2\n", " # # Here we make a better approximation for the function, which we are integrating.\n", " # for i in range(n): #new2\n", " # fxbin[dim, bins[i,dim] ] += wfun[i] #new2: just bin the function f. We saved the bin position before.\n", " \n", " w0 = sqrt(cum.sqsum*all_nsamples) # w0 = sqrt()\n", " w1 = (w0 + cum.sum)*(w0 - cum.sum) # w1 = (w0^2 - ^2) = (-^2)\n", " w = (all_nsamples-1)/w1 # w ~ 1/sigma_i^2 = (N-1)/(-^2)\n", " # Note that variance of the MC sampling is Var(monte-f) = (-^2)/N == 1/sigma_i^2\n", " cum.weightsum += w # weightsum ~ \\sum_i 1/sigma_i^2\n", " cum.avgsum += w*cum.sum # avgsum ~ \\sum_i _i / sigma_i^2\n", " cum.avg2sum += w*cum.sum**2 # avg2cum ~ \\sum_i _i^2/sigma_i^2\n", " \n", " cum.avg = cum.avgsum/cum.weightsum # I_best = (\\sum_i _i/sigma_i^2 )/(\\sum_i 1/sigma_i^2)\n", " cum.err = sqrt(1/cum.weightsum) # err ~ sqrt(best sigma^2) = sqrt(1/(\\sum_i 1/sigma_i^2))\n", " \n", " # NEW in this step3\n", " if iter>0:\n", " cum.chisq = (cum.avg2sum - 2*cum.avgsum*cum.avg + cum.weightsum*cum.avg**2)/iter\n", " \n", " print (\"Iteration {:3d}: I= {:10.8f} +- {:10.8f} chisq= {:10.8f} number of evaluations = {:7d} \".format(iter+1, cum.avg*unit_dim, cum.err*unit_dim, cum.chisq, neval))\n", " imp = Smoothen(fxbin)\n", " grid.RefineGrid(imp)\n", " \n", " cum.sum=0 # clear the partial sum for the next step\n", " cum.sqsum=0\n", " all_nsamples += nincrease # for the next time, increase the number of steps a bit\n", " if (neval>=maxeval): break\n", " \n", " cum.avg *= unit**ndim\n", " cum.err *= unit**ndim" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import *\n", "from numpy import *\n", "\n", "unit=pi\n", "ndim=3\n", "maxeval=2000000\n", "exact = 1.3932 # exact value of the integral\n", " \n", "cum = Cumulants()\n", " \n", "nbins=128\n", "nstart =100000\n", "nincrease=5000\n", "\n", "grid = Grid(ndim,nbins)\n", "\n", "random.seed(0)\n", "\n", "Vegas_step3b(my_integrant2, pi, maxeval, nstart, nincrease, grid, cum)\n", "\n", "print (cum.avg, '+-', cum.err, 'exact=', exact, 'real error=', abs(cum.avg-exact)/exact)\n", "\n", "plot(grid.g[0,:nbins])\n", "plot(grid.g[1,:nbins])\n", "plot(grid.g[2,:nbins])\n", "show()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Homework\n", "\n", "* generalize Vegas algorithm so that the integration limits are given by arbitrary numbers [a,b]. You can still assume supercube in which all dimensions have the same limit [a,b].\n", "* Test Vegas on the same example $f=1/(1-\\cos(k_x)\\cos(k_y)\\cos(k_z))/\\pi^3$ but use limits $[-\\pi,\\pi]$ instead of $[0,\\pi]$.\n", "* Speed up the part of the code that Redefines the grid RefineGrid\n", "* generalize Vegas so that it works for complex functions.\n", "* Using Vegas, evaluate the following Linhardt function:\n", "\n", "\\begin{eqnarray}\n", "P(\\Omega,\\vec{q}) = -2 \\int \\frac{d^3k}{(2\\pi)^3} \\frac{f(\\varepsilon_{\\vec{k}+\\vec{q}})-f(\\varepsilon_{\\vec{k}})}{\\Omega-\\varepsilon_{\\vec{k}+\\vec{q}}+\\varepsilon_\\vec{k}+i\\delta}\n", "\\end{eqnarray}\n", "\n", "Here $\\varepsilon_\\vec{k}=k^2-k_F^2$.\n", "\n", "We will use $k_F=\\frac{(\\frac{9\\pi}{4})^{1/3}}{rs}$ with $r_s=2$, $f(x) = 1/(\\exp(x/T)+1)$, $T=0.02 k_F^2$, $\\delta = 0.002 k_F^2$, $\\Omega=0$, $q=0.1 k_F$, integration limits can be set to $[-3 k_F,3 k_F]$.\n", "\n", "The result for $P(\\Omega=0,q<< k_F)$ should be close to $-n_F$, where $n_F = \\frac{k_F}{2\\pi^2}$.\n", "\n", "* Optional: Change the Vegas algorithm so that it computes $P(\\Omega,\\vec{q})$ at once for an array of $\\Omega$ points, such as linspace(0,0.5*kF*kF,200). Use $P(\\Omega=0)$ to redefine the grid, so that we have most efficient integration at $\\Omega=0$ (where the function is hardest to integrate)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }