Wednesday, July 13, 2016

C - Monte Carlo Method

In the following article, I'll create an application which will obtain the volume of a cylinder with the Monte Carlo method (https://en.wikipedia.org/wiki/Monte_Carlo_method). I'm not going to give a detailed explanation of the Monte Carlo method, but just how to use it in the current application. Basically, the Monte Carlo method use a random generator which gives correct and incorrect values of our problem. The amount of correct results over the incorrect will give us an approximative value our problem. In the current example, I'll first create an integer random generator, which will give points inside cube and we must check then if the points are in the cylinder or not. As this is an approximative method, we must do the test several times to have a more accurate result, thus having more samples, the result will be more accurate too.

The application can be divided in 4 files:
  • The main file (main.c): which will handle all the main methods of the following files;
  • The volume calculator (volume.c): which calculate the volumes iterating over random numbers;
  • The integer random generator file (lcg.c): it uses the LCG (Linear Congurential Generator) method to generate random numbers;
  • The plot file (plot.c) which will plot all the results to show how the each result precision increases as the samples increase.

There is no single solution for the given exercise, and there is no single way to handle it, in the current article, I'll start with the more generic method and go into the more specific (as specified in the previous list).

The main.c will use the public method of the volume calculator and the plot files. Let's assume that the volume calculator public method has the following input arguments: the size of the block (let's also assume that the block is the smallest possible to stock a cylinder inside) with the X, Y, and Z values; and number of samples with which the calcule must be made. The size of the block can be obtained by the input parameter of the main application, or be defined in our main application as constants (I choose this option). So, I want five results with the same number of samples, and 6 different number of samples (100, 10.000, 100.000, 1.000.000, 10.000.000 and 100.000.000). This tell us that two loops are needed, one inside the other. At every loop, we must store the result into a list of results (to be plotted later). And when finished the loops, let's use the plot method. The plot method will need the results of the volume calculation and the number of samples used for each result (this made a two dimensional graph), and the size of our list. Here is the code of main.c:


The volume.c will just retrieve random points, check if their inside or outside the cylinder, and calculate the volume of the cylinder from the probability a point is inside/outside the cylinder. The number of samples is given as input parameter (it was assumed in the main.c file). With this, we know that we need a public method of the random generator, a loop to iterate over the number of samples, and a private method to know if a point is inside or outside the cylinder. There is no private method in C, I only use the word "private" to tell that this method will only be used inside the volume.c file. For the loop, we only have to iterate from 0 to N (the number of samples), at each iteration, generate a single random point, and if the point is inside the cylinder, increment a counter. To generate a random point, just generate a random integer from 0 to the coordinate X, Y or Z (if we assume the integer random generator gives us a number from 0 to 1, we only have to multiply it with the large of the block). And to check if a point is inside or outside the cylinder, just check if the X and Y coordinates are inside an ellipse (the Z coordinate will always be if the first condition is true). Here is the code of volume.c:


Let's continue with the lcg.c file. This file only contains single method described before, which generates an integer between 0 and 1. To generate a random integer, I've used the Linear Congruential Generator with the values of Knuth & Lewis. And as a seed, I've used the current timestamp limited by the congruential value. Here is the code:


And finally the plot.c file. This file will only generate a graph (with gnuplot) from some given points. To plot all the points, I've set the limits of the graph 10% greater than the limits (to avoid having the maximum and minimum values in the corners). Having those limits, I've just used the gnuplot application to generate a graph. As the gnuplot application is in the system (and not an available library in C) I've first opened a stream to communicate with the application, and then send all the commands. Here is the result:


As said in the beginning, this is not the unique solution (nor the best), and the way I solve it (from the generic to the specific) isn't the unique way to do so. Finally, having all the code, to compile it and run it, I've done it with a terminal and the gcc command.


The -Wall arguments means Warning All, print all the warning message when compiling. And the -lm means to search for system library when linking.

No comments:

Post a Comment