セルオートマトンでシミュレーションをしよう(7)

疫病の伝播

疫病の伝播モデルである。モデルの詳細は「森林火災と疫病伝播」のコンテンツを参照。ここでは、具体的なプログラムを考える

各セルの状態数はいくつとするべきだろうか。健康・感染・免疫の3状態でいいだろうか。 そうではなく、感染状態と免疫状態はそれぞれnステップとmステップ続くことを表現するにはさらに感染状態をn通り、免疫状態をm通りに分けるのが簡単である。つまり、セルの状態数は1+n+mとする。状態0は健康、1からnは感染状態、n+1からn+mは免疫状態を表わす。 状態1からn+m-1までは１ステップ進むと数が1増えることにしておけば、感染状態がnステップ続いた後に免疫状態に移行し、それがmステップ続くことを表現できる。 状態n+mは次の時刻で状態0に移行することにすれば、免疫終了後に健康状態に戻る。 状態0の場合だけ、周囲のセルの状態が必要となる。以上を下の表にまとめる。

健康感染感染・・・感染免疫免疫・・・免疫
現在012・・・nn+1n+2・・・n+m
次周囲による23・・・n+1n+2n+3・・・0

状態0以外は簡単で、1を足してn+m+1で割った余りを計算すればいいだけだから 

        if(s[i][j] >= 1)snew[i][j] = (s[i][j]+1)%(n+m+1);
でよい

状態0では、以下のふたつに場合わけする

  1. 周囲にひとつも感染状態のセルがない。この場合、セルの状態は0のまま変わらず
  2. 周囲に感染状態のセルがひとつ以上ある。この場合は確率pで感染状態に移行する
感染状態のセルがひとつ以上あるかどうかの判定は、素直に書くとかなり長たらしくなる 
        else if((s[i-1][j-1] >= 1 && s[i-1][j-1] <= n) ||
                (s[i-1][j]   >= 1 && s[i-1][j]   <= n) ||
                (s[i-1][j+1] >= 1 && s[i-1][j+1] <= n) ||
                (s[i][j-1]   >= 1 && s[i][j-1]   <= n) ||
                (s[i][j+1]   >= 1 && s[i][j+1]   <= n) ||
                (s[i+1][j-1] >= 1 && s[i+1][j-1] <= n) ||
                (s[i+1][j]   >= 1 && s[i+1][j]   <= n) ||
                (s[i+1][j+1] >= 1 && s[i+1][j+1] <= n))
&&は「かつ」、||は「または」である。周囲の8セルすべてについて、値が1以上n以下であるかどうかを判定し、それがひとつでも成立すれば、「感染状態のセルがひとつ以上」ということになる。

「確率pで感染状態1に移行」を実現するには、乱数を発生してpと比較する。ただし、random()は整数乱数なので、pといきなり比較しても意味がない。実はrandom()の最大値にはRAND_MAXという名前がつけられている（stdlib.hの中で定義される）。つまり、random()は0からRAND_MAXまでの整数乱数を発生する。というわけで、random()とp*RAND_MAXとの大小を比較すればよいことになる。この部分は 

　if(random() < p*RAND_MAX) snew[i][j] = 1;
と書けばよい。初期状態を作る部分も同様にできる

