From: W. Trevor King Date: Thu, 16 Sep 2010 04:03:12 +0000 (-0400) Subject: Remove old temporary files from content/monte_carlo/. X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=73b0d6fa90d59fa9af2c668b4ec0f128f390c756;p=parallel_computing.git Remove old temporary files from content/monte_carlo/. --- diff --git a/content/monte_carlo/index.shtml-blahtex b/content/monte_carlo/index.shtml-blahtex deleted file mode 100644 index 7d2a74b..0000000 --- a/content/monte_carlo/index.shtml-blahtex +++ /dev/null @@ -1,578 +0,0 @@ - - -

Monte Carlo Techniques

- - - - -

The Monte Carlo technique to perform integrals applies best to -higher dimensional integrals for which other traditional methods -become difficult. For instance, using few hundred points to perform a -1-D integral using the trapezoidal or Simpson method is easily done. -However, the same in multi-dimension does not apply: a five -dimensional integral

- - - - I=∫dx1∫dx2∫dx3∫dx4∫dx5f(x1,x2,x3,x4,x5) - - -

would require 2005 function evaluations to implement a -method with 200 points per axis. Statistical Physics requires -integrals of much higher dimensions, often of the order of the number -of particles.

- -

The Monte Carlo approach deals with this problem by sampling the -integrand via statistical means. The method calls for a random walk, -or a guided random walk, in phase space during which the integrand is -evaluated at each step and averaged over.

- -

This problem is presented here as a challenging case to implement -via a parallel algorithm. Different implementations will require us -to use all the MPI functions we have seen so far. In the following, -the Monte Carlo method will first be explained, then it will be -implemented serially, and then the Metropolis algorithm will be -explained, and implemented serially and in parallel.

- -

Monte Carlo Integral

- -

The Monte Carlo Method is best explained via a 1-D integral

- - - - I=∫abf(x)dx - - -

Rewrite

- - - - I=(b-a)〈f(x)〉 - - -

where 〈f(x)〉 is the average of f(x) over the -interval (a,b).

- - - - 〈f(x)〉=1N∑i=1Nf(xi) - - -

The Monte Carlo method consists in choosing the sampling points -xi at random in the interval (a,b).

- - - -

The advantage over other methods in higher dimension is that better -coverage of the function can be achieved with fewer points in phase -space.

- -

The error in evaluating 〈f(x)〉 is described by the -variance

- - - - σ2=1N[〈f2〉-〈f〉2] - - -

Calculation of π

- -

The method is illustrated by integrating f(x,y)=1 in a -circle, namely

- - - - I=∫-11dx∫-11dyΘ(1-x2-y2) - - -

This corresponds to

