Monte carlo simulation in C

Question

This is my first C program. I want to optimize this simulation.

Algorithm
The simulation algorithm is:

1. The system can go from \$i\$ to \$i+1\$ with probability \$e^{-L\theta(\rho_i)}\$, where \$\rho_i=i/L\$, and \$\theta(\rho)=\rho(2b-\rho)\$.
2. When system reaches \$n=0\$, it is revieved (resurected) to position based on how much time it spent on \$n>0\$.
3. At the end we are intrested in knowing \$\langle\rho\rangle=\sum_{t}\rho_t\$.

Code
Following is the code. I believe this code can be compactified also. I think I do not understand the norms of the ANSI C standard. Feel free to correct me anywhere. I also do not understand if I am properly generating random numbers or not!.

#include<stdio.h>
#include<stdlib.h> 
#include<math.h>
#include<time.h> 

#define LATTICE 10
#define SWEEPS 100000000
#define BETA 0.5

float rho[LATTICE];

/*
Following function will make variable defined above:
    rho = 1/10, 2/10, 3/10, ... 1.
*/ 
void initialize()
{
    for(int i = 0; i < LATTICE; i++)
        rho[i] = (i+1.0)/LATTICE;
    return;
}

/*
These are the rates for going form n to n+1 for a given b.
*/
void rates(double *death_rate, float b)
{
    double theta;
    for(int i=0; i < LATTICE; i++)
    {
        theta = rho[i]*(2*b - rho[i]);
        *death_rate = exp(-LATTICE*theta);
        death_rate++;
    }
    return;
}

/*
Following function generates uniform random number
between 0 to 1.
*/
double uniform_rand()
{
    return (double)rand()/(double)RAND_MAX;
}


/*
The following function revive the system when n becomes -1
with the weights = distribution.
*/
int revive(unsigned long long *distribution, unsigned long long norm)
{
    int n = -1;
    double cumsum = 0.0, u_rand = uniform_rand();
    while(cumsum <= u_rand)
    {
        cumsum += (double) *distribution/(double)norm;
        distribution++;
        n++;
    }
    return n;
}

/*
Following function calculate the average density.
*/
double rho_avg(unsigned long long *distribution, unsigned long long norm)
{
    int i;
    double avg_density = 0.0;
    for (i=0; i<LATTICE; i++)
    {
        avg_density += (rho[i]*(*distribution))/norm;
        distribution++;
    }
    return avg_density;
}

double monte_carlo(float b, int n)
{
    unsigned long long i;
    unsigned long long distribution[LATTICE]={0};
    double death_rate[LATTICE];
    srand(time(0));
    rates(death_rate, b);
    for(i = 1; i <= SWEEPS; i++)
    {
        distribution[n] += 1;
        if (uniform_rand() < death_rate[n])
        {
            n--;
            if (n == -1)
                n = revive(distribution, i);
        }
    }
    return rho_avg(distribution, SWEEPS);
}


int main(void)
{ 
    int i;
    double avg_density;
    initialize();
    avg_density = monte_carlo(BETA, LATTICE-1);
    printf("\nAverage density is %lf\n\n", avg_density);
    return 0;
}

Answer 1

Naming things

Naming things is one of the hard things in computer science. You should make sure names of functions and variables are clear and concise. In general, use nouns for variable names and verbs for function names. Looking at your choices of names in more detail:

• LATTICE is not a lattice, it is the size of the lattice. So call it LATTICE_SIZE.
• SWEEPS is the number of sweeps to perform, so maybe it is better to call it N_SWEEPS (N is a commonly used prefix meaning "number of").
• rates() is a function, so use a verb for the function name, for example calculate_rates().
• rho_avg() is a function again, so use a verb for that as well, like calculate_rho_avg().

You should also be consistent in how you name things. Is it rho or density? Pick one and stick with it. I would also write beta instead of b, to match how you write down other greek letters like rho and theta.

Use array indexing notation where appropriate

In rates(), you are using pointer arithmetic when you could just have used standard array notation:

   for(int i = 0; i < LATTICE; i++)
   {
        theta = rho[i] * (2 * b - rho[i]);
        death_rate[i] = exp(-LATTICE * theta);
   }

Similarly, in revive(), write:

    for(n = 0; cumsum <= u_rand; n++)
    {
        cumsum += (double)distribution[n] / (double)norm;
    }

    return n - 1;

Terminology

Death rates? Revive? That sounds very morbid! Unless you are simulating some system that predicts cell death or pandemics, this is not terminology normally used in Monte Carlo simulations. Your code looks like it implements the Metropolis algorithm, where your death_rate would be equivalent to the transition probability, although I'm not sure what the equivalent of revive() would be. If it is the Metropolis algorithm you are implementing, then it doesn't look like you have detailed balance. What system are you simulating exactly? It might help to document that in the code as well.

Avoid global variables

It is good practice to avoid using global variables. That makes it easier if your program grows and becomes part of something larger. Or perhaps you want to run multiple simulations at the same time using threads. This should be easy; you already have the arrays distribution[] and death_rate[] inside monte_carlo(), just move rho[] there as well and pasas a pointer to rho to rates().

You might want to do this in a more organized way, and create a struct to hold all the information relevant to your simulation:

struct simulation_parameters {
    unsigned long long distribution[LATTICE];
    double death_rate[LATTICE];
    double rho[LATTICE];
};

...

double monte_carlo(float beta) {
    struct simulation_parameters params = {0}; // sets everything in params to zero
    calculate_rho(params.rho); // should do what initialize() did
    calculate_death_rates(params.death_rate, beta);

    for (unsinged long long i = 1; i <= SWEEPS; i++) {
        distribution[n]++;
        if (uniform_rand() < params.death_rate[n]) {
            n--;
            if (n == -1)
                n = revive(params.distribution, i);
        }
    }

    return calculate_rho_avg(distribution, SWEEPS);
}

int main(void) {
    srand(time(0)); // srand() should only be called once, so do it in main()
    printf("Average density is %lf\n", monte_carlo(BETA));
}