# Thread: 2-d monte carlo integration help

1. ## 2-d monte carlo integration help

First, I'd like to thank everyone that willingly gives their time to help, thank you!!

*NOTE: This is for school, and I am NOT wanting someone to just do this for me or give answers. I just need some tutoring and help because I am completely lost and have made 0 to little progress (I'll post what have done at the end.)

Right now I'm trying to write a program that makes a rectangle, and inside the rectangle is a circle. Virtual darts are thrown via the rand function, and the area of the circle is calculated based on hit / misses. Here are some details:

-I need to compute the area of a shape (circle) by throwing virtual darts
-Produce a tester that tells if my code is working

At the end of the program the following should be printed out:

1. Length and width of enclosing rectangle
2. Area of enclosing rectangle
3. Number of darts thrown
4. How many darts hit the shape
5. Percentage of darts that hit the shape
6. Area of the shape computed by out dart throwing
7. Actual area of shape as known from an equation
8. The absolute value difference between these quantities (i.e., the error)

In the tester we must provide the function: void GetBorders (double *xmin, double *xmax, double *ymin, double *ymax);

I'm particularly confused here, am I suppose to make a GetBorders function, then call the function in the tester function?

Here's some pseudo code that was provided:

TESTER:
(1) Write function for shape we know answer to
--> circle has area PI*r*r
--> Easy to write function for circle wt center on axis!
--> Answer 1 if distance of (x,y) from axis <= r
--> distance: sqrt(x*x + y*y) (pythagorean theorem):
** draw grid wt triangle, to show how hypotenuse is the distance
--> RETURN 1 if sqrt(x*x + y*y) <= r, ELSE 0
==> Call this function:
int DartHits(double x, double y)
(2) Have known answer of area as formula
-> circle is PI*R*R
-> make PI and R cpp macros so we can vary easily
(3) Compare known area to area computed by throwing darts!

And for the other functon:

ThrowDarts:

(1) Get how many darts to throw
--> N, read in from keyboard, store as unsigned long int
(2) Get coordinates of surrounding rectangle
--> xmin: x coordinate of left side of rectangle, which is known to be
less than smallest X value of function
--> xmax: x coordinate of right side of rectangle, >= largest X val of func
--> ymin,ymax: same for Y lower/upper of rectangle
==> types must be double, since points are continuous!
==> double xmin, xmax, ymin, ymax
==> These values must be supplied by tester, since they are tied to the
function being evaluated!
==> we will have tester supply all four values using function:
--> void GetBorders(double *xmin, double *xmax,
double *ymin, double *ymax);
(3) Compute area of rectangle:
--> Area of rect = length * width = (xmax - xmin) * (ymax - ymin)
(4) Loop over N darts thrown, and count hits
a. Generate random # between [xmin,xmax]; call this xdart (xcoord of dart)
-> call xdart = RandNumInRange(xmin, xmax);
b. Generate random # between [ymin,ymax]; call this ydart (ycoord of dart)
-> call ydart = RandNumInRange(ymin, ymax);
c. hit = DartHits(xdart, ydart);
e. if hit != 0, add one to # of hits
(5) Compute area as: rectArea*(percent hits) = rectArea*(nhits/N)
--> Make sure nhits/N computed as a double to avoid integer truncation!
** How do we get a number between xmin, xmax (more generally between w & z)
-> man -a rand brings up manpage
-> set seed using srand
-> rand(void) returns a number between [0, RAND_MAX]
-> Can change integer between [0,RAND_MAX] to real between [0,1]:
- interpret rand() and RAND_MAX as doubles, and simply divide them:
* rn = ((double)rand())/((double)RAND_MAX); /* real between [0,1] */
* make span any distance by multiplying by distance
- Distance is z-w (z known to be >= to w)
* Make start at particular value by adding starting value in
(will change ending value by same amount, so be prepared for that)
==> Can write func:
double RandNumInRange(double W, double Z)

And for what I have so far, I started the Throw Darts bit first, because I am understanding that part slightly more.

Code:
```#include <stdio.h>
#define PI 3.14159

int main (void)

{

int n; /* How many darts to throw. */

printf(" Please enter how man virtual darts you would like to throw: \n\n ");

scanf("%d", &n);

/* Need to call GetBorders function here*/ // Why are they pointers?

int rectArea = (xmax - xmin) * (ymax - ymin);

while ( 1 <= n ) //Started the loop, but I'm not sure what it actually is looping because the instructions seemed arbitrary for making this loop

{

return 0;

}```

I'm so lost at the moment I'm not even sure what else I should ask, but any advice / tips / tutoring is GREATLY appreciated. Again this is homework for me, please do not post answers. Thanks everyone!!!

2. What an interesting assignment!

I wouldn't start with code with this. A whiteboard or paper and pen would be better by far. Until you can visualize at least a portion of what your program needs to do, don't code a darn thing.

There will be a rectangle shown on the display. Your program will randomly hit a location, somewhere inside that rectangle. Some darts will hit the invisible circle inside the rectangle, others will miss the circle. Your program will receive the x and y points of the dart strikes.

Based on numerous darts hits and misses, your program will begin to learn about the position and radius (size) of the invisible circle. Of course, your program will have to know exactly where the circle is located, and what it's size is, because otherwise, it wouldn't be able to tell if the dart throw was a hit on the circle, or a miss. But for the sake of the simulation, your program should not use that information. Instead it should solve the location and size of the circle, with it's own logic, using the dart's data.

There's a little bit of geometry involved, but nothing to cringe about.

Work out how this simulation might go, by paper and pen FIRST. Throw several virtual darts, and by recording where they landed, see if they landed inside an invisible circle inside the rectangle of the paper.

If they DID hit the circle, then that spot MUST be less than or equal distance to the hypotenuse back to the center of the circle. (the radius of the circle, from the center of the circle, to the point the dart landed on).

After awhile, you'll see the logic of parts of the simulation you need. That will form the basic part of the logic for your program.

Give that a try with paper and pen, and work it through by hand first. Note whether the circle's center can be anywhere inside the rectangle, or is it always centered in the rectangle?

P.S.
Welcome to the forum, jsuite!

3. Hey, thank you so much for your time and help! I'm going to take your advice. One I finish working it through with pen and paper and by hand I'll update the code, and ask questions from there if need be. Thanks also for the warm welcome, please try to stay updated with this thread! Again, thank you!!

4. Start your update, with a list of the required functions, (just empty of code for now), from your assignment. It has very specific requirements - what variables go where, etc.

The number of dart throws should be stored in the ThrowDarts function, and be labeled N (not in main() as n). All the required stuff listed above, can then be put into the right functions, accordingly.

5. Okay, sweet! I'm slowly working on it (studying for other things aswell, very loaded this weekend) so expect an update sometime this evening. But there is one question I have, what do I do with main? That's something that was really confusing me, it never mentions main... I know for example, that if a variable need not be static, then I can put it in there... So one idea I have id that my main should have world variables, and also the final print out statements?

6. Originally Posted by jsuite
Okay, sweet! I'm slowly working on it (studying for other things aswell, very loaded this weekend) so expect an update sometime this evening. But there is one question I have, what do I do with main? That's something that was really confusing me, it never mentions main... I know for example, that if a variable need not be static, then I can put it in there... So one idea I have id that my main should have world variables, and also the final print out statements?
Programs in C usually have 3 functional parts - Input, Calculation/s, and Output. In a small program, all three can fit nicely inside main(), but on larger programs, you may need one or more functions, for each of these parts of the program.

Try thinking of main() as the starting point for your program, and the subsequent point of return. It could be just two lines of code, or several blocks of code.

Since several other functions have specific requirements, I'd get them set up first - even if the requirements are not coded, but just listed in notes:
Code:
```//?? ThrowDarts(unknown) {
/*    gets unsigned long N for number of darts to be thrown from the user */ //kind of thing.

//return??
//}```
Save the contents of main() and other non-required (by the assignment) functions, until later. You can be much more flexible with them. Use them for "glue", to make the logic flow OK.

In general, list the requirements, inside the required functions (if known), but put off all the details that are not important, until you get the core logic flowing.

Once the core logic is working, then add in the details. Saves a lot of time.

Global variables (which have program-wide scope), are listed above main(). They can be useful, but are a negative on the whole purpose of functions (which is to encapsulate a portion of the program). Use them VERY sparingly, if at all.

7. Remember not to rely on rand() or random() family of functions in C, as they are linear congruential pseudorandom number generators (Wikipedia), and are likely to introduce artefacts/errors when used in Monte Carlo stuff. I warmly recommend using Xorshift instead. In particular:
Code:
```#include <stdint.h>

static uint32_t prng_state[4] = {
123456789U,
362436069U,
521288629U,
88675123U
};

static inline uint32_t prng_u32(void)
{
uint32_t        temp;

temp = prng_state[0] ^ (prng_state[0] << 11U);
prng_state[0] = prng_state[1];
prng_state[1] = prng_state[2];
prng_state[2] = prng_state[3];
return prng_state[3] ^= (temp >> 8U) ^ temp ^ (prng_state[3] >> 19U);
}

/* 0.0 <= prng_one() < 1.0 */
static inline double prng_one(void)
{
return (double)prng_u32() / 18446744073709551616.0
+ (double)prng_u32() / 4294967296.0;
}```
To randomize the generator, set the four 32-bit unsigned integers in prng_state[] to anything (except all four to zero). In prng_one() the divisors are powers of two, so on most architectures they are either optimized away, or implemented as very fast multiplications or bit shifts. It is seriously fast, and produces excellent pseudorandom numbers.

Although Adak gave you excellent advice, I think you can make it even easier for yourself, if you think in terms of centerpoints.

Compute the centerpoint for the rectangle and the circle, the width and height of the rectangle, and the squared radius of the circle:
Code:
```    /* Assuming xmin, xmax, ymin, ymax are the rectangle,
* circle has radius radius, and the circle is centered on the rectangle. */
const double   rect_width = xmax - xmin;
const double   rect_height = ymax - ymin;
const double   rect_x = xmin;   /* Minimum x for the rectangle */
const double   rect_y = ymin;   /* Minimum y for the rectangle */
const double   circle_x = 0.5 * (xmin + xmax); /* Center at center of rectangle! */
const double   circle_y = 0.5 * (ymin + ymax); /* Center at center of rectangle! */
To generate a random point within the rectangle, you generate a point in the unit square (0,0)-(1,1), scaled and moved to cover the desired rectangle:
Code:
```    x = rect_x + rect_width * prng_one();
y = rect_y + rect_height * prng_one();```
The distance squared from point (x,y) to the center of the circle is, as Adak explained, described by the Pythagorean theorem, where the length of the two sides are x-circle_x and y-circle_y. In other words,
Code:
```    if ( (x-circle_x)*(x-circle_x) + (y - circle_y)*(y - circle_y) <= circle_squared ) {
/* Point (x,y) is within the circle.
* (If == circle_squared, then the point is on the circumference of the circle.) */
}```
For fun, I wrote a test program which takes all the parameters from the command line, checks them, and does the 2D Monte Carlo integration using the above Xorshift implementation as random number source. It works fine (as in yields the expected results at expected precision), and is just over a hundred lines of code (not counting empty lines).

8. Very interesting note about circles inside squares about half way down the page, with a pic (applet):
Approximations of

9. Okay, so I worked on it some, tried to do what I could. I must say I'm embarrassed at how little I can do thus far. I really got screwed with this course because my curriculum calls for us to start on java (which i did) then it switches us to c (what I'm in now). Long story short I chose a professor who I knew was reputable, but the changed my professor dude to an influx of computer science students this year, so they had to add and shuffle around professors. Well I get stuck with a young woman who is not very good at teaching. The problem is, the head of the department is very strict and tells her what to teach, but she never gets to everything, so when I get a project (The project is made by the head of the department) it's vey challenging, because now I have to go back and teach myself everything that she failed to even bring up. Anyways, here's what I have to far, I'm going to keep chiseling away at it.

Code:
```
#include <stdio.h>
#define PI 3.14159

int main (void)

{

//stuff

}

tester() /* tester function */
{

double distance = sqrt(x*x + y*y); /* This is the Pythagorean that i'm suppose to use, where does the x and y come from?*/
/* The circle is suppose to be centered at (0,0) so if I've chosen my radius to be 10, i can calculat x and y? */

//Here I'm suppose to draw grid with triangle, to show how hypotenuse is distance.. How do I draw it?

return 1;

else
return 0;

// call the function int DartHits(double x, double y); here
// This functions is how we know the estimated area of circle

//Now we compare actual area vs. estmated area. I think the mean print out statement?
//like this:

printf("The actual area is %f, and the estimated area is %d /n", cirAreal, DartHits);

}

ThrowDarts() /* ThrowDarts function */

{

int n;
int i;

printf("How many darts would you like to throw?");
scanf("%d", &n); /* Prompts user for how many darts to throw, and saves it n */

/* completely lost on the next step: I'm suppose to call the function void
GetBorders(double *xmin, double *xmax, double *ymin, double *ymax); and it says that
all four of these values are supplied by the tester function. So is Getborders
seperate function, or is it inside the tester function? Very confused here. */

int rectArea = (xmax - xmin) * (ymax - ymin);

for (i = 1; i <= n, i++)
// The following should go under this for loop:

/* a. Generate random # between [xmin,xmax]; call this xdart (xcoord of dart)
-> call xdart = RandNumInRange(xmin, xmax);
b. Generate random # between [ymin,ymax]; call this ydart (ycoord of dart)
-> call ydart = RandNumInRange(ymin, ymax);
c. hit = DartHits(xdart, ydart);
e. if hit != 0, add one to # of hits */

// One thing that confuses me is, in part a and b above
// It says to take a variable like xdart and ydart, then
// make them = to functions? How does that work?
// also those instructions say get a random number between xmin and xmax
// and call that number x dart
// then on the next like it says make xdart = a function???

double  estRectArea = rectArea * (nhits/n); //nhits is number of hits, generated by for loop above

/* and the last part of this function is suppose to be how they want us to calculate rand

/* ** How do we get a number between xmin, xmax (more generally between w & z)
-> man -a rand brings up manpage
-> set seed using srand
-> rand(void) returns a number between [0, RAND_MAX]
-> Can change integer between [0,RAND_MAX] to real between [0,1]:
- interpret rand() and RAND_MAX as doubles, and simply divide them:
* rn = ((double)rand())/((double)RAND_MAX);  real between [0,1]
* make span any distance by multiplying by distance
- Distance is z-w (z known to be >= to w)
* Make start at particular value by adding starting value in
(will change ending value by same amount, so be prepared for that)
==> Can write func:
double RandNumInRange(double W, double Z) */

}

return 0;

}```

10. A lot of questions I have are in the code comment.

11. Also, the circle is centered at (0,0).

12. Originally Posted by jsuite
Also, the circle is centered at (0,0).
Good to know.

In ThrowDarts(), the number of darts variable should be N, (uppercase, no idea of why), and the data type should be an unsigned long int, instead of a regular int.

GetBorders() will get from the user, where the top, bottom, left and right borders for the rectangle will be. It should be called only once when the program starts, which matches up with getting the number of darts that will be thrown. I'd put this prototype right inside ThrowDarts, itself, and call it from inside ThrowDarts.

double xmin, ymin,xmax,ymax;
xmin=1;ymin=1;
xmax=79;ymax=30 to 40;

If you write to the 80th column of the display, the cursor will automatically do a line feed, and wreck your display badly. I believe we want the max-min of both x and y, to be odd numbers. Otherwise you won't have a true center for the rectangle.

_gotoxy(xmin,ymin); printf("*"); _gotoxy(xmax,ymax); printf("*");getchar();

The above is what I can see in a console on Windows, depending on the properties and fonts I have set up for that window.

Try it and see what you can see (upper left and lower right corners only), without the screen rolling or the cursor adding a newline.

You have conio.h?

It's best if we zero in on ThrowDarts function first, and include GetBorders function since it's required for ThrowDarts to work.

Best if you make comments OUTSIDE the code. Otherwise my comments back will be inside the code, and that may look to someone like I'm writing your program.

13. Ok to address your post above: Thank you so much for the response! About the GetBorders functions and the prototype, well he wants us to use void GetBorders (double *xmin, double *xmax, double *ymin, double *ymax). I find this absurd BUT the point is to practice with pointers, it was a challenge for me t figure out how to properly call it. Then it dawned on me, this is similar to a swap (*pointer, *pointer) function.

Also, I don't think I'm actually suppose to draw the shapes? I can check about using conio., but graphics.h did nt work for me. I'm using vim editor in ubuntu 12.01 btw.

Good news though, I have made significant progress. I have a few huge questions Ill ask after I put the code here:

Code:

Code:
```#include <stdio.h>
#define PI 3.14159

