Adding David Evans notebook on Parallel Computing with updates to Pyt…

mpi/MPI_Allgather.py

@@ -0,0 +1,20 @@
+
+import numpy as np
+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+# create a 2x2 array filled with the rank of the process
+x = np.ones((2,2)) * rank 
+
+# create an empty array to store the gathered result
+X = np.empty((2*size,2))
+
+#use numpy Allgather
+                #send buffer   #recv buffer
+comm.Allgather([x, MPI.DOUBLE],[X, MPI.DOUBLE])
+
+if rank ==0:
+    print(X)

mpi/MPI_montecarlo.py

@@ -0,0 +1,16 @@
+
+import numpy as np
+from mpi4py import MPI
+
+N = 10000
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+
+x = np.random.randn(N)
+my_result = np.cos(x).mean() #my average
+
+#now average accross all processes
+results = comm.gather(my_result)
+if rank == 0:
+    print(np.mean(results))

mpi/MPI_montecarlo_allgather.py

@@ -0,0 +1,21 @@
+
+import numpy as np
+from mpi4py import MPI
+
+N = 10000
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+
+x = np.random.randn(N)
+my_result = np.cos(x).mean() #my average
+
+#now gather from all processes
+results = comm.gather(my_result)
+
+print("Gather result for process " + str(rank) +":",results)
+
+#now allgather from all processes
+results = comm.allgather(my_result)
+
+print("Allgather resuts for process " + str(rank) + ":",results)

mpi/Main.py

@@ -0,0 +1,39 @@
+
+"""
+Created on Tue Jun 25 16:50:32 2013
+
+@author: dgevans
+"""
+
+from numpy import *
+from primitives import primitives_CRRA
+import bellman
+from mpi4py import MPI
+
+rank = MPI.COMM_WORLD.Get_rank()
+
+#initialize the parameters
+Para = primitives_CRRA()
+r = 0.9995*(1/Para.beta-1)
+
+#create the grid
+agrid = hstack((linspace(Para.a_min,Para.a_min+1.,15),linspace(Para.a_min+2.,Para.a_max,15)))
+zgrid = linspace(Para.z_min,Para.z_max,10)
+Para.domain = vstack([(a,z) for a in agrid for z in zgrid])
+
+#initial value function
+def V0(a,z):
+    return Para.U(exp(z)+r*a)/(1-Para.beta)
+
+#construct the Bellman Map
+T = bellman.BellmanMap(Para,r)
+
+#iterate on the Bellman Map until convergence
+Vf,Vs = bellman.approximateValueFunction(T(V0),Para)
+while diff > 1e-6:  
+    Vfnew,Vsnew = bellman.approximateValueFunction(T(Vf),Para)
+    diff =  max(abs(Vs-Vsnew))
+    if rank == 0:
+        print diff
+    Vf = Vfnew
+    Vs = Vsnew

mpi/bellman.py

