Disease Outbreak Simulation using SIR model

Question

I have a homework where I have to write a C++ program to simulate a disease outbreak using SIR model (Susceptible, Infectious, Recover). The requirement is to use a 2D-array with 7x7 size where user will choose an X and Y coordinate to initialize an infectious person. A Susceptible person (S) will become infected (I) if there is an infected person in adjacent. Then an infected person will recover (R) if there is a Recover person in adjacent. The program will end if all people are recovered. Example output:

Day 0                          Day 1                       Day 2
s s s s s s s                  s s s s s s s               s s s s s s s
s s s s s s s                  s s s s s s s               s i i i i i s
s s s s s s s                  s s i i i s s               s i r r r i s
s s s i s s s                  s s i r i s s               s i r r r i s
s s s s s s s                  s s i i i s s               s i r r r i s
s s s s s s s                  s s s s s s s               s i i i i i s
s s s s s s s                  s s s s s s s               s s s s s s s

So far, I can only check the state in position (1,1), (1,7), (7,1), (7,7). If the next three position next to it have an infected person, it will update the state to nextDayState. Here is my code so far for two functions, SpreadingDisease and RecoverState.

    void recoverState(char currentDayState[SIZE][SIZE], char nextDayState[SIZE][SIZE], int sizeOfArray)//It will take in the currentState of Day 0. I also copy the elements in currentState to nextDayState so that it could work. 
{
    for (int i = 1; i < sizeOfArray + 1; ++i)
    {
        for (int j = 1; j <= sizeOfArray + 1; ++j)
        {
            if (currentDayState[i][j] == 'i')//If found any Infected, update it to Recover on the nextDayState array. 
            {
                nextDayState[i][j] == 'r';
            }
        }
    }

    for (int i = 1; i < sizeOfArray + 1; ++i)
    {
        for (int j = 1; j <= sizeOfArray + 1; ++j)
        {
            currentDayState[i][j] = nextDayState[i][j];
            //After all people are recover, update the currentState and output it to terminal. 
        }
    }
}
void spreadDisease(const char currentDayState[SIZE][SIZE], char nextDayState[SIZE][SIZE], int sizeOfArray, int day = 1)
{
    for (int i = 1; i < sizeOfArray + 1; ++i)
    {
        for (int j = 1; j <= sizeOfArray + 1; ++j)
        {
            if (currentDayState[i][j] == 's')
            {
                if (i == 1 && j == 1)
                {
                    if (currentDayState[1][2] == 'i' || currentDayState[2][1] == 'i' || currentDayState[2][2] == 'i')
                    {
                        nextDayState[1][1] = 'i';
                    }
                }
                if (i == 1 && j == 7)
                {
                    if (currentDayState[1][6] == 'i' || currentDayState[2][6] == 'i' || currentDayState[2][7] == 'i')
                    {
                        nextDayState[1][7] = 'i';
                    }
                }
                if (i == 7 && j == 1)
                {
                    if (currentDayState[6][1] == 'i' || currentDayState[6][2] == 'i' || currentDayState[7][2] == 'i')
                    {
                        nextDayState[7][1] = 'i';
                    }
                }
                if (i == 7 && j == 7)
                {
                    if (currentDayState[6][6] == 'i' || currentDayState[7][6] == 'i' || currentDayState[6][7] == 'i')
                    {
                        nextDayState[7][7] = 'i';
                    }
                }
            }
        }
    }
}

I figure out that If I can somehow get the X and Y coordinate from the user, then I can use that coordinate to update the state of the next day. Unfortunately, I don't know how to assign the X and Y coordinate into the function to start with it.

P/S: Thank you for all of your answers. I very appreciate your kindness. However, I should have mentioned the requirement of my assignment before. Since I only study till the User-Defined Functions part, I am not allowed to use anything else beyond that. So I am limited to use 2D-array, If-else, Looping only to solve this problem. Map and Vector is far beyond my knowledge right now xD.

score 1 · Answer 1 · answered Oct 18 '18 at 07:36

