This page was last updated: March 27, 2024
Monte Carlo simulation is used to determine the range of outcomes for a series of parameters, each of which has a probability distribution showing how likely each option is to happen. In this project, you will take a scenario and develop a Monte Carlo simulation of it, determining how likely a particular output is to happen.
Clearly, this is very parallelizable -- it is the same computation being run on many permutations of the input parameters. You will run this with OpenMP, testing it on different numbers of threads (at least 1, 2, 4, 6, and 8).
This project is inspired by this video. Be sure to watch it before reading on (it's only 41 seconds).
A golf ball is set in motion at the top of a ski jump at height BeforeY.
It rolls until its height is AfterY at which point it shoots out horizontally.
It continues until it hits the ground.
A hole of radius RADIUS exists at a horizontal distance DistX from the end of the ski jump.
Will the golf ball land inthe hole?
If we knew the exact measurements, we could figure this out exactly, but we don't have exact measurments. Several people have taken measurements and all have come up with different numbers. Their mean measurements and tolerances (±) are shown in the table below:
| Quantity | Mean Value | ± Value |
|---|---|---|
| BeforeY | 80.f | 5.f |
| AfterY | 20.f | 1.f |
| DistX | 70.f | 5.f |
Given all of this uncertainty of the dimensions, what is the probability that the golf ball will actually end up in the cup?
OK, you physics hounds, knock it off. Yes, we are going to neglect depth effects, moment of inertia, air resistance, etc, etc. Let's just have fun with this problem as stated.
Run this for some combinations of trials (NUMTRIALS) and threads (NUMT). Do timing for each combination. Like we talked about in the Project Notes, run each experiment some number of times (NUMTIMES) and record just the peak performance.
Produce a rectangular table and two graphs. The two graphs need to be:
Choosing one of the runs (one of the ones with the maximum number of trials would be good), tell me what you think the actual probability is.
Choosing one of the runs, compute Fp, the Parallel Fraction, for this computation.
Given this Parallel Fraction, what would the maximum speedup be if you could throw hundreds of cores at it?
The code below is printed in the handout to make it easy to look at and discuss.
Note: if you are on Windows, take the ", stderr" out of the #pragma line!
#include <stdio.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
// print debugging messages?
#ifndef DEBUG
#define DEBUG false
#endif
// setting the number of threads:
#ifndef NUMT
#define NUMT 2
#endif
// setting the number of trials in the monte carlo simulation:
#ifndef NUMTRIALS
#define NUMTRIALS 50000
#endif
// how many tries to discover the maximum performance:
#define NUMTIMES 20
#define CSV
// ranges for the random numbers:
#define GRAVITY 32.2f
const float BEFOREY = 80.f;
const float AFTERY = 20.f;
const float DISTX = 70.f;
const float RADIUS = 3.f;
const float BEFOREYDY = 5.f;
const float AFTERYDY = 1.f;
const float DISTXDX = 5.f;
float BeforeY[NUMTRIALS];
float AfterY[NUMTRIALS];
float DistX[NUMTRIALS];
float
Ranf( float low, float high )
{
float r = (float) rand(); // 0 - RAND_MAX
float t = r / (float) RAND_MAX; // 0. - 1.
return low + t * ( high - low );
}
// call this if you want to force your program to use
// a different random number sequence every time you run it:
void
TimeOfDaySeed( )
{
struct tm y2k = { 0 };
y2k.tm_hour = 0; y2k.tm_min = 0; y2k.tm_sec = 0;
y2k.tm_year = 100; y2k.tm_mon = 0; y2k.tm_mday = 1;
time_t timer;
time( &timer );
double seconds = difftime( timer, mktime(&y2k) );
unsigned int seed = (unsigned int)( 1000.*seconds ); // milliseconds
srand( seed );
}
int
main( int argc, char *argv[ ] )
{
#ifdef _OPENMP
#ifndef CSV
fprintf( stderr, "OpenMP is supported -- version = %d\n", _OPENMP );
#endif
#else
fprintf( stderr, "No OpenMP support!\n" );
return 1;
#endif
TimeOfDaySeed( ); // seed the random number generator
omp_set_num_threads( NUMT ); // set the number of threads to use in parallelizing the for-loop:`
// better to define these here so that the rand() calls don't get into the thread timing:
// fill the random-value arrays:
for( int n = 0; n < NUMTRIALS; n++ )
{
BeforeY[n] = Ranf( BEFOREY - BEFOREYDY, BEFOREY + BEFOREYDY );
AfterY[n] = Ranf( AFTERY - AFTERYDY, AFTERY + AFTERYDY );
DistX[n] = Ranf( DISTX - DISTXDX, DISTX + DISTXDX );
}
// get ready to record the maximum performance and the probability:
double maxPerformance = 0.; // must be declared outside the NUMTIMES loop
int numSuccesses; // must be declared outside the NUMTIMES loop
// looking for the maximum performance:
for( int times = 0; times < NUMTIMES; times++ )
{
double time0 = omp_get_wtime( );
numSuccesses = 0;
#pragma omp parallel for default(none) shared(BeforeY, AfterY, DistX, stderr) reduction(+:numSuccesses)
for( int n = 0; n < NUMTRIALS; n++ )
{
// randomize everything:
float beforey = BeforeY[n];
float aftery = AfterY[n];
float distx = DistX[n];
float vx = ?????
float t = ?????
float dx = ?????
if( ????? <= RADIUS )
numSuccesses++;
} // for( # of monte carlo trials )
double time1 = omp_get_wtime( );
double megaTrialsPerSecond = (double)NUMTRIALS / ( time1 - time0 ) / 1000000.;
if( megaTrialsPerSecond > maxPerformance )
maxPerformance = megaTrialsPerSecond;
} // for ( # of timing tries )
float probability = (float)numSuccesses/(float)( NUMTRIALS ); // just get for last NUMTIMES run
#ifdef CSV
fprintf(stderr, "%2d , %8d , %6.2f , %6.2lf\n",
NUMT, NUMTRIALS, 100.*probability, maxPerformance);
#else
fprintf(stderr, "%2d threads : %8d trials ; probability = %6.2f ; megatrials/sec = %6.2lf\n",
NUMT, NUMTRIALS, 100.*probability, maxPerformance);
#endif
}
Print out: (1) the number of threads, (2) the number of trials, (3) the probability of a golf ball going into the cup and (4) the MegaTrialsPerSecond. Printing this as a single line with commas between the numbers but no text is nice so that you can import these lines right into Excel as a CSV file. That's what the CSV #define is for.
To choose a random number between two floats, use:
#include <stdlib.h>
float
Ranf( float low, float high )
{
float r = (float) rand(); // 0 - RAND_MAX
float t = r / (float) RAND_MAX; // 0. - 1.
return low + t * ( high - low );
}
// call this if you want to force your program to use
// a different random number sequence every time you run it:
void
TimeOfDaySeed( )
{
struct tm y2k = { 0 };
y2k.tm_hour = 0; y2k.tm_min = 0; y2k.tm_sec = 0;
y2k.tm_year = 100; y2k.tm_mon = 0; y2k.tm_mday = 1;
time_t timer;
time( &timer );
double seconds = difftime( timer, mktime(&y2k) );
unsigned int seed = (unsigned int)( 1000.*seconds ); // milliseconds
srand( seed );
}
Put the following lines into a file called Makefile:
proj01: proj01.cpp
g++ proj01.cpp -o proj01 -lm -fopenmp
Run it as:
make proj01
./proj01
You can save yourself a ton of time by setting this up to run from a script.
Check the Project Notes to see how to do that in bash, C-shell, or Python.
If you've never done this before, learn it now!
You will be surprised how much time this saves you throughout this class, and other classes.
Here it is as a bash script:
#!/bin/bash
for t in 1 2 4 6 8
do
for n in 1 10 100 1000 10000 100000 500000 1000000
do
g++ proj01.cpp -DNUMT=$t -DNUMTRIALS=$n -o proj01 -lm -fopenmp
./proj01
done
done
Run it as:
bash loop.bash >& proj01.csv
Diverting it to a CSV file sets you up to import it right into Excel.
Your turnin will be done at http://teach.engr.oregonstate.edu and will consist of:
| Feature | Points |
|---|---|
| Proper rectangular data table | 10 |
| Good graph of performance vs. number of trials | 20 |
| Good graph of performance vs. number of threads | 20 |
| Good estimate of the probability | 15 |
| Good estimate of Fp, the Parallel Fraction (show your work) | 15 |
| Commentary | 20 |
| Potential Total | 100 |