LoginSignup
11
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@tommyecguitar(とみー)

社内セミナーで取り上げたMATLAB Tips 【感染拡大シミュレーションを例に】

Posted at

先日弊社社内でMATLABセミナーの講師をする機会があり、テーマ設定から内容(補足)をまとめておきます。
開催にあたって、参加者層はこんな感じ。

  • MATLABの「M」の字くらいは知っている。使用歴は全くなしからプロ級まで様々。
  • 主に制御がしたい人、データ分析がしたい人、趣味に使いたい人が多かった。
  • 役職もバラバラ

要するにいろんな人に参加いただきました。
セミナーの題目としては「趣味のMATLAB」、「MATLAB x 線形代数」を設定しています。一つの処理内容をお題としてライブコーディングをし、MATLABの文法の寛容さ、初学者でも扱いやすいハードルの低さをアピールします。

一応前回ちくっと言われてムッとしたのではっきり申しておきますと
この記事に社内機密は一切ありませんので安心してください！！

具体的な内容は以下の@higoさんの記事を参考にしました。

コロナウイルスの拡散シミュレーション

タイムリーだしMATLABの得意な線形代数が使えると思ったので、これをアレンジ(ほぼ丸コピ)してセミナーに使いました。記事中にソースコードもあったのでそれを見ながらシンプルに書き換えてみました。その時に使ったTipsとセミナーのまとめをメモっておきます。

シミュレーション概要

シミュレーションの内容は元の記事を参考にしていただきたいのですが、ざっと引用してくると以下の条件でシミュレーションします。

まず、このシミュレーションの時間単位は日です。つまり、１ステップで１日が経過します。次に、300×300の格子を考え、各格子点で人を表します。各人は、上下左右と斜めの人（合計8人）と接触しています。保菌者と接触している人は、一定の確率で接触した保菌者からウイルスが移り感染します。保菌者は、規定の日数が経過すると、回復します。回復した人は抗体を獲得し、二度と感染しません。各パラメータは、次のように設定しました。
・感染確率：3%（１日あたり3%の確率で、接触した保菌者から感染）
・回復までの経過日数：10日

動作の様子は元記事をご覧ください。

材料

上記の条件のうちセミナーの題材に取り上げたのは、「あるマスの周囲8マスにいる保菌者数をカウントする」部分です。ここの処理方法を、参加者の皆様でアイデアを出していただきました。
これを求めるために与えられた材料は以下の配列です。

V_{i,j} = \left\{
\begin{array}{l}
0 & Non\ Infected\\
1 & Infected
\end{array}
\right.

MATLAB上ではtrue or falseの論理値の配列ですが、doubleの行列と演算するとdoubleとして動作します。便利。心配ならdouble(X)など使ってdoubleに変換しましょう。

ちなみにセミナーのテーマが「MATLAB x 線形代数」なので、ツールボックスの使用は禁じました。
元記事にもあるように以下の条件も与えています。

  • マスの最外周の人は感染しない。
  • 最初の感染者は回復しない。

どちらも後段で処理します。

出て来たアイデア

セミナーでいただいた処理アイデアを示します。

% InfectIdx = 感染者フラグ、感染ならtrue(=1)、非感染ならfalse(=0)。

AroundInfectSum = zeros(300);

for irow = 2:299
    for icol = 2:299
        AroundInfectSum(irow,icol) = InfectIdx(irow-1,icol-1)+InfectIdx(irow-1,icol)+InfectIdx(irow-1,icol+1)+...
                                     InfectIdx(irow,icol-1)+InfectIdx(irow,icol+1)+...
                                     InfectIdx(irow+1,icol-1)+InfectIdx(irow+1,icol)+InfectIdx(irow+1,icol+1);
    end
end

基本の書き方ですね。自身のマスの上下左右を足算するのをforループで回す方法です。意外にもそこそこ動作が速くてびっくりしましたが綺麗ではないですね。

他には上の例をsumを使ってやる方法も出ました。sum関数の方が記述が綺麗ですが、配列の全成分を参照しに行くので遅かったです。惜しい！

チャット欄にはこんな意見も。。。

forループの呪縛から抜けられないのか。。。

このリアルタイムで考えていただいてる感じ、主催者冥利につきますね！ありがたい！

限られた時間で最適解を出すのは非常に難しいですが、みなさん頭を使って本気で取り組んでくださいました。ありがとうございます。

多分最もシンプルな処理

ツールボックスなしでの最もシンプルな処理式はこれです。$C$が周囲8マスの感染者数を要素に持った配列です。

S_{i,j} = \left\{
\begin{array}{l}
1 & i = j-1, j, j+1\\
0 & Otherwise
\end{array}
\right.\\
C = SVS-V

$S$をちゃんと書くとこういう行列です。

S =
\begin{bmatrix}
1&1&0&0&\cdots&0&1\\
1&1&1&0&\ddots&0&0\\
0&1&1&1&\ddots&0&0\\
0&0&1&1&\ddots&0&0\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\
1&0&0&0&\cdots&1&1
\end{bmatrix}

MATLABだとこれでかけます。

size = 100;
S = circshift(eye(size),1,1) + eye(size) + circshift(eye(size),-1,1);

発想の元になったのは置換行列です。例えば以下の行列と計算例をみてください。

S_1 =
\begin{bmatrix}
0&1&0&\cdots&0\\
0&0&1&\ddots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&0&0&\cdots&0
\end{bmatrix}\\

AS_1 =
\begin{bmatrix}
a_{1,N}&a_{1,1}&a_{1,2}&\cdots&a_{1,N-1}\\
a_{2,N}&a_{2,1}&a_{2,2}&\ddots&a_{2,N-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
a_{N,N}&a_{N,1}&a_{N,2}&\cdots&a_{N,N-1}\\
\end{bmatrix}\\

S_1A =
\begin{bmatrix}
a_{2,1}&a_{2,2}&a_{2,3}&\cdots&a_{2,N}\\
a_{3,1}&a_{3,2}&a_{3,3}&\cdots&a_{3,N}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
a_{1,1}&a_{1,2}&a_{1,3}&\cdots&a_{1,N}\\
\end{bmatrix}\\

$S_1$は単位行列を右もしくは上に一つシフトしたものです。これを対象の行列の後ろから作用させると、行列の要素が右に一つずれます。前から作用させると上に一つずれます。同じように$S_2$を以下のように定義すると、左と下へのシフト計算ができます。

S_2 =
\begin{bmatrix}
0&0&\cdots&0&1\\
1&0&\cdots&0&0\\
\vdots&\ddots&\ddots&\ddots&\vdots\\
0&0&\ddots&0&0\\
0&0&\cdots&1&0
\end{bmatrix}

今回は上下左右のマスの値の合計を求めたいので、$V$を上下左右にシフトさせて和をとっていけばいいことになります。最外周は後段で値を書き換えるので1マスずらす分なら何も考慮しなくていいです。

C=S_2VS_1+S_2V+S_2VS_2\\
+VS_1+VS_2\\
S_1VS_1+S_1V+S_1VS_2\\
=(S_1+I+S_2)V(S_1+I+S_2)-V\\
=SVS-V

$S=(S_1+I+S_2)$です。これで一発で周囲8マスの値の合計が取れました。元記事の方ではGPU計算や他の言語への移植も考慮に入れてarrayfunctionを使用していましたが、MATLABであればこれが最もローコストな演算だと思います。

その他紹介したTips

描写について

もう一つ細かいTipsですが、ループでFigureを更新して動画やアニメーションを保存したり、Figureの更新で変化の仕方をリアルタイムで観察したい場合は、plotやscatterなどの描画関数を逐一呼び出すよりも、グラフのデータバッファを直接書き換えた方が更新が速いです。
plotなら

%% 初期値
[x,y] = % 何かしらのデータ;
p = plot(x,y);

while 1 % forでもいい

% xとyの更新

% 描画更新
p.XData = x;
p.YData = y;
drawnow

end

という流れです。

論理値配列を用いたインデックス

MATLABを使う方でもあまり馴染みがない使い方かもしれませんが、コードのシンプル化に対して非常に有用です。
詳しくはMATLABのドキュメントを参照して欲しいのですが、よく使うのは以下のようなコード。

A = % 何かの配列;
B = % 何かの配列;

Thresh = 5;

% B(i,j) < 5の時、A(i,j) = 0にする。
idx = B < Thresh;

A(idx) = 0;

idxが論理インデックスで、中身はtrueとfalseです。B<Threshとなる要素がtrue、そのほかがfalseになっています。MATLABではこのidxをそのままAのインデックスにでき、記述がシンプルになります。
また、idx1, idx2,・・・と複数のインデックスのandやorをとることもでき、

idx = (idx1 & idx2) | idx3

のようにすることで、より細かいフィルタリングができます。ライブスクリプトで値を動的に選択したり、クリッピングなどに使えます。

ツールボックス

話を戻して周囲8マスカウントですが、Image Processing Toolboxのimfilterを用いることで同じこと、いや、もっと汎用性の高い使い方ができます。今回の例なら

C = imfilter(V,ones(3,3));

で終了です。しかもこのimfilterは周囲8マスだけでなく、もっと広範囲の和も取れますし、各マスに対して重み付けすることもできます。Image Processingと銘打ってますが、波動の時間離散解析に使える(計算中に畳み込みフィルタリングが登場する)こともあって、音、光、電磁気、熱などの解析屋さんには重宝しそうです。

GPUArray

この畳み込みフィルタリングですが、上にあげたような物理シミュレーションや、ディープラーニングのニューロン間の畳み込みに応用できそうです。しかし、両者ともかなりのメモリを食う処理になりそうですね。そうなると当然計算速度は遅くなりますから、MATLABのPrallel Computing ToolboxのgpuArray()の紹介をしました。

この関数は非常に便利で、入力として配列を与えるだけでGPU上で勝手に並列演算してくれます。これだけで、普段ならめんどくさいGPUプログラミングができます。ただし、使用する際はメモリのリフレッシュなどしっかりケアしましょう。

最近は物理解析でもGPU演算が主流になってきた印象があるので、MATLABでも有効活用しましょう。

ソースコード

最後に、元記事のコードをアレンジしたものを載せておきます。

clear
%%
% 人口の平方根
PopulX = 300;

% 感染期間
Period = 10;

% 感染確率
Prop = 1/30;

% 個々人の感染状況
% 0->未感染
% 1~period -> ウイルス保有中
% period -> 回復済
State = zeros(PopulX);

% 初期感染者のインデックス
StateInitIdx = floor(PopulX/2);

% 初期感染者はパラメータ1(感染1日目)
State(StateInitIdx,StateInitIdx) = 1;

% 日数
DateNum = 1;

%% ウィンドウ
f = figure('Position',[100,300,1500,500]);
tile = tiledlayout(1,3);
ftitle = title(tile, [num2str(DateNum) '日目']);

nexttile
gNowInfect = imagesc(zeros(PopulX));
axis equal
axis tight
title('現在の感染者')

nexttile
gTotalInfect = imagesc(zeros(PopulX));
axis equal
axis tight
caxis([0 Period])
title('既感染者')

nexttile
pNowInfect = plot(0,1);
title('現在の感染者数')
axis square

%% ループ処理
% 【条件】
% ・param = 1~periodの人は周囲のparam = 0の人に一定確率で感染させる。
% 処理の簡単さのため最外周の人は感染しない。
% 初期段階で収束しないように最初の1人は回復しない。
% ツールボックスを使用しない。

%% 事前準備
% 方式切り替え
Method = 3;

% シフト配列
E = eye(PopulX);
Shift = circshift(E,1,1) + E + circshift(E,-1,1);

%% ループ
% for dateNum = 1:100
while 1
    if strcmp(get(f,'currentcharacter'),'q')
        close(f)
        clc
        break
    end

    % 感染者論理インデックス
    InfectIdx = State > 0 & State < Period;
    
    % 感染者数
    InfectNum = length(InfectIdx(InfectIdx == 1));
    
    % 描画更新
    ftitle.String = [num2str(DateNum) '日目'];
    gNowInfect.CData = InfectIdx;
    gTotalInfect.CData = State;
    pNowInfect.XData = [pNowInfect.XData DateNum];
    pNowInfect.YData = [pNowInfect.YData InfectNum];
    drawnow
    
    %% 周囲8マスにいる感染者の数カウント
    % 感染者の論理インデックスbioIdxから、各マスの周囲8マスにいる感染者数をカウントする。
    
    % 行列で対処。case 1を1行に収めたもの。
    % 人口x2分のメモリでOK
    AroundInfectSum = Shift*InfectIdx*Shift - InfectIdx;
    
    %% 
    % 非感染者インデックス
    HealthyIdx = State == 0;
    
    % ランダム確率行列
    SpreadPIdx = (1-Prop).^AroundInfectSum < rand(PopulX);
    
    % 次回感染者
    NextInfectIdx = SpreadPIdx & HealthyIdx;
    
    InfectIdx = InfectIdx | NextInfectIdx;
    InfectIdx([1,end],:) = false;
    InfectIdx(:,[1,end]) = false;
    
    % 感染者の日数をインクリメント
    State(InfectIdx) = State(InfectIdx)+1;
    State(StateInitIdx,StateInitIdx) = 1;
    
    DateNum = DateNum + 1;
end
11
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
11
9

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?