使用 SIR 模型的疾病爆发模拟

Question

我有一个作业，我必须编写一个 C++ 程序来使用 SIR 模型（易感性、传染性、恢复）来模拟疾病爆发。要求是使用 7x7 大小的 2D 数组，用户将在其中选择 X 和 Y 坐标来初始化感染者。如果附近有感染者，易感者 (S) 将被感染 (I)。如果附近有 Recover 人，则感染者将恢复 (R)。如果所有人都康复，该计划将结束。 示例输出：

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

到目前为止，我只能检查位置 (1,1), (1,7), (7,1), (7,7) 的状态。如果它旁边的下三个位置有感染者，它会将状态更新为 nextDayState。 到目前为止，这是我的两个函数的代码，SpreadingDisease 和 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';
                    }
                }
            }
        }
    }
}

我发现如果我能以某种方式从用户那里获得 X 和 Y 坐标，那么我可以使用该坐标来更新第二天的状态。不幸的是，我不知道如何将 X 和 Y 坐标分配到函数中以开始它。

P/S：感谢您的所有回答。我非常感谢你的好意。但是，我应该在之前提到我的任务要求。因为我只学习到用户定义的函数部分，所以我不能使用除此之外的任何东西。所以我仅限于使用 2D-array、If-else、Looping 来解决这个问题。地图和矢量现在远远超出我的知识范围 xD。

Answer 1

我猜在这种情况下迭代可能不起作用，我建议您使用带有数组边界值的递归作为停止递归的条件。 希望有道理

Answer 2

这项任务让我想起了我在大学的日子（那是很久以前的事了）。 它似乎是Conway's Game of Life 的一个变体，我在一年级时得到的作业。因此，我无法抗拒……

之前的一些笔记：

1. 二维数组在 C++ 中有点不方便。如果不使用某种 new[]（或不符合标准的 g++ VAL 扩展），则必须使用常量大小或调整它们的大小是不可能的。更好的选择通常是std::vector。不是嵌套std::vectors，而是可以通过适当的运算符重载来“伪造”这两个维度。幸运的是，我手头有一个最小的工作版本，来自我最近对Multi-threading benchmarking issues 的另一个回答。

2. 关于模拟步骤 i，我得出以下逻辑：
   如果患者 X 是

   • 's'：检查她/他周围的所有邻居是否有人被感染（'i'）。如果是这样，感染患者 X。
   • 'i'（前一天被感染）：让她/他康复（'r'）。
   • 'r'（已康复）：对他什么也不做，即让她/他康复（'r'）。

   请注意，当前不同案例的测试可以在板的所有行/所有列的一次迭代中完成 - 无需在单独的函数中执行此操作。

3. 最有趣的案例是's'。对于 [i][j] 处的患者 X，必须检查所有邻居。这些是 [i + iP][j + jP] 的患者，iP 在 [-1, 1] 和 jP 在 [-1, 1]。当 iP == 0 和 jP == 0 时，对这 9 个值进行迭代将检查患者 X 本身。可以检查这种特殊情况，但我忽略了它（根据上述逻辑）患者无法感染自己。这在最内层循环中节省了对 iP 和 jP 的额外检查，恕我直言。

4. 仔细一看，如果 i == 0 或 i == 行数 - 1 或 j == 0 或 j，[i + iP][j + jP] 可能会导致无效坐标== 列数 - 1。这将需要大量额外的测试来授予有效索引，但我使用了另一个技巧：我分别使板子更大以提供周围的边框。我不使用它来写作，但这为我提供了安全的读取访问。我必须承认的是，从这些边界单元格中读取数据不会篡改我的模拟逻辑。我用's' 初始化整个板子，包括边界单元格。由于永远不会写入边界单元格（初始化时除外），因此它们永远不会被感染符合我的概念。

所以，这是我的模拟步骤：

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

我不明白为什么user10522145 认为没有递归就无法做到这一点。 （顺便说一句，我相信相反：每次递归都可以变成可能累积或堆叠中间结果的迭代。）考虑到 OP 已经为当前和新计划了单独的板，我实际上不知道在哪里需要递归状态（这大大简化了事情）。

9×9 板的模拟输出：

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

最后（剧透警告）完整的源代码：

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

Answer 3

你使用的是C++，所以尽量使用标准库...

神奇优化的疾病模拟功能：

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

这个函数的作用是它同时给你每天的地图，现在你可以一个一个地迭代它们......

注意： 另外，不要忘记以#include <algorithm>开头为std::adjacent_find()...

示例：

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

编辑： Live on coliru（使用以中心为感染点的 9x9 数组）

好吧，看看它是否能提供你想要的输出......

亲切的问候，

鲁克斯。