Code Review Asked by Kartik Chhajed on December 8, 2020
This is my first C program. I want to optimize this simulation.
Algorithm
The simulation algorithm is:
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 %lfnn", avg_density);
return 0;
}
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
.
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;
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.
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 %lfn", monte_carlo(BETA));
}
Correct answer by G. Sliepen on December 8, 2020
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP