人口密集領域の階層化

著者：畝見達夫 作成: 令和4年10月4日，編集: 10月5日．

SimEpidemic 個体ベース感染シミュレータ で用いる， 人口密集領域の階層化方法について説明する。

目的

人口統計などから２次元平面上の人口分布を得ることができるが， 感染症の伝播をモデル化する際には，人口密度は参考に過ぎず， 伝播経路となる人と人の接触をより適切に考えるには，多階層建築物内の階層を考慮する必要がある。 どんなに２次元空間上の密度が高くても，異なるフロアにいる人々の間で感染が広がるとは考えにくい。

正確にモデル化するには，地域ごとの居住可能な建屋の階層数の分布に関する情報や， さらには部屋割りの情報なども必要だが，これを利用可能なデータとしてまとめたものの入手は困難なため， ここでは便宜的に，限界密度をパラメータとして定め，その値を超えないよう階層化する。

アルゴリズムの概要

分布に従う位置の決定アルゴリズム を用いて 指定された密度分布に従う２次元平面上の点の集合を生成する。 この時点では，すべての点は同じ平面つまり地上階にあるものとする。 上述の限界密度を２点間の最小距離 \(\underline{d}\) として与える。 各点 \(a\) について，周囲の点 \(b\) との距離 \(d\) を計算し，\(d < \underline{d}\) なら点 \(b\) を1つ上の階に移動する。 上の階に移動した各点について同様の処理を実行し， 同じ階にあるすべての２点間の距離が \(\underline{d}\) 未満となるまで繰り返す。

距離を調べる順序により，結果として各点に与えられる階の分布は異なる。

空間分割による高速化

都市人口などを扱う場合には計算対象となる点の数が10万を超えることになる。 全体の計算時間を合理的な時間内に収めるには，上記の概要に示した２点間の距離を調べる計算について， 網羅的な計算では \(N\) を点の数とすると \(\mbox{O}(N^2)\) の計算量が必要となり， 合理的な時間内に計算が終了しない可能性が出てくる。

ここでは２次元平面を固定された均等な正方形の部分領域に分割し， 同じ区画および周囲の８近傍の区画に属する２点だけを調べることで高速化を図る。 対象を3×3の9区画に抑えるため，区画にあたる正方形の辺の長さを最小距離 \(\underline{d}\) 以上とする。 干渉する可能性のある区画の組み合わせについて網羅的に計算することになるが， 区画の数を \(M\) とすると１区画あたりの点の数の平均は \(N/M\)， 2区画間の隣接関係は，上下，左右，斜め×2 の合計4とおりあり， それに区画内の2点の関係を含めるた 5とおりについて，すべての関係を計算する必要があるから， 計算量は \(\mbox{O}(N^2\cdot 5/M)\) となる。 たとえば \(50\times 50\) の格子を用いれば \(M=2,500\) であるから， 計算量を \(M/5=500\) 分の一に抑えられる。 互いに干渉しない区画間での計算は並列に実行できるので， 複数のスレッドを用いて応答時間の短縮を図ることもできる。

Processing (JAVA) のサンプルコード

int N = 100000; // The number of points.
int NMesh = 50;
float MinDist = 3;
float minx = 1e10, maxx = -1e10, miny = 1e10, maxy = -1e10;
int flmax = -1;
MyPlace[] pts = new MyPlace[N];
class MyPlace {
  float x, y;
  int fl;
}
void findSpanXY(MyPlace[] pts) {
  for (MyPlace p : pts) {
    if (minx > p.x) minx = p.x;
    if (maxx < p.x) maxx = p.x;
    if (miny > p.y) miny = p.y;
    if (maxy < p.y) maxy = p.y;
  }
}
void setup() {
  // *** setup the 2D coordinates for all members in pts.
  findSpanXY(pts);
  int fl = makeStack(pts);
}
int make_stack(MyPlace[] pts) {
  MeshWorld mWorld = new MeshWorld(NMesh);
  mWorld.setupData(pts, pts.length, minx, maxx, miny, maxy);
  MakeStack ms = new MakeStack();
  do { ms.checkDistances(mWorld); } while (ms.n > 0);
  return ms.fl;
}