ほかの部分はほとんどGame of lifeそのままでよいから、プログラム全体は以下のようになるだろう。 
#include<stdio.h>
#include<eggx.h>
#include<stdlib.h>
#define NX 100
#define NY 100
#define EDGE 4
main()
{ 
  int s[NX+2][NY+2], snew[NX+2][NY+2];
  int i,j,k;
  int step,rseed;
  int n,m;
  double p,q,r;
  int win;

/* preparation for graphics */
  win=gopen(NX*EDGE,NY*EDGE);
  window(win,1,1,NX,NY);
  gsetbgcolor(win,"white");
  layer(win,0,1);

/* input step and rseed */
  scanf("%d%d",&step,&rseed);
  srandom(rseed);

/* set parameters */
  p = 0.5;
  q = 0.01;
  r = 0;
  n = 5;
  m = 10;

/* set random initial state */
  for(i=1;i<NX+1;++i){
   for(j=1;j<NY+1;++j){
     if(random()< q*RAND_MAX) s[i][j]= 1;
     else if(random() < (q+r)*RAND_MAX) s[i][j] = n+1;
     else s[i][j] = 0;
  }}

/* transfer boundary */
  for(i=1;i<NX+1;++i){
   s[i][0] = s[i][NY];
   s[i][NY+1] = s[i][1];
  }
  for(j=0;j<NY+2;++j){
   s[0][j] = s[NX][j];
   s[NX+1][j] = s[1][j];
  }

/* display initial state */
   gclr(win);
   newcolor(win, "red");
   for(i=1;i<NX+1;++i){
     for(j=1;j<NY+1;++j){
      if(s[i][j] >= 1 && s[i][j] <=n)fillrect(win,i,j,1,1);}}
   newcolor(win, "green");
   for(i=1;i<NX+1;++i){
     for(j=1;j<NY+1;++j){
      if(s[i][j] >= n+1 && s[i][j] <=n+m) fillrect(win,i,j,1,1);}}

   copylayer(win,1,0);

  for(k=0; k<step; ++k){
/* start one step */
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
/* one site */
        if(s[i][j] >= 1)snew[i][j] = (s[i][j]+1)%(n+m+1); 
        else if((s[i-1][j-1] >= 1 && s[i-1][j-1] <= n) ||
                (s[i-1][j]   >= 1 && s[i-1][j]   <= n) ||
                (s[i-1][j+1] >= 1 && s[i-1][j+1] <= n) ||
                (s[i][j-1]   >= 1 && s[i][j-1]   <= n) ||
                (s[i][j+1]   >= 1 && s[i][j+1]   <= n) ||
                (s[i+1][j-1] >= 1 && s[i+1][j-1] <= n) ||
                (s[i+1][j]   >= 1 && s[i+1][j]   <= n) ||
                (s[i+1][j+1] >= 1 && s[i+1][j+1] <= n))
                 { if(random() < p*RAND_MAX) snew[i][j] = 1;}
      }}
/* transfer snew[][]  to s[][] */
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j)s[i][j] = snew[i][j];
     }
/* transfer boundary */
     for(i=1;i<NX+1;++i){
      s[i][0] = s[i][NY];
      s[i][NY+1] = s[i][1];
     }
     for(j=0;j<NY+2;++j){
      s[0][j] = s[NX][j];
      s[NX+1][j] = s[1][j];
     }

/* display one step */
     gclr(win);
     newcolor(win, "red");
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
       if(s[i][j] >= 1 && s[i][j] <=n)fillrect(win,i,j,1,1);}}
     newcolor(win, "green");
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
       if(s[i][j] >= n+1 && s[i][j] <=n+m) fillrect(win,i,j,1,1);}}

     copylayer(win,1,0);
     msleep(10);
/* end of one step */
}
}

セルの色は健康を背景色の白とし、感染状態を赤、免疫状態を緑で表わしている。盤面を表示するための 

     gclr(win);
     newcolor(win, "red");
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
       if(s[i][j] >= 1 && s[i][j] <=n)fillrect(win,i,j,1,1);}}
     newcolor(win, "green");
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
       if(s[i][j] >= n+1 && s[i][j] <=n+m) fillrect(win,i,j,1,1);}}
の部分では、まず背景色を塗ったのち、赤のセルを全部塗ってしまい、次に緑のセルを全部塗っている。なぜ、セル一個ごとに状態を判定して 
     gclr(win);
     for(i=1;i<NX+1;++i){
      for(j=1;j<NY+1;++j){
       if(s[i][j] >= 1 && s[i][j] <=n)
         { newcolor(win, "red"); fillrect(win,i,j,1,1);}
       else if(s[i][j] >= n+1 && s[i][j] <=n+m) 
         {  newcolor(win, "green"); fillrect(win,i,j,1,1);}
      }}
としなかったかというと、後者は表示が非常に遅くなるためである。newcolor命令を使う回数は極力減らしたほうがいいようだ。

次へ