@@ -0,0 +1,115 @@
+
+"""
+Created on Tue Jun 25 13:40:30 2013
+Holds all of the code for the bellman equation
+@author: dgevans
+"""
+
+from scipy.optimize import minimize_scalar
+
+from numpy import * 
+from mpi4py import MPI
+from scipy.interpolate import SmoothBivariateSpline
+from numpy.polynomial.hermite import hermgauss
+
+#compute nodes for gaussian quadrature.
+zhatvec,z_weights = hermgauss(10)
+z_weights /= sum(z_weights)
+
+
+def approximateValueFunction(TV,Para):
+    '''
+    Approximates the value function over the grid defined by Para.domain.  Uses
+    mpi. Returns both an interpolated value function and the value of TV at each point
+    in the domain.
+    '''
+    comm = MPI.COMM_WORLD
+    #first split up domain for each process
+    s = comm.Get_size()
+    rank = comm.Get_rank()
+    n = len(Para.domain)
+    m = n//s
+    r = n%s
+    #let each process take a slice of the domain
+    mydomain = Para.domain[rank*m+min(rank,r):
+                           (rank+1)*m+min(rank+1,r)]
+
+    #get the value at each point in my domain
+    myV = hstack(map(TV,mydomain))
+    #gather the values for each process
+    Vs = comm.gather(myV)
+
+    if rank == 0:
+        #fit 
+        Vs = hstack(Vs).flatten()
+        Vf = SmoothBivariateSpline(Para.domain[:,0],Para.domain[:,1],Vs)
+    else:
+        Vf = None
+    return comm.bcast(Vf),comm.bcast(Vs)
+
+
+
+
+
+class BellmanMap(object):
+    '''
+    Bellman map object.  Once constructed will take a value function and return
+    a new value function
+    '''
+
+    def __init__(self,Para,r):
+        '''
+        Initializes the Bellman Map function.  Need parameters and r
+        '''
+        self.Para = Para
+        self.r = r
+        self.w = 1.
+
+    def __call__(self,Vf):
+        '''
+        Given a current value function return new value function
+        '''
+        self.Vf = Vf
+
+        return lambda x: self.maximizeObjective(x[0],x[1])[0]
+
+
+    def maximizeObjective(self,a,z):
+        '''
+        Maximize the objective function Vf given the states a and productivity z.
+        Note the state for this problem are assets and previous periods productivity
+        return tuple (V,c_policy,a_policy)
+        '''
+
+        #total wealth for agent is easy
+        W = (1+self.r)*a + exp(z)*self.w
+        beta = self.Para.beta
+        Vf = self.Vf
+
+
+        def obj_f(aprime):
+            c = W-aprime
+            return -(self.Para.U(c) + beta*self.E(lambda zprime:Vf(aprime,zprime).flatten(),z))
+
+        a_min = self.Para.a_min
+        a_max = min(self.Para.a_max,W-0.00001)
+
+        res = minimize_scalar(obj_f,bounds=(a_min,a_max),method='bounded')
+
+        return -res.fun,W-res.x,res.x
+
+    def E(self,f,z_):
+        '''
+        Compute the exepected value of a function f of zprime conditional on z
+        '''
+        #Pz_min,Pz_max,z_pdf = self.get_z_distribution(z_)
+        z_min = array([self.Para.z_min])
+        z_max = array([self.Para.z_max])
+        #Compute the nodes for Gaussian quadrature
+        zvec = self.Para.rho * z_ + self.Para.sigma_e*zhatvec*sqrt(2)
+
+        #lower bound for z
+        zvec[zvec<z_min] = z_min 
+        #upper bound for z
+        zvec[zvec>z_max] = z_max
+        return z_weights.dot(f(zvec))

mpi/mpi_helloworld.py

@@ -0,0 +1,8 @@
+
+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD #retreive the communicator module
+rank = comm.Get_rank() #get the rank of the process
+size = comm.Get_size() #get the number of processes
+
+print("Hello world from process "+str(rank)+" of " + str(size))

mpi/primitives.py

@@ -0,0 +1,54 @@
+
+"""
+Created on Tue Jun 25 13:29:54 2013
+
+@author: dgevans
+"""
+from numpy import *
+
+class primitives(object):
+    '''
+    Class that holds the basic parameters of the problem
+    '''
+
+    beta = 0.95#discount factor
+
+    rho = 0.9#persistence of idiosyncratic productivity
+
+    sigma_e = 0.1#standard deviation of idiosyncratic productivity
+
+    a_min = 0.0 #minimum assets
+
+    a_max = 60.0 #maximum assets
+
+    z_min = -2. #minimum productivity
+
+    z_max = 2. #maximum productivity
+
+
+
+class primitives_CRRA(primitives):
+    '''
+    Extension of primitives class holding the preference structure
+    '''
+
+    sigma = 2.
+
+    def U(self,c):
+        """
+        CRRA utility function
+        """
+        sigma = self.sigma
+        if sigma == 1.:
+            return log(c)
+        else:
+            return (c)**(1-sigma)/(1-sigma)
+
+    def Uc(self,c):
+        """
+        Derivative of the CRRA utility function
+        """
+        sigma = self.sigma
+        return c**(-sigma)
+
+