世界を正方格子に分割するオブジェクト

class MeshWorld {
  ArrayList<MyPlace> mesh[][];
  float minX, maxX, minY, maxY;
  int meshSz;
  MeshWorld(int nMesh) {
    mesh = new ArrayList[nMesh][];
    for (int i = 0; i < nMesh; i ++)
      mesh[i] = new ArrayList[nMesh];
    meshSz = nMesh;
  }
  private float indexF(float v, float minV, float maxV) {
    return (v - minV) / (maxV - minV) * meshSz;
  }
  private int index(float v, float minV, float maxV) {
    int idx = int(indexF(v, minV, maxV));
    return (idx < 0)? 0 : (idx >= meshSz)? meshSz - 1 : idx;
  }
  float indexXF(MyPlace p) { return indexF(p.x, minX, maxX); }
  float indexYF(MyPlace p) { return indexF(p.y, minY, maxY); }
  int indexX(MyPlace p) { return index(p.x, minX, maxX); }
  int indexY(MyPlace p) { return index(p.y, minY, maxY); }
  void setupData(MyPlace[] pts, int n, float x1, float x2, float y1, float y2) {
    minX = x1; maxX = x2; minY = y1; maxY = y2;
    for (int i = 0; i < n; i ++) {
      MyPlace p = pts[i];
      int idxX = indexX(p);
      int idxY = indexY(p);
      if (mesh[idxX][idxY] == null)
        mesh[idxX][idxY] = new ArrayList();
      mesh[idxX][idxY].add(p);
    }
  }
}

２次元平面上に分布する点を階層化するオブジェクト

class MakeStack {
  int n = 0, fl = 0;
  void checkDistance(MyPlace a, MyPlace b) {
    if (a.fl != fl || b.fl != fl) return;
    float d = dist(a.x, a.y, b.x, b.y);
    if (d < MinDist) { b.fl ++; n ++; }
  }
  private void checkDistancesInLists(ArrayList<MyPlace> A, ArrayList<MyPlace> B) {
    if (A != null && B != null)
      for (MyPlace a : A) for (MyPlace b : B) checkDistance(a, b);
  }
  void checkDistances(MeshWorld mWorld) {
    n = 0;
    for (ArrayList<MyPlace>[] listList : mWorld.mesh)
      for (ArrayList<MyPlace> list : listList) {
        if (list != null) {
          int sz = list.size();
          for (int k = 0; k < sz - 1; k ++) {
            MyPlace p = list.get(k);
            for (int l = k + 1; l < sz; l ++) checkDistance(p, list.get(l));
          }
        }
      }
    for (int i = 1; i < mWorld.meshSz; i ++)
      for (int j = 0; j < mWorld.meshSz; j ++) {
        checkDistancesInLists(mWorld.mesh[j][i-1], mWorld.mesh[j][i]);
        checkDistancesInLists(mWorld.mesh[i-1][j], mWorld.mesh[i][j]);
      }
    for (int i = 1; i < mWorld.meshSz; i ++)
      for (int j = 1; j < mWorld.meshSz; j ++) {
        checkDistancesInLists(mWorld.mesh[j-1][i-1], mWorld.mesh[j][i]);
        checkDistancesInLists(mWorld.mesh[j][i-1], mWorld.mesh[j-1][i]);
      }

    for (int i = 0; i < mWorld.meshSz; i ++)
      for (int j = 0; j < mWorld.meshSz; j ++) if (mWorld.mesh[i][j] != null) {
        ArrayList<MyPlace> list = mWorld.mesh[i][j];
        for (int k = list.size() - 1; k >= 0; k --)
          if (list.get(k).fl <= fl) list.remove(k);
        if (mWorld.mesh[i][j].size() == 0) mWorld.mesh[i][j] = null;
      }
    fl ++;
  }
}

ここでは明示的な並列化を行なっていない。