{ "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 }