void GetBorders(double *xmin, double *xmax, double *ymin, double *ymax)
{

double a, b, c, d;

a = 27;
b = 57;
c = 33;
d = 67;

*xmin = a;
*xmax = b;
*ymin = c;
*ymax = d;
}

Tester()
{
double xmax, xmin, ymax, ymin;

GetBorders(&xmin, &xmax, &ymin, &ymax);
}

ThrowDarts()
{
double n;
printf("How many darts would you like to throw?");
scanf("&d", &llu);

double xmax, xmin, ymax, ymin;
GetBorders(&xmin, &xmax, &ymin, &ymax);
double rectArea = (xmax -xmin) * (ymax - ymin);

double xdart;
double ydart;
int i = 1;

//loop goes here:
while (i <= n)
double nhit;

xdart = rand(xmin,xmax); //generate random number between xmin and xmax
xdart = RandNumInRange(xmin,xmax);

ydart = rand(ymin,ymax); //generate random number between ymin and ymax
ydart = RandNumInRange(ymin,ymax);

nhit = DartHits(xdart, ydart);
if (nhit != 0)
nhit++;

double estRectArea = rectArea * (nhit/n);
}

double RandNumInRange (double w, double z)
{

}

int main (void)
{
}

```

Questions:

The RandNumInRange - I don't know where to start. It needs to pick a random number in a range, and use a seed. Tried googling about seeds, but they are all explained relative to time.h. I'm not suppose to use time.h. So what are seeds, and can you explain what I should be doing when making this function?

*Note: things are all out of order, I know. Once I get the RanNumInRange function figured out, ill be fine tuning it, which I'm sure will lead me to more questions. Also the reason why I was super super lost at the very beginning is that I was missing out on some super super fundamental rules of C. (lol) I'll spare you the details though.

14. The random number generator in C is designed to be able to repeat the sequence at will (using the same compiled program in the same OS using the same libraries). This means that you can reproduce the data anytime you want.
Basically it does this by using a number for each sequence.
you have sequence 1: it starts with 674532, then 5421, then 9811111, ..., ... (numbers made up)
you have sequence 2: it starts with 3145, then 3312967, then 7777, ..., ... (numbers made up)

If you don't specify a sequence for the random number generator it uses sequence 1.
You can basically pretend to specify a random sequence by using the number of seconds since sometime in 1970 (the return value of time(NULL)) by using srand(time(NULL)).

srand() instructs the random number generator to use the numbered sequence; rand() retrieves the next number from the "active" sequence.

Basically you use srand() once for the whole application (whether with a predefined sequence number or a "random" one is irrelevant) and rand() as often as needed.

15. This program shows a couple ways to get random numbers. The first one gives a double variable, an random number in the range you want. The double always ends with .0, which is to say, it's always a whole number.

The second way it shows to get a random number, is to give a real double, in the range of 0-1. That can obviously be boosted, but it's one of the requirements for your program.

Code:
```#include <stdio.h>
#include <stdlib.h>

int main(void) {
double xmin, xmax, ymin, ymax,rnumd;
int i, rnum,xmaxi;

srand(1);

for(i=0,xmax=30;i<300;i++) {
/* gives even doubles between 1 and xmax */
//rnumd = (double) (rand() % (int) (xmax +1));   //run me too, but comment

/* gives random doubles between 0 and 1 */       //this line out when you run it.
rnumd = ((double)rand()/(double)RAND_MAX);

printf("%5.2f   ",rnumd);
}
printf("\n");
return 0;
}```
srand(someInteger), "seeds" the random number generator with a string of digits. Each seed number gives a different string of random numbers.

If you give the same seed number (or no seed number), or you don't call srand() before calling rand(), your string of random numbers will always be the same (very good for testing purposes). You should call srand() only ONCE each time the program starts. NEVER repeat the call, thinking you'll get a more random number - you won't.

rand() is not a great random number generator, but it's MORE than what you need for this program.