Thread: 2-d monte carlo integration help

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!

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! */
    const double   circle_squared = radius * radius;

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).

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
#define RADIUS 10


int main (void)


{


//stuff


}






tester() /* tester function */
{


double  cirArea = PI * RADIUS * RADIUS;
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?


if (distance <= RADIUS)
   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;


}

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.

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.

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.