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

我有一个家庭作业，我必须编写一个C++程序，使用SIR模型(易感、传染性、恢复)来模拟疾病爆发。要求使用7x7大小的2D阵列，用户可以选择X和Y坐标来初始化感染者。如果附近有感染者，易感人群(S)将被感染(I)。如果附近有康复者，感染者将康复(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数组，如果其他，循环只能解决这个问题。Map和Vector远远超出了我目前所知的xD。

这项作业让我想起了在大学的日子(那是很久以前的事了)。这似乎是康威的《生命的游戏》的变体，我在大一时收到了这本书。因此，我无法抗拒。。。

之前的一些注意事项：

二维数组在C++中有点不方便。要么必须使用常量大小，要么在不使用某种new[](或g++的VAL扩展，它不符合标准)的情况下无法调整大小。更好的选择通常是std::vector。可以通过适当的运算符重载来"伪造"这两个维度，而不是嵌套std::vector。幸运的是，我手头有一个最低限度的工作版本，这是我最近对多线程基准测试问题的另一个回答。
关于模拟步骤i，我得出了以下逻辑：
如果患者X是
- 's'：检查她/他周围的所有邻居是否有人感染('i')。如果是，感染患者X
- 'i'(前一天感染)：让她/他康复('r')
- 'r'(康复)：对他不做任何事，即保持她/他康复('r')
请注意，不同当前情况的测试可以在板的所有行/所有列的一次迭代中完成–没有必要在单独的功能中这样做。
最有趣的案例是's'。对于[i][j]的患者X，必须检查所有邻居。这些是[i+iP][j+jP]的患者，iP在[-1，1]，jP在[-1、1]。当iP==0和jP==0时，对这9个值进行迭代将检查患者X本身。这种特殊情况可以检查，但我忽略了它(按照上面的逻辑)，病人不能感染自己。这在最内部的循环中为iP和jP节省了额外的检查，这是IMHO欢迎的。
仔细一看，您会发现，如果i==0或i==行数-1或j==0或j==列数-1，[i+iP][j+jP]可能会导致无效坐标。这将需要大量额外的测试来授予有效的指数，但我使用了另一个技巧：我分别将板变大，以提供一个边界。我不使用它进行写作，但这为我提供了安全的阅读访问。我所要承认的是，从这些边界单元格中读取不会篡改我的模拟逻辑。我用'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.

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

您正在使用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;
}

这个函数的作用是同时为您提供每天的地图，现在您可以一个接一个地迭代它们。。。

注意：此外，不要忘记在std::adjacent_find()的开头使用#include <algorithm>

示例：

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

编辑：在coliru上直播(使用以中心为感染点的9x9阵列)

好吧，看看它是否能给出你想要的输出。。。

问候，

Ruks。

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

相关内容

最新更新

热门标签：