This assignment remembered me to my days at the University (and that's quite long ago). It seems like a variant of Conway's Game of Life which I got as assignment when I was first year's student. Hence, I couldn't resist...

Some notes before:

Two dimensional arrays are a bit inconvenient in C++. Either you have to use constant size or resizing them is not possible without using some kind of new[] (or g++'s VAL extension which is not standard conform). The better alternative is usually std::vector. Instead of nesting std::vectors, the two dimensions can be "faked" by appropriate operator overloads. For my luck, I had a minimal working version at hand from another recent answer of mine to Multi-threading benchmarking issues.
Concerning the simulation step i, I came to the following logic:
If patient X is
- 's': check all neighbours around her/him whether somebody is infected ('i'). If so, infect patient X.
- 'i' (infected in day before): let her/him recover ('r').
- 'r' (recovered): do nothing with him i.e. keep her/him recovered ('r').
Please, note that the tests of the different current cases can be done in one iteration of all rows/all columns of the board – no necessity to do this in separate functions.
The most interesting case is 's'. For patient X at [i][j], all neighbours have to be checked. These are the patients at [i + iP][j + jP] with iP in [-1, 1] and jP in [-1, 1]. Iterating over these 9 values will check the patient X itself when iP == 0 and jP == 0. This special case could be checked but I ignored it (as by the above logic) a patient cannot infect itself. This saves an extra check for iP and jP in the most inner loop which is IMHO welcome.
At closer glance, you will realize that [i + iP][j + jP] might result in invalid coordinates if i == 0 or i == number of rows - 1 or j == 0 or j == number of columns - 1. This would require a lot of additional tests to grant valid indices but I use another trick: I make the board respectively larger to provide a border around. I don't use it for writing but this provides me safe read-accesses. All I have to grant is that reading from these border cells will not tamper my simulation logic. I initialize the whole board including border cells with 's'. As border cells are never written (except in initialization) they are never infected what matches my concept.

So, this is my simulation step:

void doSimStep(const Board &board, Board &board1)
{
  assert(board.getNumRows() == board1.getNumRows());
  assert(board.getNumCols() == board1.getNumCols());
  for (size_t i = 1, nRows = board.getNumRows() - 1; i < nRows; ++i) {
    for (size_t j = 1, nCols = board.getNumCols() - 1; j < nCols; ++j) {
      const char person = board[i][j];
      char person1 = person;
      switch (person) {
        case 's': { // search for infection in neighbourhood
          bool infect = false;
          for (int iP = -1; !infect && iP <= 1; ++iP) {
            for (int jP = -1; !infect && jP <= 1; ++jP) {
              infect = board[i + iP][j + jP] == 'i';
            }
          }
          person1 = infect ? 'i' : 's';
        } break;
        case 'i': // infected -> recover
          // fall through
        case 'r': // recovered: stable state
          person1 = 'r';
          break;
        default: assert(false); // Wrong cell contents!
      }
      board1[i][j] = person1;
    }
  }
}

I don't get why user10522145 believes this cannot be done without recursion. (Btw., I believe the opposite: every recursion can be turned into an iteration which may accumulate or stack intermediate results.) I actually don't know where the recursion would be necessary considering that the OP already planned separate boards for current and new state (which simplifies things much).

Output of a simulation with a 9×9 board:

Init.:
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s

Day 0:
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s i s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s

Day 1:
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s i i i s s s
 s s s i r i s s s
 s s s i i i s s s
 s s s s s s s s s
 s s s s s s s s s
 s s s s s s s s s

Day 2:
 s s s s s s s s s
 s s s s s s s s s
 s s i i i i i s s
 s s i r r r i s s
 s s i r r r i s s
 s s i r r r i s s
 s s i i i i i s s
 s s s s s s s s s
 s s s s s s s s s

Day 3:
 s s s s s s s s s
 s i i i i i i i s
 s i r r r r r i s
 s i r r r r r i s
 s i r r r r r i s
 s i r r r r r i s
 s i r r r r r i s
 s i i i i i i i s
 s s s s s s s s s

Day 4:
 i i i i i i i i i
 i r r r r r r r i
 i r r r r r r r i
 i r r r r r r r i
 i r r r r r r r i
 i r r r r r r r i
 i r r r r r r r i
 i r r r r r r r i
 i i i i i i i i i

Day 5:
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r
 r r r r r r r r r

No further progress detected on day 6.
Done.

Live Demo on coliru

And finally (spoiler alert) the complete source code:

#include <cassert>
#include <iomanip>
#include <iostream>
#include <vector>

template <typename VALUE>
class MatrixT; // forward declaration

template <typename VALUE>
void swap(MatrixT<VALUE>&, MatrixT<VALUE>&); // proto

template <typename VALUE>
class MatrixT {
  friend void swap<VALUE>(MatrixT<VALUE>&, MatrixT<VALUE>&);
  public:
    typedef VALUE Value;
  private:
    size_t _nRows, _nCols;
    std::vector<Value> _values;
  public:
    MatrixT(size_t nRows, size_t nCols, Value value = (Value)0):
      _nRows(nRows), _nCols(nCols), _values(_nRows * _nCols, value)
    { }
    ~MatrixT() = default;
    MatrixT(const MatrixT&) = default;
    MatrixT& operator=(const MatrixT&) = default;

    size_t getNumCols() const { return _nCols; }
    size_t getNumRows() const { return _nRows; }
    const std::vector<Value>& get() const { return _values; }
    Value* operator[](size_t i) { return &_values[0] + i * _nCols; }
    const Value* operator[](size_t i) const { return &_values[0] + i * _nCols; }
};

template <typename VALUE>
void swap(MatrixT<VALUE> &mat1, MatrixT<VALUE> &mat2)
{
  std::swap(mat1._nRows, mat2._nRows);
  std::swap(mat1._nCols, mat2._nCols);
  std::swap(mat1._values, mat2._values);
}

typedef MatrixT<char> Board;

bool operator==(const Board &board1, const Board &board2)
{
  return board1.getNumRows() == board2.getNumRows()
    && board1.getNumCols() == board2.getNumCols()
    && board1.get() == board2.get();
}

std::ostream& operator<<(std::ostream &out, const Board &board)
{
  for (size_t i = 1, nRows = board.getNumRows() - 1; i < nRows; ++i) {
    for (size_t j = 1, nCols = board.getNumCols() - 1; j < nCols; ++j) {
      out << ' ' << board[i][j];
    }
    out << '\n';
  }
  return out;
}

void doSimStep(const Board &board, Board &board1)
{
  assert(board.getNumRows() == board1.getNumRows());
  assert(board.getNumCols() == board1.getNumCols());
  for (size_t i = 1, nRows = board.getNumRows() - 1; i < nRows; ++i) {
    for (size_t j = 1, nCols = board.getNumCols() - 1; j < nCols; ++j) {
      const char person = board[i][j];
      char person1 = person;
      switch (person) {
        case 's': { // search for infection in neighbourhood
          bool infect = false;
          for (int iP = -1; !infect && iP <= 1; ++iP) {
            for (int jP = -1; !infect && jP <= 1; ++jP) {
              infect = board[i + iP][j + jP] == 'i';
            }
          }
          person1 = infect ? 'i' : 's';
        } break;
        case 'i': // infected -> recover
          // fall through
        case 'r': // recovered: stable state
          person1 = 'r';
          break;
        default: assert(false); // Wrong cell contents!
      }
      board1[i][j] = person1;
    }
  }
}

int main()
{
  size_t nRows = 9, nCols = 9;
#if 0 // disabled for demo
  std::cout << "N Rows: "; std::cin >> nRows;
  std::cout << "N Cols: "; std::cin >> nCols;
  /// @todo check nRows, nCols for sufficient values
#endif // 0
  // init board
  std::cout << "Init.:\n";
  Board board(nRows + 2, nCols + 2);
  std::fill(board[0], board[nRows + 2], 's');
  std::cout << board << '\n';
  // infect somebody
  size_t i = nRows / 2 + 1, j = nCols / 2 + 1;
#if 0 // disabled for demo
  std::cout << "Patient 0:\n";
  std::cout << "row: "; std::cin >> i;
  std::cout << "col: "; std::cin >> j;
  /// @todo check i, j for matching the boundaries
#endif // 0
  board[i][j] = 'i';
  // simulation loop
  for (unsigned day = 0;;) {
    std::cout << "Day " << day << ":\n";
    std::cout << board << '\n';
    // simulate next day
    ++day;
    Board board1(board);
    doSimStep(board, board1);
    if (board == board1) {
      std::cout << "No further progress detected on day "
        << day << ".\n";
      break; // exit sim. loop
    }
    // store data of new day
    swap(board, board1);
  }
  // done
  std::cout << "Done.\n";
  return 0;
}

It didn't exceed the limits of an SO answer. May be, things could be done shorter changing the requirements or making things less readable. What do you would change? — Scheff's Cat, Oct 18 '18 at 09:34

Ruks · Answer 2 · 2018-11-04T12:54:44.657

You are using C++, so use the standard library to the maximum...

The magically optimized disease simulation function:

/*
 *-----------------------
 * Key:
 * ----------------------
 * 0 - Susceptible person
 * 1 - Infected person
 * 2 - Recovered person
 * 
 * @param init_infect_x Person to infect at x position...
 * @param init_infect_y Person to infect at y position...
 * @param map_size_x Width of the map...
 * @param map_size_y Height of the map...
 */
std::vector<std::vector<std::vector<int>>> disease_simulator(size_t const init_infect_x = 0u,
                                                             size_t const init_infect_y = 0u,
                                                             size_t const map_size_x = 7u, size_t const map_size_y = 7u)
{
    if (map_size_x == 0u || map_size_y == 0u || init_infect_x + 1 > map_size_x || init_infect_x + 1 < 0 || init_infect_y
        + 1 > map_size_y || init_infect_y + 1 < 0) // Well, we can't create a map which is empty...
        return std::vector<std::vector<std::vector<int>>>();
    std::vector<std::vector<std::vector<int>>> map_list;
    std::vector<std::pair<int, int>> spread_pos;
    std::vector<std::vector<int>> map(map_size_y, std::vector<int>(map_size_x, 0));
    map[init_infect_y][init_infect_x] = 1;
    map_list.emplace_back(map);
    while (std::adjacent_find(map.begin(), map.end(), std::not_equal_to<>()) != map.end())
    {
        for (auto i = 0; i < signed(map.size()); i++)
            for (auto j = 0; j < signed(map[i].size()); j++)
                if (map[i][j] == 1)
                {
                    map[i][j] = 2;
                    spread_pos.emplace_back(std::make_pair(j, i));
                }
        for (auto const pos : spread_pos)
        {
            if (pos.second - 1 >= 0 && map[pos.second - 1][pos.first] == 0) // Up...
                map[pos.second - 1][pos.first] = 1;
            if (pos.first - 1 >= 0 && map[pos.second][pos.first - 1] == 0) // Left...
                map[pos.second][pos.first - 1] = 1;
            if (pos.second - 1 >= 0 && pos.first - 1 >= 0 && map[pos.second - 1][pos.first - 1] == 0) // Up left...
                map[pos.second - 1][pos.first - 1] = 1;
            if (pos.second - 1 >= 0 && pos.first + 2 <= signed(map_size_x) && map[pos.second - 1][pos.first + 1] == 0)
                // Up right...
                map[pos.second - 1][pos.first + 1] = 1;
            if (pos.second + 2 <= signed(map_size_y) && map[pos.second + 1][pos.first] == 0) // Down...
                map[pos.second + 1][pos.first] = 1;
            if (pos.first + 2 <= signed(map_size_x) && map[pos.second][pos.first + 1] == 0) // Right...
                map[pos.second][pos.first + 1] = 1;
            if (pos.second + 2 <= signed(map_size_y) && pos.first + 2 <= signed(map_size_x) && map[pos.second + 1][pos.
                first + 1] == 0) // Down right...
                map[pos.second + 1][pos.first + 1] = 1;
            if (pos.second + 2 <= signed(map_size_y) && pos.first - 1 >= 0 && map[pos.second + 1][pos.first - 1] == 0)
                // Down left...
                map[pos.second + 1][pos.first - 1] = 1;
        }
        map_list.emplace_back(map);
        spread_pos.clear();
    }
    return map_list;
}

What this function does is that it gives you the map of each day simultaneously, now you can just iterate over them one by one...

Note: Also, don't forget to #include <algorithm> at the beginning for std::adjacent_find()...

Example:

int main()
{
    auto days_map = disease_simulator();
    for (auto i = 0u; i < days_map.size(); i++)
    {
        std::cout << "Day " << i << ":" << std::endl;
        for (auto elem2 : days_map[i])
        {
            for (auto elem3 : elem2)
                switch (elem3)
                {
                case 0:
                    std::cout << "s ";
                    break;
                case 1:
                    std::cout << "i ";
                    break;
                case 2:
                    std::cout << "r ";
                    break;
                default:
                    std::cout << ' ';
                    break;
                }
            std::cout << std::endl;
        }
        std::cout << std::endl;
    }
    std::cout << "All people have recovered!" << std::endl;
    return 0;
}

Edit: Live on coliru (Using 9x9 arrays with center as the infect point)

Well, see if it gives your desired output...

Kind regards,

Ruks.

Can you demonstrate that it works (e.g. on [**coliru**](http://coliru.stacked-crooked.com/))? — Scheff's Cat, Oct 18 '18 at 09:43
I edited it to add demonstration with Coliru there, check it out! — Ruks, Oct 18 '18 at 09:47

score 0 · Answer 3 · answered Oct 18 '18 at 05:32

0

I guess iteration might not work in this case ill suggest you use recursion with the array boundary values as the condition to stop recursion. Hope it made sense

answered Oct 18 '18 at 05:32

Lets assume infection starts at 1,1 your code will recover it on day one but if you try doing it on paper you will see it should recover on day 2. Appart from that the code only made for corner cases as you told, I am not suggesting that it is not at all possible with iteration but it might get complex and inefficient appart from that this kind of problem as you solve on paper is naturally recursive as you might have figured. – Oct 18 '18 at 06:06

Disease Outbreak Simulation using SIR model

3 Answers3