- - - < !--annotation encoding="LaTeX"> - I = \int_{-1}^{1} dx \int_{-1}^{1} dy - \begin{cases} - 1 & \text{if $x^2+y^2 \leq 1$} \\ - 0 & \text{if $x^2+y^2 \gt 1} - \end{cases} - - I = \int_{-1}^{1} dx \int_{-1}^{1} dy - \left[ \begin{array}{ll} - 1 & if x^2+y^2 \leq 1 \\ - 0 & if x^2+y^2 > 1 - \end{array} \right. - - I=∫-11dx∫-11dy[1ifx2+y2≤10ifx2+y2>1 - - -

< ! - - MATH: $ --"$
- -The Monte Carlo approach is implemented by enclosing the circle -in a square of side 2 and using random numbers in the interval -x = (-1,1) and y = (-1,1). The area of the circle -is π 12 and the area of the square is -22. The ratio of these areas is then π/4.

- -

Random numbers, taken as x - y coordinates of points in a -plane evidently fill in the plane uniformly. Therefore, the -numbers of random numbers that fall within the circle and the -square are proportional to their areas. Therefore

- -

    π = 4 # hits in circle  /  -# random numbers

- -

The program random_demo.c -generates x and y coordinates in the interval (0,1) via calls to -RAND( ), the build-in pseudo random-number generator function in -the "C" compiler. The code is serial and can run on any -computer.

- -

random_demo 2000 > data_file

- -

runs the random_demo to produce 2000 random points which are -piped in a file data_file. The latter can be plotted -using gnuplot. It produces the following graph.

- -

 
-The program monte_pi_serial.c calculates -π by the Monte-Carlo method. It uses random numbers in the -interval [0,1].

- -

The error in calculating π follows from the standard -deviation calculation in the following way:

- -

< ! - - MATH: $\sigma^2 = \frac{1}{N} [ < f^2 > - < f >^2 ]$ - - > --"$
- --"$ -
-  -"$
-  -"$

-
-  Back to top of page -
- -

Applicability of -the Method

- -

The random sampling of a function clearly works best if the -function is a constant.  (Why don't you then get the exact -value of π in the example above?) If the function varies -rapidly, or is locally peaked, the Monte Carlo method will not -converge as rapidly.  What is required is a "guidance" of -the random sampling that emphasizes sampling where the function -is large or varies rapidly.  Let's illustrate the concept in -1-D.

- -

< ! - - MATH: $I = \int_0^1 f(x) dx$ - - > -$ I = \int_0^1 f(x) dx $
-Introduce a weight function w(x)

- -

< ! - - MATH: $I = \int_0^1 w(x) \frac{f(x)}{w(x)} dx$ - - > --"$
-where w(x) is required to satisfy

- -

< ! - - MATH: $\int_0^1 w(x) dx = 1$ - - > -$ \int_0^1 w(x) dx = 1 $
-Change variable to y

- -

< ! - - MATH: $y(x) = \int_0^{x'} w(x') dx'$ - - > --"$
-such that

- -

< ! - - MATH: $dy = w(x) dx$ - - > -dy = w(x) dx, y(0) = -0, and y(1) = 1.

- -

-< ! - - MATH: $I = \int_0^1 dy \frac{ f( x(y) ) } { w( x(y) }$ - - > --"$
-The weight function w(x) is chosen such that [f/w] -is aproproximately constant. The Monte Carlo sampling of -[f/w] with uniform random sampling in y variable -will then be accurate.

- -

-< ! - - MATH: $I = \frac{1}{N} \sum_{i=1}^{N} \frac{ f( x(y_i) ) } { w( x(y_i) }$ - - > --"$
- -Koonin and Meridith (Computational Physics, 1989) give the -example

- -

< ! - - MATH: $I = \int_0^1 dx \frac{1}{1+x^2}$ - - > --"$
-with

- -

< ! - - MATH: $w(x) = \frac{1}{3} ( 4 - 2 x )$ - - > --"$
-In this simple case, the inversion x(y) is possible -analytically.

- -

< ! - - MATH: $y = \frac{1}{3} x ( 4 - x )$ - - > --"$
-  x = 2 - (4 -3 y )

- -

In general, this inversion x(y) is difficult or impossible to -perform analytically. An alternate approach must be used.

- -

The relation  dy = w(x) -dx  implies another interpretation of the -method.  Random numbers uniformly distributed in y -will be distributed according to w(x) in x. -Therefore, the guiding action of w(x) can be achieved by -generating random numbers according the the probability profile -w(x). This holds in phase spaces of arbitrary dimension. -The problem translates into that of finding random distribution -of arbitrary probability distributions.

- -

Uniform random numbers in interval (0,1) have an average value -of  0.5.  A binning in intervals dx will show an -average number of random number

- -

number / bin =  total # of random numbers / number of -bins.

- -

The number of random numbers in the bins will fluctuate around -the average with a Gaussian distribution.

- -

Random numbers distributed according to the distribution -w(x) will bin according to w(x).

-
-  Back to top of page -
- -

Metropolis Algorithm

- -

The Metropolis algorithm generates a random walk guided by the -weight function w(x).  The algorithm in 2-D involves -the following steps.  Starting from a point -(xi, yi).

- -
    -
  1. choose  η (step size)
    2. - -
  2. generate 2 uniformly distributed random numbers, - R1 and R2, in the interval - (-1,1)
    3. - -
  3. calculate a new trial point according to
    - < ! - - MATH: $x_{i+1}^T = x_i + \delta x_R$ - - > - -
    - < ! - - MATH: $y_{i+1}^T = y_i + \delta y_R$ - - > - -
    - where the steps in x and y are η - R1 and η R2 - respectively
    4. - -
  4.  evaluate the ratio of the weight function at the - trial point to the one at the old point
    - < ! - - MATH: $r = \frac{ w( x_{i+1}^T, y_{i+1}^T ) } { w( x_{i}, y_{i} -) }$ - - > - -
    5. - -
  5. if r>1, accept the new point - (xi+1,yi+1) = - (xi+1t,yi+1t) - and loop back to step (2)
    6. - -
  6. if r<1, accept the new point with probability - r.  That is, generate R, a random number in - the interval (0,1), and accept - (xi+1t,yi+1t) - as the new point (xi+1,yi+1) if - r > R and loop back to step (2)
    7. - -
  7. otherwise, reject it, and keep the previous point as the - new one (xi+1,yi+1) = - (xi,yi)
    8. - -
  8. loop back to step (2)
    9. -
- -

Koonin and Mederith (Computational Physics, 1989) provide a -simple proof that the algorithm produces the desired random -number distribution w(x).

- -

The step η  in the Metropolis algorithm is -arbitrary.  If chosen small, the acceptance ratio will be -high (a point in the neighborhood of an acceptable point will -also likely be acceptable).  But overall coverage of the -phase space will be very slow.  If η is large, -the phase space will be efficiently covered, but a new point will -seldom be accepted by the Metropolis criteria.  In either -case ineficiency follows. An empirical rule is to choose -η such as to have an acceptance ratio about 0.5 - -0.6.

- -

A concern about the algorithm is that the random numbers -generated by the Metropolis algorithm are correlated.   A -way to guard against this (Koonin and Mederith , Computational -Physics, 1989)  is to calculate the correlation -functions

- -

-< ! - - MATH: $c(k) = \frac{ - ^2 } { - ^2 }$ - - > --"$
- -and select each kth random number -generated by the algorithm with k selected to produce a -small value of Ck.

-
-  Back to top of page -
- -

Guided Monte Carlo -Integral

- -

Lets consider the integral
-  -"$
-The integral is peaked (torus shape) around the -origin.  
-
-The Gaussian function

- -

< ! - - MATH: $w(x,y) = \frac{1}{\pi} e^{ - (x^2+y^2) }$ - - > --"$
-can be chosen to guide the random number selection.

- -

The program Metropolis_Gaussian_Demo.c -  illustrates a random numbers generation according to -w(x,y). Note that arbitrary values for the step size -(η=0.2) and k (k=1) were chosen; no effort was -made to chose them optimally.

- -

This code, with 2000 points, produced the following random -number distribution.

- -


-The program Metropolis_Gaussian_Integral.c -  performs the integral.

- -

Metropolis_Gaussian_Integral -p > data_file
-    or
-Metropolis_Gaussian_Integral

- -

This code prompts the user for more points to be added to the -Monte Carlo integral. The calculation of the error may show the -calculation not to be accurate given that the choices of -η and k were not optimized.

-
-  Back to top of page -
- -

Parallel Monte Carlo -Algorithm

- -

Random Numbers Generation

- -

Pseudo random numbers are often generated by a linear -congruent algorithm.  The integer random number -xi+1 is generated from the previous one -by
-  xi+1 = ( -xi a + b ) mod -c  
-The coefficients a, b and c are chosen to -generate statistically good random numbers in a specific integer -range (Numerical Recipes).  Often the random numbers are -generated using different linear congruence relation, and mixed -via random selection as in RAND1 and RAND2 (Numerical -Recipes).

- -

The procedure to generate random numbers is -deterministic.  Two different nodes on a parallel machine -will produce exactly the same sequence if "seeded" in the same -way. The question then arises as to how to generate "good" random -numbers on the different nodes of the parallel computer that will -remain meaningful when used as a single set.

- -

One solution is to generate sequences on each node based on -different seeds.  The code random.demo_2_series.c -simulates this approach by generating 1 or 2 sequences of random -numbers that can be analyzed for their goodness.  The idea -is to look for the goodness of a sequence in which the first and -second half of the series are generated with different seed -versus sequences generated from a single seed.

- -

random_demo_2_series 2 1000 > data_file

- -

The following filters can be used to acquire some feelings for -the goodness of these random numbers

- - - -

Clearly the random numbers sequences generated from sequences -using different seeds have different properties than sequences -generated from single seed. 

- -

We use here the approach that Random numbers should be -generated on one node and distributed to the nodes which need -them.  This avoids the use of sequences generated from -different seeds.  Realistic applications often require -time-consuming function evaluations. Therefore the inefficiency -inherent in our approach is dwarfed by the large function -evaluation times.  This is also the approach used by  -Gropp et al. in "Using MPI".

-
-  Back to top of page -
- -

Parallel -Monte Carlo Codes

- -

A parallel implementation of the Monte Carlo integration must -have the random sequence distributed to the different nodes, -either by communication or by a local generation of the -sequence.  We implement here three versions of the code.

- -

Master-Workers paradigms

- -

This code generates the random numbers on nodes 0 and passes -them onto the worker nodes.  The latter compute the Monte -Carlo sub-sums and pass these sub-sums to node 0 which finishes -the sums.  The MPI_Reduce() function can be used to expedite -the global sums

- -

- -

Metropolis_Gaussian_Parallel.c -implements this approach.

- -

Client-Server paradigm

- -

This code implements the concept of a random number server; -this is a node reserved to serve random numbers with weight -w(x,y) to whichever node requests it.  The code sets node 0 -as a master node to direct the calculation and interface with the -user.  Node (numprocs-1) is the random numbers -server, and nodes 1 to (numprocs-2) are the worker -nodes which complete the Monte Carlo sub-sums.

- -

The code Metropolis_CS.c implements a version -of of this approach in which the functionality of the different -nodes is relegated to C functions. This illustrates the fact that -variables are local to the nodes. - -

Gropp et al, in "Using MPI" implements this version of a Monte -Carlo integration to calculate π.

- -

Full Workers -version

- -

In this version, node 0 administers the ranges of random -numbers to be generated for each workers, but leaves the actual -generation of the random numbers to the nodes themselves.  -Upon receiving a range, i.e., random numbers range 1001 to 2000, -the worker nodes generate these numbers and calculate the Monte -Carlo sub-sums.  This implies that they must generate -and discard those random numbers prior to their range -selection.  This is a waste of computing cycles, but given -the speed of computing versus the network speed it maybe the most -efficient approach after all.

- -

The scheme is therefore

- - - -

It is left to the reader to implement this approach.

- - diff --git a/content/monte_carlo/x b/content/monte_carlo/x deleted file mode 100644 index 773d413..0000000 --- a/content/monte_carlo/x +++ /dev/null @@ -1,394 +0,0 @@ -

Applicability of -the Method

- -

The random sampling of a function clearly works best if the -function is a constant.  (Why don't you then get the exact -value of π in the example above?) If the function varies -rapidly, or is locally peaked, the Monte Carlo method will not -converge as rapidly.  What is required is a "guidance" of -the random sampling that emphasizes sampling where the function -is large or varies rapidly.  Let's illustrate the concept in -1-D.

- -

< ! - - MATH: $I = \int_0^1 f(x) dx$ - - > -$ I = \int_0^1 f(x) dx $
-Introduce a weight function w(x)

- -

< ! - - MATH: $I = \int_0^1 w(x) \frac{f(x)}{w(x)} dx$ - - > --"$
-where w(x) is required to satisfy

- -

< ! - - MATH: $\int_0^1 w(x) dx = 1$ - - > -$ \int_0^1 w(x) dx = 1 $
-Change variable to y

- -

< ! - - MATH: $y(x) = \int_0^{x'} w(x') dx'$ - - > --"$
-such that

- -

< ! - - MATH: $dy = w(x) dx$ - - > -dy = w(x) dx, y(0) = -0, and y(1) = 1.

- -

-< ! - - MATH: $I = \int_0^1 dy \frac{ f( x(y) ) } { w( x(y) }$ - - > --"$
-The weight function w(x) is chosen such that [f/w] -is aproproximately constant. The Monte Carlo sampling of -[f/w] with uniform random sampling in y variable -will then be accurate.

- -

-< ! - - MATH: $I = \frac{1}{N} \sum_{i=1}^{N} \frac{ f( x(y_i) ) } { w( x(y_i) }$ - - > --"$
- -Koonin and Meridith (Computational Physics, 1989) give the -example

- -

< ! - - MATH: $I = \int_0^1 dx \frac{1}{1+x^2}$ - - > --"$
-with

- -

< ! - - MATH: $w(x) = \frac{1}{3} ( 4 - 2 x )$ - - > --"$
-In this simple case, the inversion x(y) is possible -analytically.

- -

< ! - - MATH: $y = \frac{1}{3} x ( 4 - x )$ - - > --"$
-  x = 2 - (4 -3 y )

- -

In general, this inversion x(y) is difficult or impossible to -perform analytically. An alternate approach must be used.

- -

The relation  dy = w(x) -dx  implies another interpretation of the -method.  Random numbers uniformly distributed in y -will be distributed according to w(x) in x. -Therefore, the guiding action of w(x) can be achieved by -generating random numbers according the the probability profile -w(x). This holds in phase spaces of arbitrary dimension. -The problem translates into that of finding random distribution -of arbitrary probability distributions.

- -

Uniform random numbers in interval $[0,1]$ have an average value -of  0.5.  A binning in intervals dx will show an -average number of random number

- -

number / bin =  total # of random numbers / number of -bins.

- -

The number of random numbers in the bins will fluctuate around -the average with a Gaussian distribution.

- -

Random numbers distributed according to the distribution -w(x) will bin according to w(x).

-
-  Back to top of page -
- -

Metropolis Algorithm

- -

The Metropolis algorithm generates a random walk guided by the -weight function w(x).  The algorithm in 2-D involves -the following steps.  Starting from a point -(xi, yi).

- -
    -
  1. choose  η (step size)
    2. - -
  2. generate 2 uniformly distributed random numbers, - R1 and R2, in the interval - (-1,1)
    3. - -
  3. calculate a new trial point according to
    - < ! - - MATH: $x_{i+1}^T = x_i + \delta x_R$ - - > - -
    - < ! - - MATH: $y_{i+1}^T = y_i + \delta y_R$ - - > - -
    - where the steps in x and y are η - R1 and η R2 - respectively
    4. - -
  4.  evaluate the ratio of the weight function at the - trial point to the one at the old point
    - < ! - - MATH: $r = \frac{ w( x_{i+1}^T, y_{i+1}^T ) } { w( x_{i}, y_{i} -) }$ - - > - -
    5. - -
  5. if r>1, accept the new point - (xi+1,yi+1) = - (xi+1t,yi+1t) - and loop back to step (2)
    6. - -
  6. if r<1, accept the new point with probability - r.  That is, generate R, a random number in - the interval (0,1), and accept - (xi+1t,yi+1t) - as the new point (xi+1,yi+1) if - r > R and loop back to step (2)
    7. - -
  7. otherwise, reject it, and keep the previous point as the - new one (xi+1,yi+1) = - (xi,yi)
    8. - -
  8. loop back to step (2)
    9. -
- -

Koonin and Mederith (Computational Physics, 1989) provide a -simple proof that the algorithm produces the desired random -number distribution w(x).

- -

The step η  in the Metropolis algorithm is -arbitrary.  If chosen small, the acceptance ratio will be -high (a point in the neighborhood of an acceptable point will -also likely be acceptable).  But overall coverage of the -phase space will be very slow.  If η is large, -the phase space will be efficiently covered, but a new point will -seldom be accepted by the Metropolis criteria.  In either -case ineficiency follows. An empirical rule is to choose -η such as to have an acceptance ratio about 0.5 - -0.6.

- -

A concern about the algorithm is that the random numbers -generated by the Metropolis algorithm are correlated.   A -way to guard against this (Koonin and Mederith , Computational -Physics, 1989)  is to calculate the correlation -functions

- -

-< ! - - MATH: $c(k) = \frac{ - ^2 } { - ^2 }$ - - > --"$
- -and select each kth random number -generated by the algorithm with k selected to produce a -small value of Ck.

-
-  Back to top of page -
- -

Guided Monte Carlo -Integral

- -

Lets consider the integral
-  -"$
-The integral is peaked (torus shape) around the -origin.  
-
-The Gaussian function

- -

< ! - - MATH: $w(x,y) = \frac{1}{\pi} e^{ - (x^2+y^2) }$ - - > --"$
-can be chosen to guide the random number selection.

- -

The program Metropolis_Gaussian_Demo.c -  illustrates a random numbers generation according to -w(x,y). Note that arbitrary values for the step size -(η=0.2) and k (k=1) were chosen; no effort was -made to chose them optimally.

- -

This code, with 2000 points, produced the following random -number distribution.

- -


-The program Metropolis_Gaussian_Integral.c -  performs the integral.

- -

Metropolis_Gaussian_Integral -p > data_file
-    or
-Metropolis_Gaussian_Integral

- -

This code prompts the user for more points to be added to the -Monte Carlo integral. The calculation of the error may show the -calculation not to be accurate given that the choices of -η and k were not optimized.

-
-  Back to top of page -
- -

Parallel Monte Carlo -Algorithm

- -

Random Numbers Generation

- -

Pseudo random numbers are often generated by a linear -congruent algorithm.  The integer random number -xi+1 is generated from the previous one -by
-  xi+1 = ( -xi a + b ) mod -c  
-The coefficients a, b and c are chosen to -generate statistically good random numbers in a specific integer -range (Numerical Recipes).  Often the random numbers are -generated using different linear congruence relation, and mixed -via random selection as in RAND1 and RAND2 (Numerical -Recipes).

- -

The procedure to generate random numbers is -deterministic.  Two different nodes on a parallel machine -will produce exactly the same sequence if "seeded" in the same -way. The question then arises as to how to generate "good" random -numbers on the different nodes of the parallel computer that will -remain meaningful when used as a single set.

- -

One solution is to generate sequences on each node based on -different seeds.  The code random.demo_2_series.c -simulates this approach by generating 1 or 2 sequences of random -numbers that can be analyzed for their goodness.  The idea -is to look for the goodness of a sequence in which the first and -second half of the series are generated with different seed -versus sequences generated from a single seed.

- -

random_demo_2_series 2 1000 > data_file

- -

The following filters can be used to acquire some feelings for -the goodness of these random numbers

- - - -

Clearly the random numbers sequences generated from sequences -using different seeds have different properties than sequences -generated from single seed. 

- -

We use here the approach that Random numbers should be -generated on one node and distributed to the nodes which need -them.  This avoids the use of sequences generated from -different seeds.  Realistic applications often require -time-consuming function evaluations. Therefore the inefficiency -inherent in our approach is dwarfed by the large function -evaluation times.  This is also the approach used by  -Gropp et al. in "Using MPI".

-
-  Back to top of page -
- -

Parallel -Monte Carlo Codes

- -

A parallel implementation of the Monte Carlo integration must -have the random sequence distributed to the different nodes, -either by communication or by a local generation of the -sequence.  We implement here three versions of the code.

- -

Master-Workers paradigms

- -

This code generates the random numbers on nodes 0 and passes -them onto the worker nodes.  The latter compute the Monte -Carlo sub-sums and pass these sub-sums to node 0 which finishes -the sums.  The MPI_Reduce() function can be used to expedite -the global sums

- -

- -

Metropolis_Gaussian_Parallel.c -implements this approach.

- -

Client-Server paradigm

- -

This code implements the concept of a random number server; -this is a node reserved to serve random numbers with weight -w(x,y) to whichever node requests it.  The code sets node 0 -as a master node to direct the calculation and interface with the -user.  Node (numprocs-1) is the random numbers -server, and nodes 1 to (numprocs-2) are the worker -nodes which complete the Monte Carlo sub-sums.

- -

The code Metropolis_CS.c implements a version -of of this approach in which the functionality of the different -nodes is relegated to C functions. This illustrates the fact that -variables are local to the nodes. - -

Gropp et al, in "Using MPI" implements this version of a Monte -Carlo integration to calculate π.

- -

Full Workers -version

- -

In this version, node 0 administers the ranges of random -numbers to be generated for each workers, but leaves the actual -generation of the random numbers to the nodes themselves.  -Upon receiving a range, i.e., random numbers range 1001 to 2000, -the worker nodes generate these numbers and calculate the Monte -Carlo sub-sums.  This implies that they must generate -and discard those random numbers prior to their range -selection.  This is a waste of computing cycles, but given -the speed of computing versus the network speed it maybe the most -efficient approach after all.

- -

The scheme is therefore

- - - -

It is left to the reader to implement this approach.