patterncppModerate

Monte Carlo virus infection simulation

Submitted by: @import:stackexchange-codereview·Mar 10, 2026·

Viewed 0 times

montesimulationviruscarloinfection

Problem

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
#include
#include
#include
#include

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 3) {
printf("usage: maine [simulations=1000]\n");
exit(EXIT_SUCCESS);
}

double transmission_rate = strtod(argv[1], NULL);
if (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 = 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_

Solution

@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 xrange_t {
    T start;
    T stop;
public:
    xrange_t(T start, T stop) : start(start), stop(stop) {}

    class iterator : public std::iterator {
        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 
xrange_t xrange(T start, T stop) {
    return xrange_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 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.

Code Snippets

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);
}

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;
};

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

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

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;
}

Context

StackExchange Code Review Q#55432, answer score: 11

Revisions (0)

No revisions yet.