Monte Carlo virus infection simulation

Question

I am learning C++ from a C background, and I was wondering if I could convert this old C program I'd written into C++ as a personal exercise. The original C program is a simple Monte Carlo simulation that simulates a virus spreading throughout a school population.

Currently, my C++ implementation is almost identical to the C version, despite a few minor changes. Given this is a somewhat large block of algorithmically-focused code, I am not looking for specific details, but I was wondering if any of these C constructs could be replaced with more idiomatic C++ versions.

For the record, I am aware of the minor differences, such as the possibility of using iostreams instead of printf for the I/O, but I often find C's IO to be more expressive in this situation, anyway. I'm more curious about more general problems.

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <climits>
#include <ctime>

static long NUM_SIMULATIONS;
static int INFECTION_THRESHOLD;

static int uninfected_students = 0;
static int infected_students[] = {0, 0, 0, 0, 0};

void simulate_day();
int run_simulation();

int main(int argc, char** argv) {
    srand(time(NULL));

    if (argc < 2 || argc > 3) {
        printf("usage: maine <rate> [simulations=1000]\n");
        exit(EXIT_SUCCESS);
    }

    double transmission_rate = strtod(argv[1], NULL);
    if (transmission_rate <= 0 || transmission_rate > 1) {
        fprintf(stderr, "maine: error: rate must be a number in the range (0, 1]\n");
        exit(EXIT_FAILURE);
    }
    INFECTION_THRESHOLD = (int) (RAND_MAX * 0.04 * transmission_rate);

    if (argc > 2) {
        NUM_SIMULATIONS = strtol(argv[2], NULL, 10);
        if (NUM_SIMULATIONS <= 0) {
            fprintf(stderr, "maine: error: simulations must be a number in the range (0, %li]\n", LONG_MAX);
            exit(EXIT_FAILURE);
        }
    } else {
        NUM_SIMULATIONS = 1000;
    }

    double infected_students = 0;
    int successful_infections = 0;
    printf("Running simulation..."); fflush(stdout);

    clock_t start_time = clock();

    for (int i = 0; i < NUM_SIMULATIONS; i++) {
        int simulation = run_simulation();
        infected_students += simulation;
        if (simulation >= 10)
            successful_infections++;
        printf("\rRunning simulation... (%06.2f%%)", (double) i / NUM_SIMULATIONS * 100.0);
        fflush(stdout);
    }

    clock_t end_time = clock();
    clock_t elapsed_time = end_time - start_time;

    //        Running simulation... (000.00%)
    printf("\rSimulation complete.  (%6.2fs)\n", (double) elapsed_time / CLOCKS_PER_SEC);
    infected_students /= NUM_SIMULATIONS;

    printf("\n");
    printf("    Average number of students infected:    %8.4f (%05.2f%%)\n", infected_students, infected_students / 300.0 * 100.0);
    printf("      Infection probability per student:    %15.2f%%\n", (infected_students - 1) / 299.0 * 100.0);
    printf("      Number of 'successful infections':    %8i (%05.2f%%)\n", successful_infections, (double) successful_infections / NUM_SIMULATIONS * 100.0);
    printf("\n");

    return 0;
}

int run_simulation() {
    uninfected_students = 299;
    infected_students[0] = 1;
    memset(infected_students + 1, 0, sizeof(*infected_students) * 4);

    while (uninfected_students > 0 && (infected_students[0] > 0 || infected_students[1] > 0 || infected_students[2] > 0 || infected_students[3] > 0 || infected_students[4] > 0)) {
        simulate_day();
    }

    return 300 - (uninfected_students < 0 ? 0 : uninfected_students);
}

void simulate_day() {
    int contagious_students = infected_students[3] + infected_students[4];
    int newly_infected_students = 0;
    for (int i = 0; i < contagious_students; i++) {
        int current_uninfected = uninfected_students;
        for (int j = 0; j < current_uninfected; j++) {
            if (rand() >= INFECTION_THRESHOLD) continue;
            uninfected_students--;
            newly_infected_students++;
        }
    }

    memmove(infected_students + 1, infected_students, sizeof(*infected_students) * 4);
    infected_students[0] = newly_infected_students;
}

Answer 1

@Nobody has already commented on a number of issues at the level of individual lines of code. While these changes would help the code make a little better use of C++, most of them apply about equally to C.

To make noticeably better use of C++, you probably want to look at the basic structure of the program. It seems reasonable to me that a "simulation" would be an object. That class would then include member functions to execute a simulation, and to execute a single day of the simulation. Initialization would (of course) happen in the constructor.

The other obvious point would be to make better use of the standard library. The standard collections and algorithms (for two obvious points) can contribute quite a bit to making the code simpler and more readable.

Personally, I'd start with a simple range class that lets me iterate over a range of integers (or other whatever) a little more cleanly:

template <class T>
class xrange_t {
    T start;
    T stop;
public:
    xrange_t(T start, T stop) : start(start), stop(stop) {}

    class iterator : public std::iterator<std::forward_iterator_tag, T> {
        T current;
    public:
        iterator(T t) : current(t) {}
        T operator *() { return current; }
        iterator &operator++() { ++current; return *this; }
        bool operator!=(iterator const &other) const { return current != other.current; }
        bool operator==(iterator const &other) const { return current == other.current; }
    };

    iterator begin() { return iterator(start); }
    iterator end() { return iterator(stop); }
};

template <class T>
xrange_t<T> xrange(T start, T stop) {
    return xrange_t<T>(start, stop);
}

This is a fair amount of code, but I keep it around in a header, so in most normal files, using it just requires one extra line of code. From there, I'd write a simulation class something on the general order of:

class simulation {
public:
    simulation(int total_students, double threshold, double rate) 
        : infections(5),
        uninfected_students(total_students-1),
        total_students(total_students),
        threshold(threshold),
        rate(rate)
    {
        infections[0] = 1;
    }

    int operator()() {
        while (uninfected_students > 0 && 
            std::all_of(infections.begin(), infections.end(), [](int i) { return i>0; }))
        {
            simulate_day();
        }
        return total_students - std::max(uninfected_students, 0);
    }

    friend std::ostream &operator<<(std::ostream &os, simulation const &s) {
        std::cout << "infected: " << 100.0 * s.infections[0] / s.total_students << "%\n";
        std::cout << "P = " << 100.0 * (s.infections[0] -1) / s.total_students;
    }

private:
    void simulate_day() {
            // ...
    }

    int uninfected_students;
    const int total_students;
    double threshold;
    double rate;
    std::deque<int> infections;
};

I haven't written the code for simulate_day (nor, obviously, tested the code for simulation), but I think the general idea starts to become apparent. It would probably use xrange for most of the loops, so they'd end up something like for (auto i : xrange(0, contagious_students).

The standard algorithms and collections would simplify parts of the code quite a bit though (IMO, anyway). Just for one simple example, your:

memmove(infected_students + 1, infected_students, sizeof(*infected_students) * 4);
infected_students[0] = newly_infected_students;

...could end up something like:

infections.pop_back();
infections.push_front(newly_infected_students);

At least to me, this seems quite a bit more readable--and there's even a decent chance that it'll be faster (it's \$O(1)\$ instead of \$O(N)\$, but in this case N is pretty small, so it probably won't make a big difference either way).

Anyway, using that would look something like this:

int main(){
    static const int num_simulations = 1000;

    double rate = 0.04;
    double threshold = RAND_MAX * 0.04 * rate;

    int students_infected = 0;
    int successful_infections = 0;

    simulation sim(0, threshold, rate);

    for (auto s : xrange(0, num_simulations)) {
        int simulations = sim();
        students_infected += simulations;
        successful_infections += (simulations >= 10);
    }       

    std::cout << sim;
    std::cout << "students infected : " << students_infected
              << "\n\"successful\" infections: " << successful_infections;
}

I reemphasize: this is really just a sketch, not finished code by any means. Especially in main, I'm pretty sure I missed at least some parts of what the code does, but I think this is at least close to the general idea of how a simulation can/would work. At least for the moment, I'm more concerned with overall structure than with the details of each part.

Answer 2

Reduce global variables

It might be only a toy application but for the sake of future growth (and multithreading) it is a good idea to avoid global variables.

Even for C code I find it inappropriate to have them. I would place these variables in an own class Simulation and create an instance of this class for each run. That way you can have multiple Simulations in parallel without them disturbing each other.

Use a commandline parser

Your command line interface is very small so it is feasible to do this by hand but I would recommend using a library that does the parsing for you. This way the interface is easier to extend/maintain.

Seed parameter

It is more of a usability issue but I would allow for the user to (optionally) pass the seed for the random number generator so results can be repeated (which is nice for debugging).

Use a better RNG

If you can use C++11 (or boost) I would recommend using the std::mt19937 rng and the appropriate distributions for your usecase.

Avoid magic numbers

There are several place where you use numbers in your program without explanation of their meaning. Make them named constants and the program will become much more understandable. One very suspicious is:

uninfected_students = 299;
//...a
return 300 - (uninfected_students < 0 ? 0 : uninfected_students);

Which looks like the total number of students is 300. Now imagine you want to have 400 students. You must remember to change all places where the number300 occurs but you might miss the 299. Instead you should have something like:

// we start with one student being infected
uninfected_students = total_number_of_students - 1;
//...
return total_number_of_students - (uninfected_students < 0 ? 0 : uninfected_students);

printf

There are definitely times when printf is more readable than the iostreams but that does not mean you have to use printf. You should have a look into boost.format which offers a very similar interface to printf and interoperates with iostreams.