ウイルス進化のシミュレーションプログラムの説明(ソースコード無料配布)
はじめに
本記事は、C言語で書かれたウイルスの感染と進化の確率的なシミュレーションを行うプログラム「COVID19_VR.exe」の解説です。
・実効再生産数(RT)や人口、感染確率、感染者の体質、マスク人口比、マスクによる感染・被感染確率、自粛する人間比率、緊急事態宣言の発動条件、クラスターの数、クラスター間感染率、などを設定できます
・与えられた条件下で指定日数(デフォルトでは4年間)シミュレーションを行う（世界線を一通り経過させる）
・世界線を指定回数（デフォルトでは100回）通し、特定事象（変異ウイルスのパンデミックや医療崩壊）が起こる回数をカウントする
・変異ウイルスが出現する確率、変異ウイルスが市中感染に及ぶ確率を表計算形式（CSV）で出力
・特定の世界線における感染を度数分布グラフで出力
などが可能です。
記事の文字数は12万8400文字ほどですが、ほとんどソースコードの文字数です。実質的な文字数は4000文字程度なのでびびらないでください。
沿革
新型コロナウイルスが中華人民共和国の武漢で発生した最初の記録は2019年12月8日です。この記事を書いている2021/02/20においても未曾有の世界的な感染が今も続いています。
しかしついに、ワクチンという最大の武器がついに完成しました。ワクチン完成に尽力された研究者の皆様には、最大限の賛辞と感謝を申し上げたいと思います。ワクチンがもたらされた今、コロナ対策は次の段階を考えていく段階に入ったと思います。
・なぜ新型コロナウイルスはこれほどまでに感染性が高いのか
・なぜ新型コロナウイルスはこれほどすばやく進化するのか
・今後、同様の感染症を抑えるには何が必要か
新型コロナウイルスの感染力の要素については別記事
「新型コロナウイルスはなぜ感染力が高いのか考察」
https://note.com/yunyun_ja/n/nf5c1773d7157
をご覧ください。
新型コロナウイルスの進化についての考察と仮説については別記事
「新型コロナウイルスはなぜ進化が早いのか仮説＆考察」
https://note.com/yunyun_ja/n/n66c54cf909ee
をご覧ください。
本プログラムでわかったことと考察は別記事
「シミュレーションでわかったこと」
https://note.com/yunyun_ja/n/nfd9f0317ed9e
をご覧ください。
本記事では、ウイルス進化をプチ変異を仮定して複数回シミュレーションするプログラムの解説を行います。ソースコードも記事内にあるので、使用方法を読んで使ってみてください。
シミュレーションの目的は、別記事
「今後、同様の感染症を抑えるには何が必要か」
https://note.com/yunyun_ja/n/n69781a80ea82
なども含めて検証するためでもあります。本プログラムは、様々な統計データを設定し、感染者数や感染グラフなどの統計情報を得る簡単なシミュレーションプログラムです。C言語で書かれ、移植性を高めるためグラフィック等はなく、出力は全てテキストデータです。故に、適切なCコンパイラがあればLinux、Windows、Macいずれでもほぼ何もいじることなくコンパイルし、実行できるはずです。
最初は単純な統計的仮定をおくことで、極めて高速に感染シミュレーションが行えることから1000行～2000行程度のシミュレーションプログラムのつもりでしたが、実用的になるには結局体内ウイルス量管理や、治癒・重症化・致死ステージなども必要であるとわかり、コメント含めて4000行超のプログラムとなってしまいました。
バイナリ配布も考えましたが、昨今の生き馬の目を抜くようなネット上で開発者不詳のプログラムなど怖くて使えないであろうという点、そしてなにより、興味のある人にはしっかりいじって欲しいし、より興味深いデータを引き出してくれる優秀な人にいじってもらえれば、という思いから、ソースコードを公開します。
なお、これは失敗なのですが、初心者でもコンパイルできるようにと「すべて単一ソースコード」にしています。ヘッダファイルにするべき内容は先頭にまとめてはいますが、べた書きです。さらに本当は4種類のソースコードに分けているのですが、掲載に際してまとめています。
出典について
先に書いておきますが、公開時に出典を明記してくだされば、本プログラムは改造、公開自由です。ただし、僕は他の記事を見てもらえばわかるとおり不正を嫌います。本プログラムは大学生のしょぼい卒論程度にはなりうると自負しており、その点から卒論に関して真面目にやっている学生をこのプログラムを転用することで出し抜くことがない様、一応次の規定を設けます。不正が嫌いなだけなので、その文脈においての規定です。不正と同じくらい余計な束縛も嫌いな性分なので。
・無断転載は公開日数×5万円申し受けます
・本プログラムで得られた結果を根拠としてワクチンや感染症対策の邪魔をした人からは、著作者人格権の侵害として、迷惑をかけた（機関数OR部門数）のうち多い方をaとして、a×20万円を申し受けます
・卒論・修論・博論・その他学習機関におけるレポートに、本プログラムを主とした研究と称した内容を記載した方からは、20万円申し受けます
こんなことは書きたくないんですが、書かないと理解できない人も多いのでごめんなさい。
コンパイル方法
標準関数しか使っていないので、Cコンパイラならどれでも使えると思います。ただし、ウイルス数で兆単位の数値を使うので、16bitOSは流石に無理でしょう。32bitのWindows7なら余裕です。とりあえず、WindowsでMingw等のgcc系コンパイラが入っている場合を想定します。
①．本プログラムのソースコードを特定のフォルダ・ディレクトリに文字コード「UTF－8」で保存します。ここでは、「Windows」で「D:￥work」に保存したとします。ファイル名は「COVID19_VR_tenkai.c」です。
②．「D:￥work」でコンパイルコマンドを打ち込みます。日本語の書かれたCSVファイルを出力するためUNICODE対応でコンパイルします。ここでは、「gcc」を使う場合のコマンドを例示します。
gcc -DUNICODE -D_UNICODE COVID19_VR_tenkai.c
ファイル名を指定する場合は、
gcc -DUNICODE -D_UNICODE -o COVID19_VR.exe COVID19_VR_tenkai.c
です。
特殊なフラグ
コンパイル時に次のマクロを与えるとコンパイル設定が変化し、プログラムの挙動が変わります。
ShutdownEnd
シミュレーションが終わったら、システムをシャットダウンします。シミュレーションに数時間かかる場合、寝ている間にやってもらうことになりますが、つけっぱなしはいや、という場合はこのマクロを有効にしてコンパイルしてください。ただし、シャットダウンコマンドはOSごとに異なるので「int myPCshutdown(void)」という関数処理を書き換える必要があるかもしれません。コマンドは
gcc -DUNICODE -D_UNICODE -DShutdownEnd -o COVID19_VR.exe COVID19_VR_tenkai_ShutdownEnd.c
です。
MERSENNE_TWISTER
乱数生成がデフォルトからメルセンヌツイスタに切り替わります。
http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT/index-jp.html
僕が、乱数の性質がアルゴリズムに影響しないか調べるために定義したものです（デフォルトがgccの線形合同法の場合の挙動は変化ありませんでした）。本プログラムのアルゴリズムにおける乱数の使用方法であれば、線形合同法でもメルセンヌツイスタでも影響がないことは確認済みです。このマクロは外部ファイルを要求するので使用しないでください。使ってもリンクに失敗します。
Details
デバッグのために、一回の感染ごとに感染者の体内のウイルス数まで全てログとして記録します。あっというまにGB単位でHDを圧迫する上、ありえないほどパフォーマンスが低下します。オススメしません、というか開発者以外は使いません。デフォルトはオフなのでそのままで大丈夫です。
プログラムの使い方
コンパイルが終了すれば、「a.exe」または「COVID19_VR.exe」ができていますので、これを実行します。実行すると、デフォルトでは
Simulation_Run5_ex3( fp_log, fp_glp, seed );
が実行され、4系統のファイルが生成されていきます。この関数はクラスター間の感染率を変化させて、医療崩壊が起こる確率を算出します。
出力結果はテキストファイルで、それぞれ
統計情報_世界000線○○人○○クラスタ(開始時刻).log
世界線を指定回数（デフォルトでは100回、000線～099線）繰り返したときの各世界の統計情報です。世界線数だけ生成されます。人口当たりの免疫特性などは毎回同じですが、緊急事態宣言の発効と解除は毎回異なります。同じシミュレーションを検証するための、乱数の種もここに記録されています。
度数分布_世界000線○○人○○クラスタ(開始時刻).log
一目見てわかる、テキスト形式の感染者数の度数分布です。世界線数だけ生成されます。原種タイプのコロナウイルスを1、プチ変異タイプを2,3、プチ変異が遺伝的に混じって超変異した株を4としています。もっと色々な変異ウイルスを作りたければ、マクロ定義
#define VIRUS_N (3+1)
#define VIRUS_0 -1
etc.
と関連する箇所をいじれば一応作成可能です。でも面倒ですよ。
度数分布は各日付ごとに最も多いタイプを左に寄せます。Simulation_Run5_ex3関数では、市中に突然デルタ株のような実効再生産数の高い変異株が集団内に4株現れたときにどのように広がるかをシミュレーションしています。デフォルトでは、原種の1.8倍の実効再生産数でシミュレートしています。
複数回シミュレーション結果(開始時刻).log
世界線ごとの特定事象を神視点で記録したファイルです。つまり、各条件ごとに100回ずつ世界を作り直して試し、世界が滅んだ回数を記録するようなそんな視点です。デフォルトでは医療崩壊が起きた回数を、条件ごとに記録します。このファイルをCSV形式に貼り付けるかCSVとして開くと、グラフを作成できる表計算データになります。
複数回シミュレーショングラフ用データ(開始時刻).csv
世界線ごとの特定事象を神視点で記録したファイルです。Simulation_Run5_ex3関数ではあまり意味を持ちません。もっと詳細なデータが重要なシミュレーションで、各世界線を俯瞰して記録します。
といったファイルが生成されます。
各関数の取り扱いは別稿で。
バッチファイル
シミュレーションには時間がかかる場合があります。特に、人口が増えると感染処理や病状の進行処理、ウイルス増加処理などが増えるので数時間～数十時間かかるかもしれません。そういった場合は、家でPCを使っているときじゃなく、寝ている間や仕事に行っている間にシミュレーションしてほしいですよね？
そういう時は、古きよきバッチファイル（笑）が地味にオススメです。
あらかじめフォルダごとに別々の設定でコンパイルした実行ファイルを入れておき、バッチファイルで実行します。こうすると、後で整理するときに便利です。
最後に実行するシミュレーションプログラムを「ShutdownEnd」フラグ付でコンパイルしておけば、最後のタスク終了後にPCがシャットダウンされます。
ソースコード
コイツをUTF-8形式のテキストファイルとして保存してください。この説明記事に則れば、ファイル名は「COVID19_VR_tenkai.c」です。
#include <tchar.h>
#include <stdio.h>
#include <time.h>
#include <sys/time.h>
#include <conio.h>
#include <windows.h>
#include <stdio.h>
#include <wchar.h>
#include <locale.h>
// 2021/01/16 新型コロナウイルスSARS-CoV-2による感染症COVID-19 に関する
// 数値シミュレーションを行うため開発開始。
// 本プログラムの目的は、コロナウイルスが変異株を生み出す条件を確率的に
// シミュレートし、変異が”広まる”確率をなるべく下げる手法を考察する事である。
// 僕の考えでは、コロナウイルスは極めて厄介な形質を既に獲得している。
// これまでの感染予防対策では、コロナウイルスを押さえ込めないかもしれない。
// 有効な手段を早期に見極める必要があるのだ。
// 当初は無向グラフを使った感染シミュレーションを考えていたが、
// 感染と変異のみに着目するならば、グラフによる固定的な人間関係よりも、
// むしろ大量の配列データに対するランダム感染の方が実態に沿っていることに
// ランニング中に気付いた。コレならばかなり計算量が少なくてすむ。
// C言語であれば、万単位の人間データもPC処理能力で計算できるだろう。
//
// 本プログラムにおける仮定・用語
// ・ウイルスの進化には、形質として現れている変異の前に、隠れ変異が起きていることが多い
// ・マスクはウイルスにとって物理的で無慈悲でほぼ無差別なボトルネックである
// ・マスクや隔離によって感染関与ウイルスの数が減れば、大数の法則が働かず適応度を無視した絶滅が起こりうる
// ・遺伝子プールとなる密なクラスターは進化を促す。とりわけ、コロナはクラスターの密度が濃い
// ・感染が確立した患者同士であっても、ウイルスの相互交換はなるべく避けるべき。
// ・進化速度を鈍化するには、クラスターの密度を下げるのが有効かも
// ・感染において少なくとも四つのボトルネック、障壁があり、障壁に適応度があればそのウイルスは個体数を増やす
// ・四つの障壁とは、「体内における増殖競争」「体外排出ルート」
// 「対人空間（距離や風なども含む）」「クラスター跳躍」
// ・実効再生産数Rtは、ある条件下で、一人の感染者が何人の感染者を新たに生み出すかという大まかな期待値と
// いうことができる。その条件は無数にあり、t、つまり感染が始まってからの時間も当然関係する。
// 実効再生産数は、実務的にはデータから導き出す理論値であり、未来に対して常に不確かである。
// 本プログラムでは、実効再生産数はREPTNUM_VIRUS というマクロで定義している。
// この値は、「ウイルスにとっての、感染の際のチャレンジ引換券」のような扱いであり、基本的には
// この値分だけ人間に感染を仕掛けることができる。
// 2.4なら、2.4人に対して感染する(実際には確率的な2.4人。もっと厳密には、成功率50%の感染を4回、成功率40%の感染1回試行する)。
// この値に対して、さらにマスクや体質や体内ウイルス数が加味され、感染成立がシミュレートされる。
// この扱いは、どちらかというと「無防備な集団に対して感染するときのウイルス種固有の値」である
// 基本再生産数R0("あーる・のー"と読む。0は下付けの小さいゼロだが、ここでは表現できない)に近い。
// なぜ実効再生産数としているかというと、一言で言えば
// 「計算上はR0として組み込んでいるが、理想のRtを出すための調整値」であるからといえる。
// 本プログラムでは実効再生産数は感染パラメーターのひとつであり、重要なのはウイルス変異の動向である。
// Rtは事前のシミュレートで実際のコロナウイルスの動向を反映できる程度に「手動で調整」している。
// つまり、式としてはRtではないのだが、Rtとして機能するように値を調整しているのである。
//
//
// これらの条件を色々変更し、対策を練ることとする
//
// 本プログラムにおける、新型コロナウイルスSARS-CoV-2による感染症COVID-19 の一般的性質
// ・感染してから発症までは5日、感染してから他者に感染すことができるようになるまで3日
// 軽症者は感染から治癒するまでに11日かかる
// ・感染後の死亡にいたるまでの日数はバラバラである。ただし、超重症だと発症から死亡までの平均日数は. 17.8日間
// 程度の目安はある。重症化した人からの感染リスクは、本人の能動的活動が不可能になることで大きく変化し、
// 十分に対策が取れていれば概ね0になるとしても問題ない。
// 本シミュレーションでは、人間の群れのサイズを減少させる要因としてのみ重要とする。
// 厳密に言えば、医療の行き届かない地域においてはより大きな意味を持つ。
// 医療の届かない人たち、例えば保険に加入できず自宅で療養するしかない、あるいは家すらない貧困者は、
// コミュニティ内で複数系統のウイルスに感染し、遺伝子を体内でミックスする
// 変異の試験場となってしまうという皮肉なリスクがある。社会的弱者を放っていては、いずれ手痛い教訓を得ることになるだろう
// 医療の行き届かない社会的弱者に対して、時に社会的強者は残酷にも切り捨てるような真似をするが、
// それは「人間側の遺伝子プールの減少」を招くため、長期的に見て圧倒的愚策である。
//
//
// 2021/02/05
// 現時点で理論的に考えられるウイルスへの対抗策は以下の要点を押さえるべきと考える
// ・分断・・・ウイルス同士の遺伝子交換がなるべく起こらない環境が望ましい。
// ・反グローバル・・・人体間の移動が容易になってきている。なるべくローカルに留めるべき。
// ・不条理な差別、人事・・・能力ではなく偶然や生まれに左右されるような選別をするべし。マスクは神。
// ・体外に出るウイルス量自体を減らすべきなので、既に発症している人に咳止めを処方するのは良いかも
// ・多様性を潰す。現在無能に見える変異も長期的には集団に貢献する可能性がある。
// 2021/03/08
// 夢の中で重要な要素に気付いた。変異株が、変異の母体となった人間の体内で少数派であるうちは、他者感染のリスクは高くない。
// そして、感染リスクは体内ウイルス量に依存するとすると、ウイルスの増殖ステージが後半になるほど、
// 変異株にもそれなりの感染チャンスが生まれることになる。
// つまり、発症の後期においても症状があまり出ない不顕性感染をした人物が、体内でウイルス変異を起こしていると
// これは変異株のウイルス集団にとって僥倖である。
//
// 2021/03/18
// int fprintf_daily_statistics2()等の統計関数の集計処理をget_statistics_flock()に集約
// 当日の感染者数と累計感染者数の違いから、累計の方は関数内で独自に集計処理をしていたが、
// 統計データを集約する構造体STATISTICSの拡充により全てまとめられるようになった。
// 2021/04/11
// 集団をクラスタに分けて、クラスタ間感染率を変化させて実験データを取ったところ、まったく変化がないことが判明
// 恣意的に値を求めるのは間違っているが、さすがにこれはおかしいので詳細を出力したところ原因おそらく判明
// 理論的な値は後述するが、緊急事態宣言の条件が現実に即していなかったため、同時感染者数が最大で200人も許される事態
// であったのが最大要因と思われる。実際には歯止めをかけるには数日を要するのでもっと多いであろう。
// 2021/04/22
// いよいよ集団をクラスタに分けて行うシミュレーション。ところがクラスタを極限まで絞り、なおかつクラスタ間感染を
// 絞っているのに感染爆発が起きる問題が起きた。いくらなんでもおかしいので見直してみると
// 感染を行う daily_infection2() で法の処理を行う際に、
// roll_the_dice(PER_DEF, mod_wait)をして、判定がはずれだったときにRtを減じていないという
// 初歩的なミスだった。これでは、いくらダイスロールがはずれても感染が成功するまで
// 感染試行をしてしまう。
// daily_infection2_sub_Rt_Decrease_sars_v_n( sars_v );
// を追加して解決。このミスのために2日かけてこつこつ算出した880MBものデータがゴミ屑と化した・・・
//
//
//
//
//
// :::::::::::重要！！:::::::::::
// 本プログラムで境界バグが起こりやすいのは、体内ウイルス個体数に関わる計算である。
// 具体的には、SARS_V.Population の値である。ウイルス数は指数関数的に増大するので
// 値はかなり大きいものになる。一方で、格納する変数の大きさは十分とは限らない。
// 32bit環境ではintは2147483647(21億)、long long は9223372036854775807(922京)だろうか？
// これらは環境依存なので事前には確定できない。
// 体内ウイルス個体数は、ウイルスの増殖時および変異ウイルスの発生判定に用いられる。
// ・ウイルスの増殖時は整数オーバーフローに注意する必要がある。
// ・変異確率の計算時には、大きすぎる値との計算で値が切り落とされないかチェックすべきである。
// 特に変異の計算値は目に見えにくい稀な現象であるため、計算にバグがあっても発覚しにくいのである。
//
// :::::::::::改善！！:::::::::::
// パフォーマンス上のボトルネックについて
// 現在(2021/03/26)に把握している重い関数はget_statistics_flock()である。
// 人口3000人、4年間の世界線を20回シミュレートすると実に実行時間の50%近くを占め、58238回呼び出されている。
// 次点で重いのはdayAdvanceだが、これはメモリのコピーを行っているだけで、アルゴリズムの余地はない。
// 詳細は20210326_082624_COVID19_VR.prof.log、プロファイル生データは20210326_082618_gmon.out。
// 改善済み
// 調べてみると、緊急事態宣言を発令・解除するdaily_event()が毎回患者数のカウントを行うために使っていることがわかった。
// 本来であれば統計情報はdaily_statistics_calc()が再計算が必要なときだけget_statistics_flock()を呼び出すように
// 考えていたが、データの更新は一日の中でもウイルス変異や増殖など刻一刻と変化している。
// つまり、再計算が必要かどうかの管理は意外と粒度が細かく、ちょっと厄介な問題を抱えている。特に、
// 「市民のデータに変更があったら必ず再計算フラグを立てる」必要があり、これはかかわりのある全関数に徹底させる必要があるが、同時に
// 「市民データに統計的な変更がないなら再計算フラグを立てない」必要もある。
// でないと、変化を生じさせうる関数が呼び出されただけで変更ありとみなされかねない。
// そうなれば、結局関数呼び出しのたびに再計算が必要とみなされるのでまったく改善にならない。
// さらに、変化系関数daily_infection2()、daily_mutation2()、daily_proliferation()などはそれぞれ
// 与える変化も異なるのでチェックの要諦も異なるのである。どのループ内でどの処理が行われたら変化があったとみなすか。
// シミュレーションでデータの取得ミスは致命的なので、このリスクは高すぎる。
// 一見余計な実装であるが、以上の理由から、daily_event()で呼び出される、緊急事態宣言用の2つの関数
// daily_event_sub_EmergencyOnCheck()、daily_event_sub_EmergencyOffCheck()が使うため、
// 統計用関数とは別の患者数のみカウント用関数、daily_event_sub_PatientCount()を別個実装する。
//
//
//
// 現時点でわかっている、感染を抑える方策一覧
//
// 集団の人数を減らしていく
// クラスターあたりの人数を減らしていく
// 異なるクラスター同士の接触（感染確率）を減らす
// 自粛する人数の割合を増やしていく
// 緊急事態宣言を早めに出す
// マスク装着率を増やす
// 現時点でわかっている、変異ウイルスの発生率を抑える方策一覧
//
// 医療崩壊を防ぐ
// 集団の人数を減らしていく
// クラスターあたりの人数を減らしていく
// 異なるクラスター同士の接触（感染確率）を減らす
// 自粛する人数の割合を増やしていく
// 緊急事態宣言を早めに出す
// マスク装着率を増やす
// 乱数生成は基本的にC標準のrand()とsrand(seed)を使う。
// 本プログラムは誰にでも試せることを旨とするので、外部ライブラリを意図的に排している。
// しかし、シミュレーションにおいて乱数の質は大事なので、マクロによる切り替えで
// メルセンヌツイスタも使えるようにしている。
//
#ifdef MERSENNE_TWISTER
#include "rand.h"
#define irand sfmt_rand_n
#define roll_the_dice roll_the_dice_mt
#define srand sfmt_srand
#endif
//ファイル名の最大値
#define FNAME_MAX (256*2)
// クラスターに含まれる人間数のデフォルト値
#define MAN_N 100
// この実効再生産数については注意が必要。感染症対策において使わる実効再生産数Rtというよりは
// 考え方としては基本再生産数R0に近い。しかし、この値はシミュレーションの条件に合わせて「調整される」
// 意味合いが強く、その意味で実効再生産数である。詳しくは本ファイル上部の「用語」で解説。
// 通常、日本におけるコロナウイルスの実効再生産数は3.4程度と目されているが、本シミュレーションでは
// マスクあり80%、自粛する人間80%、集団10万人クラスタサイズ10000としたときに3前後のRtを示す値として
// として調整している
// 目安として、最初の感染者が1人で実効再生産数2なら
// 初期値 5日目 10日目 15日目 20日目 25日目 30日目
// 当日感染者数 1 2 4 8 16 32 64
// 感染中の人数 1 3 6 12 24 48 96
// 累計感染者数 1 3 7 15 31 63 127
// 4人で実効再生産数3.4なら
// 初期値 5日目 10日目 15日目 20日目 25日目 30日目
// 当日感染者数 4 13.6 46.2 157.2 534.5 1817.4 6179.2
// 感染中の人数 4 17.6 59.8 203.4 691.7 2351.9 7996.6
// 累計感染者数 4 17.6 63.8 221.0 755.5 2573.0 8752.2
//
#define REPTNUM_VIRUS_1 ( 12.0 ) //実効再生産数 Effective reproduction number のデフォルト値新型コロナは通常2.4
#define REPTNUM_VIRUS_4 ( 1.8 * REPTNUM_VIRUS_1 ) //実効再生産数 Effective reproduction number の変異型の値
#define ONEDAY_TICKET (REPTNUM_VIRUS_4*10/9/10)
// 感染 確率probability 粒度 particle
// 人人感染を起こすときの確率試行の細かさ。
#define PROBABILITY_PARTICLE (0.5)
//1日に突然変異が起こる確率の逆数。100なら1日にウイルス1個体あたり100分の1の確率。
#define MUTATION_PROBABILITY 10000
// 感染者の体内で変異ウイルスが発生するタイムリミット。
// これがないと、治癒しない患者の体内でいずれ変異ウイルスが発生する
#define MUTATION_DAY_LIMIT 14
//1日に増殖する倍率
#define PROLIFERATION 1000
//体内に存在できるウイルス数の限界値。とりあえず1兆
#define VIRUS_MAX 1000000000000
// 本人の遺伝的形質
#define CONS_STRONG 1 // STRONG 感染しない無敵個体
#define CONS_HIDE 2 // HIDE 感染しても症状がないが、5日後に他人に感染すことが出来る隠れ感染者
#define CONS_MIDDLE 3 // MIDDLE 普通に発症する体質。感染すると発症し、5日後に他人に感染すことが出来る
#define CONS_SEVERE 4 // Severe 重症化する体質。発症から1週間程度で重症化
#define CONS_DEATH 5 // DEATH 重症化後死亡する。深刻SERIOUSな体質。
// 本人の状態
#define CONDITION_ALIVE 1 // 健康に動いている
#define CONDITION_INCUB 2 // 発症している
#define CONDITION_SEVEE 3 // 重症
#define CONDITION_IMMOV 4 // 不動 Immovable
#define CONDITION_DEATH 5 // 死亡
// 本人の行動特性(性格)。性格によっては、感染チャンスを複数回獲得する。つまりRtに強化値がかかる
#define PERSONALITY_INFLUENCER 1
#define PERSONALITY_ALONE 2
#define PERSONALITY_MIDDLE 3
// この値は注意。本シミュレーションでは、「感染」に重点を置いているので、基本的な多くのイベントは
// 「感染日」が起算となる。
// ただし、「他者へ感染可能になるタイミング」は、体内ウイルス量が重要であり、排出するウイルス量も
// 体内ウイルス量にある程度依存するべきである。よって、以下のタイミングは次のようにして決められる
//
// 上記を補足すると、TIME_LIMIT_CAN_INFECTは計算に用いるだけで実際に何らかの条件になるわけではない。
#define TIME_LIMIT_CAN_INFECT 3 // 他者へ感染可能になるタイミング。感染後3日
#define TIME_LIMIT_INCUBATION ( TIME_LIMIT_CAN_INFECT + 2 ) // Incubation period 潜伏期間。発症するタイミング
#define TIME_LIMIT_RECOVERY ( TIME_LIMIT_INCUBATION + 7 ) // 普通に発症し重症化しない遺伝的要因をもつ人が治癒するタイミング
#define TIME_LIMIT_SEVERE ( TIME_LIMIT_INCUBATION + 7 ) // 重症化する遺伝的要因をもつ人が重症化するタイミング
#define TIME_LIMIT_DEATH ( TIME_LIMIT_INCUBATION + 17 ) // 重症化した人が死亡するまでのタイムリミット。発症後17日程度とする
// 発症スケジュール(感染した日は余裕を持つため0起算、それ以外は1起算)
// 一般的に、多くの感染症は潜伏期間(症状が出るまでの期間)と他者へ感染すことができない期間が一致するが
// 新型コロナウイルスは、潜伏期間である無症状期が終わる2日ほど前から他者へ感染すことができる。
// このずれが、活発に動くバイオハザードな2日間という悪夢を生む。本人は無自覚である。
//
// 軽症の場合(感染後5日で発症、その2日前から他者へ感染は可能、発症後7日で治癒、発症後14日までは他者へ感染可能)
// ┏━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━
// ┃ 無症状┃発症 ┃治癒
// ┣潜伏期間 ┏━━━┻━━━━━━━━━━━━━┻━━━━━━━━━━━━━┳━━━━━━
// ┃ ┃感染可 ┃感染の恐無し
// ┣━┳━┳━╋━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━╋━┳━┳━┳
// ┃ 0┃ 1┃ 2┃ 3┃ 4┃ 5┃ 6┃ 7┃ 8┃ 9┃10┃11┃12┃13┃14┃15┃16┃17┃18┃19┃20┃21┃
// ┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃
// ┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃
// ┗━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻
// 重症の場合(感染後5日で発症、その2日前から他者へ感染は可能、発症後7日で重症、予後が悪いと発症後17日で死亡)
// ┏━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━
// ┃ 無症状┃発症 ┃重症化 ┃死亡
// ┣潜伏期間 ┏━━━┻━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━┻━━
// ┃ ┃感染可 強い感染力
// ┣━┳━┳━╋━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳━┳
// ┃ 0┃ 1┃ 2┃ 3┃ 4┃ 5┃ 6┃ 7┃ 8┃ 9┃10┃11┃12┃13┃14┃15┃16┃17┃18┃19┃20┃21┃
// ┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃日┃
// ┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃目┃
// ┗━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻━┻
//マスクをつけるかどうか(HUMANS構造体用)
#define MASK_ON 1
#define MASK_OFF 0
#define MASK_DEFENSE (0.7) // 感染させられる側がマスクをつけることによるウイルス減少率
#define MASK_WEAKNESS (0.7) // 感染させる 側がマスクをつけることによるウイルス減少率
//自粛するかどうか(HUMANS構造体用)
#define SELF_RESTRAINT_ON 1
#define SELF_RESTRAINT_OFF 0
// 免疫を獲得しているかどうか(HUMANS構造体用)
#define IMMUNITY_ON 1
#define IMMUNITY_OFF 0
// コロナウイルス感染において、次のウイルスが増殖してくるには約6時間という結果が出ている。
// 仮に100倍になるとすると6時間後に100個、半日で1万、一日で100000000個(1億)個である。
// https://www.jmedj.co.jp/journal/paper/detail.php?id=14354
// SARSコロナウイルスは，6時間程度で増殖し，10^5〜6／mL程度のウイルスを産生する
// https://www.jmedj.co.jp/journal/paper/detail.php?id=14278
// コロナウイルスの体内ウイルス数についてはきちんとした調査結果があまりないため、
// 次のデータはノロウイルスのデータを最悪パターンとして一時的に採用する。
// ・10〜１００個ほどで感染が起こる
// ・便１グラム中には、１００万〜１兆個のウイルス
typedef struct SARS_V{// SARS-CoV-2のプロパティ
int type;// ウイルス種別従来型は1
long long Population;// 個体数 Population オーダーは×POPULATION_ORDER (たぶん100)。つまり、Populationが1なら100個体として扱い、感染成立。
double Rt;//実効再生産数
}SARS_V;
// 感染するときのウイルス個体数。対象が無防備であることを想定したデフォルト値。マスク等で減る。
#define Default_VIRUS_INVASION_N 10
// 感染したウイルス個体数に下駄を履かせる値。Populationにこの値をかける。つまりPopulationが1なら100個体。ダースみたいな基数単位
#define POPULATION_ORDER 100
// 人人感染を起こすための体内ウイルス個体数の最低値。これ以下だと人人感染できない。
// 潜伏期間を5日とするなら、この値はPROLIFERATION^4 が妥当。1000なら1000*1000*1000*1000
#define INVASION_LIMIT (long long)( (long long)(PROLIFERATION) * (long long)(PROLIFERATION) * (long long)(PROLIFERATION) * (long long)(PROLIFERATION) -1 )
#define INVASION_LIMIT_FOOT_CUT (INVASION_LIMIT/1000000) //daily_infection_sub_virus_invasion()でダイスロールできる値まで足きりする値
//#define INVASION_LIMIT (PROLIFERATION * PROLIFERATION * PROLIFERATION -1 )
// 人人感染を起こすとき、体内ウイルス量はINVASION_LIMIT。他者感染はここからある程度絞られて起こる。
// この値で大事なのは、従来型がINVASION_LIMITまで増殖したとき、変異型が少量ながらも感染可能性をもつ値であることである。
// この値はdaily_infection_sub_virus_invasion( )が感染判定するときに使う。この値とINVASION_LIMITの比で
// 感染ウイルス量が算出される。最初の変異型が現れたときに変異型が従来型の1/10でも感染可能とするなら
// INVASION_LIMIT ÷ INVASION_LIMIT_RATIO = 10とする
#define INVASION_LIMIT_RATIO (INVASION_LIMIT/10)
// 本プログラムで扱うウイルスの変異種数+1（親株）
#define VIRUS_N (3+1)
// 本プログラムで扱うウイルスの変異
// 単純な数値だけど、後で変更するかもしれないので一元管理しとく
#define VIRUS_0 -1 //デフォルト値。この値のときはまだウイルスが登録されていない(この枠に未感染)
#define VIRUS_1 1 //従来型。
#define VIRUS_2 2 //潜在的な変異1。まだ特に有利な形質ではない
#define VIRUS_3 3 //潜在的な変異2。まだ特に有利な形質ではない
#define VIRUS_4 4 //危険な変異型。変異1と変異2が同時感染した人間内でRNA交換を通して発生する。
#define SUM_POPULATION( a ) (a.sars_v[0].Population + a.sars_v[1].Population + a.sars_v[2].Population + a.sars_v[3].Population)
typedef struct HUMANS{// 人間のプロパティ
int no;// 個人ID
int cons;// 体質 Constitution 感染後の予後を決める
int count;//感染後の経過日数。感染した段階で0になる。通常は-1。ウイルスごとにカウントしないのは、免疫ができるまでの日数は初期感染日時に依存すると考えられるため
SARS_V sars_v[ VIRUS_N ];
int virusn;//何種類に感染しているか
int condition;// 元気か、動けるか、生きているかどうか
int mask;//マスクをつけているか
int self_restraint;//行動自粛できるか
int personality;//インフルエンサー等行動特性(性格)
int immunity;//免疫があるか
}HUMANS;
typedef struct FLOCK{// 人の群れ(FLOCK)情報を保持する。元々はGroupを使いたかったが、winsock2.hと衝突。
HUMANS *humans;
int n;//人口population。HUMAN構造体配列の個数。
}FLOCK;
// 緊急事態宣言について、その数値スケールに問題が発覚したのでその理論を修正する。
// 元は、ACT_SELF_RESTRAINT_LIMIT を 2とし、100分率スケールで扱っていたのであるが、この値は緩すぎると考えられる。
// 東京都は、緊急事態宣言の解除条件として、「10万人あたりの新規感染者が25人を下回る」ことを目安にしている。
// 本プログラムでは、「ウイルス保持者数」に注目したシミュレーションを行う。
// そこで、10万人あたり25人をウイルス保持者数に換算する。
// 新規感染者が25人になった日をn日とする。このとき、n日までにも感染者が出ていて、
// その感染者数は25以上ということになる。感染者は快癒する場合は感染から14日程度で快癒する。
// ということは、n日までの感染者数の総数はn-14日までの感染者をカウントすることになる。
// その総数は25×14=350より多いと考えられる。例えば
// n-14日目 81
// n-13日目 74
// n-12日目 68
// n-11日目 62
// n-10日目 57
// n-09日目 52
// n-08日目 48
// n-07日目 44
// n-06日目 40
// n-05日目 37
// n-04日目 34
// n-03日目 31
// n-02日目 29
// n-01日目 27
// n 日目 25
// とすると、合計は709となる（25から逆算して、=INT(C18*1.1)みたいに換算。減衰率は適当）。
// つまり、1.1減衰するとすると、10万人あたり15日間で709人程度の感染者と計算できる
// これは0.709%であり、1万人あたりでも70人までしか感染者を許さないかなり厳しい値である。
// 現実にはもっと厳しいだろう（減衰率1.1は少ない）から、70人より少なくても発動されるかもしれない。
// ちなみに、東京都は1000万人規模なので、7万人程度が感染した時点で緊急事態、
// 渋谷だけなら住人は20万人超なので1400人までである。
// よって、2021/04/14のアップデート以降、全て万分率に変更する。
// 移行時の処理もれがないよう、百分率が関わるメンバ変数名は末尾に_mを付加し、
// 全てのコンパイルが通ってから外すこととする。そして、分率計算には以降マクロ PER_DEF を用いる
#define PER_DEF 10000 //100分率、1000分率、万分率のような、特定の値に対する比で表現する場合の値
#define ACT_SELF_RESTRAINT_M 8000 //ACTIONS.self_restraint のデフォルト値。緊急事態宣言で自粛する人の割合。0-100の万分率
#define ACT_SELF_RESTRAINT_LIMIT_M 20 //ACTIONS.self_restraint_limit のデフォルト値。緊急事態宣言発動条件。0-10000の万分率
// 本構造体は、グループの細かいパラメータを設定するものである。
// 全員が均質な動きをするわけではないので、ランダムにこれらの値を設定する。
// 人群の行動傾向もウイルス拡散に大きく影響する。例えば、ウイルス感染における実効再生産数は
// 特定のグループの感染傾向から「推測する」ことになるが、当然別グループでは異なってくる。
// 布マスクをつければ、ウイルス自体の能力は変わらなくても、ウイルス元からは7割減、感染対象も1割減とされている。
// こういった行動の割合を一元管理するのがこの構造体である。この構造体で割り当てられた確率的値に基づいて
// 各人間のパラメータや感染時の確率変動を計算する。
// 基本的に10000分率を使い、整数値で0～10000の範囲で指定する。
// 端数が生じる場合は、「特定の方へ切り上げ」。
// 例えば、maskメンバの「マスクをつけるかどうか」は、マスクをつける方に切り上げる。
// 全体人数が10でこの値(ACTIONS.mask)が3000(万分率)なら、30%がマスクをつける。これは10人中3.3人なので4に切り上げ。
// consで始まる体質シリーズは、全体が10000になるように指定する。
// consで始まる体質シリーズは、全体が10000になるように指定する。
// self_restraint は特定条件で外出自粛をする人の割合。感染者の人口に占める割合がself_restraint_limit以上を条件とする
// personality_ で始まる性格シリーズは、平均的な人たちに対して行動が特殊なため実効再生産数が異なる人たちのパラメータ
typedef struct ACTIONS{//グループ全体の行動傾向について、まとめて指定するための構造体。
int mask_m;//マスクをつける確率。万分率で0～10000で指定する。
int cons_strong_m;//感染しない無敵個体の、全人口に対する割合(最大10000の万分率)
int cons_hide_m;//隠れ感染者の、全人口に対する割合(最大10000の万分率)。最も恐ろしい、無症状感染者
int cons_middle_m;//通常感染者の、全人口に対する割合(最大10000の万分率)
int cons_severe_m;//重症化する人間の、全人口に対する割合(最大10000の万分率)。感染後10日目に重症化し、感染させる能力は落ちる
int cons_death_m;//感染後14日目に死ぬ人間の、全人口に対する割合(最大10000の万分率)
int self_restraint_m;//行動を自粛する割合。全人口のself_restraint_limit%が感染したときに行動を自粛する割合(最大10000の万分率)
int self_restraint_limit_m;//行動を自粛する閾値。(最大10000の万分率)。全人口のどれくらいの割合が感染したときに自粛に入るか。
int personality_Influencer_m;//活動的な人物の割合。(最大10000の万分率)。感染においてRt×3。1ターン3回攻撃みたいな感じ
int personality_alone_m;//インドアな人物の割合。(最大10000の万分率)。感染においてRt×0.2。自身が感染してもほぼ他人に感染させない。偉い。
int personality_middle_m;//平均的人物の割合。(最大10000の万分率)。
int personality_Influencer_weight_m;//インフルエンサーの能力強化値。万分率で指定。通常の人を10000としてRt値を増やす。高いほど感染に貢献する。デフォルト30000
int personality_alone_weight_m;//インドアな人物の能力弱体化値。万分率で指定。Rt値を減ずる。デフォルト2000
int immunity_m;//免疫を持つ人の割合。万分率で指定。
int modulo_M;//クラスタ分けのための法。通常、合同算術では法nと表記するが、本ライブラリでは人数nなので法はmと表記する
int mod_risk_wait_m;//異なるクラスタ間で感染するための通過率。万分率で指定。10000で通常感染、0で完全遮蔽。
}ACTIONS;
typedef struct STATISTICS{//統計データを格納する
int day;
int n;//グループの人数
int living;//健康な生者の人数
int patient;//患者の人数(厳密には感染者、言い換えるとウイルス保有者の人数)
int infect_new;//当日の新規感染者の人数
int onset_new;//当日の新規発症者の人数
int dead;// 死者数
int sars_v[ VIRUS_N ];//各ウイルスの感染者数。0なら発生していない、1以上なら発生している
int infect_total;//これまで感染した人の総計
int sars_v_total[ VIRUS_N ];//各ウイルスのこれまでの感染者数総計。
}STATISTICS;
#define EMERGENCY_FLG_ON 1
#define EMERGENCY_FLG_OFF 0
#define EMERGENCY_RATIO (0.7) //緊急事態宣言によって自粛した場合のRtの減り方。この値をかける
// 「医療崩壊のリスク定義」
// 医療崩壊のリスクもシミュレーション結果であるが、その算出は難しい。
// 医療崩壊する条件は、「人員確保」「病床確保」「機材確保」などがリスク低減要因として、
// 「市民の外出」「患者の暴言・わがまま」「医療スタッフの辞職」「重症者の増加」「患者数の絶対数の増加」などがリスク増大要因
// として働く。これらの条件は都道府県市区町村によって違うため、精緻に計算することはそもそもナンセンスである。
// 本シミュレーションでは、医療崩壊のリミットとして、「患者数の総数」を基準とする。
// 患者数の総数はどの条件よりもはっきりしており、かつどの要因よりも危機的である(急激に増大しうる)からである。
// 全人口に対する患者数がどれくらいになると医療崩壊を起こしうるかは大阪の例を参考にする。
// 大阪府民の数は882.3万人、うち重症病床248床、軽症中等症病床1781床、宿泊療養施設3059床を確保している。
// 2021年04月18日の時点で、これらはそれぞれ244床(98.4%)、1387床(77.9%)、1013床(33.1%)、が使用されている。
// 重症者のリミットは逼迫しており、人工心肺装置ECMOの数に限界が来たときにトリアージが始まり、宿泊療養施設の一部か
// 一般ホテルなどの大型宿泊施設が重症者用に転用されるのではないかと思われる
// 最初のトリアージはECMOが必要なほどの患者は切り捨てることになると思われる。
// このラインは大阪特有であり、ECMOの総数がより多い県では大阪より危険度は低く、総数が少なければもっと早く医療崩壊する。
// また、重症化は本人の体質以外に一度の感染で侵入したウイルスの絶対数も影響するため、
// 無防備な県民性があれば、それが重症者の増加につながり、医療崩壊が高リスクである。
// そこで、次のマクロでは、大阪府民数に対する「感染者数」で医療崩壊ラインを推定することにする。
// 当然、他の県では条件が異なるのであるが、重要なのは「患者数が多ければいずれ崩壊する」ということであり、
// 感染者数が爆発的に増大すれば、たとえ重症者が少なくとも医療スタッフが潰れてしまうためこの計算がもっとも
// 折衷案として妥当であろう。
// 医療崩壊ラインを決める割合。算出の背景については前記「医療崩壊のリスク定義」参照。
// 感染者数がこの割合（万分率）を超えたときに医療崩壊と判定する
// 2021/04/18 4111
// 2021/04/17 4802
// 2021/04/16 4511
// 2021/04/15 4515
// 2021/04/14 5482
// 2021/04/13 3449
// 2021/04/12 2113
// 七日間の感染者数合計 28983
// 大阪府民 8823000
// 人口比 0.00328
// 万分率 32.85
// 市民1万人あたり33人の感染者が出たら医療崩壊
#define MEDICALCAPACITYOVER 33
#define MEDICALCAPACITYOVER_FLG_ON 1
#define MEDICALCAPACITYOVER_FLG_OFF 0
// シミュレーション対象集団の構造体FLOCKは2つ確保する。
// これは、n日目の状態はn-1日目から生成するが、n-1日目のデータに直接書き込むと
// 二重評価や異常な書き換えを防ぐことが困難だからである。
// 片方yesterdayは読み取り専用、もう片方todayは書き込み専用である
// 翌日の処理は、todayの内容を丸ごとyesterdayへコピーしてから、感染等の実行を演算し、todayへ書き込んでいく
typedef struct DAILYCALENDAR{//2日間の状態(昨日と今日)を保持し、昨日の状態から今日の状態を生成する
FLOCK *today;
FLOCK *yesterday;
int day;//経過日数。シミュレーション開始時点で0とし、1日経過ごとに加算。
int emergency_flg;//緊急事態宣言の有無
int medicalover_flg;
ACTIONS actions;
STATISTICS statistics;//todayの統計情報。参照前に必ずDAILYCALENDAR.dayとSTATISTICS.dayが一致していることを確認
}DAILYCALENDAR;
HUMANS * humans_make(int n);
int humans_free(HUMANS *humans);
FLOCK * flock_make(int n);
int flock_free(FLOCK *flock);
DAILYCALENDAR * DailyCalendar_make( int n, ACTIONS *actions );
int DailyCalendar_free(DAILYCALENDAR *dailycalendar );
int dayAdvance(DAILYCALENDAR *dailycalendar);
int actions_property_set( ACTIONS *actions, int mask, int cons_strong, int cons_hide, int cons_middle, int cons_severe, int cons_death, int self_restraint, int self_restraint_limit, int personality_Influencer, int personality_alone, int personality_middle, int personality_Influencer_weight, int personality_alone_weight, int immunity, int modulo_M, int mod_risk_wait );
int flock_property_set( FLOCK *flock, ACTIONS *actions );
int flock_property_print( FLOCK *flock, FILE *fp);
int actions_property_print(DAILYCALENDAR *dailycalendar, FILE *fp);
#ifndef MERSENNE_TWISTER
int irand(int n);
int roll_the_dice(int denominator, int numerator);
int roll_the_dice2(int denominator, int numerator);
#endif
int find_infection(FLOCK *flock, int *i);
int find_Virus_type(HUMANS *humans, int type);
int is_Virus_Rt(SARS_V *sars_v);
int is_Virus_Population(SARS_V *sars_v);
int daily_condition(DAILYCALENDAR *dailycalendar);
int daily_condition_sub_human_recovery(HUMANS *humans);
FILE * now_date_fopen(TCHAR *fname, TCHAR *ext);
FILE * now_date_fopen_n(TCHAR *fname, TCHAR *ext);
int get_yymmddhhmmss_int_now( int *year, int *month, int *day, int *hh, int *mm, int *ss );
FILE * now_date_micro_fopen(TCHAR *fname, TCHAR *ext);
int get_yymmddhhmmssms_int_now( int *year, int *month, int *day, int *hh, int *mm, int *ss, int *ms );
int File_Existence(TCHAR *fname);
int fprintf_dailycalendar(FILE *fp, DAILYCALENDAR *dailycalendar);
int fprintf_dailycalendar_flock(FILE *fp, FLOCK *flock);
int fprintf_dailycalendar_humans(FILE *fp, HUMANS *humans);
int fprintf_daily_statistics(FILE *fp, DAILYCALENDAR *dailycalendar);
int fprintf_humans_statistics(FILE *fp, DAILYCALENDAR *dailycalendar);
int get_statistics_flock(DAILYCALENDAR *dailycalendar, STATISTICS *statistics);
int fprintf_daily_virus_Frequency_distribution(FILE *fp, DAILYCALENDAR *dailycalendar);
int fprintf_daily_statistics2(FILE *fp, DAILYCALENDAR *dailycalendar);
int fprintf_daily_statistics3(FILE *fp, DAILYCALENDAR *dailycalendar);
int bubbleSort(int *index, int *sars_v, int len);
int fprintf_daily_virus_Frequency_distribution_Sorting(FILE *fp, DAILYCALENDAR *dailycalendar);
int fprintf_daily_graph(FILE *fp, DAILYCALENDAR *dailycalendar, int world_line);
int fprintf_daily_graph2(FILE *fp, DAILYCALENDAR *dailycalendar, int world_line);
int fprintf_daily_Calendar_setting(FILE *fp, int world_line, DAILYCALENDAR *dailycalendar);
int fprintf_daily_Calendar_setting_act(FILE *fp, int world_line, ACTIONS *actions, int n);
int fprintf_daily_Calendar_setting2(FILE *fp, int world_line);
int myPCshutdown(void);
//SARS_V * SARS_v_make(int type);
int virus_infection(HUMANS *humans, SARS_V *sars_v);
int daily_infection(DAILYCALENDAR *dailycalendar);
int daily_infection_sub_viruscopy(SARS_V *sars_v1, SARS_V *sars_v2, int mask);
int daily_infection2(DAILYCALENDAR *dailycalendar);
int daily_infection2_sub_viruscopy(SARS_V *sars_v1, SARS_V *sars_v2, int mask);
int daily_infection2_sub_Rt_Decrease_n(HUMANS *humans);
int daily_infection2_sub_Rt_Decrease_sars_v_n(SARS_V *sars_v);
int daily_infection_sub_virus_invasion(HUMANS *humans, SARS_V *sars_v);
int daily_infection_sub_Rt_0(HUMANS *humans);
int daily_infection_sub_Rt_n(HUMANS *humans, double n);
int daily_mutation(DAILYCALENDAR *dailycalendar);
int daily_mutation2(DAILYCALENDAR *dailycalendar);
int daily_proliferation(DAILYCALENDAR *dailycalendar);
int daily_event( DAILYCALENDAR *dailycalendar );
int daily_event_sub_EmergencyOnCheck(DAILYCALENDAR *dailycalendar, int patient);
int daily_event_sub_EmergencyOffCheck(DAILYCALENDAR *dailycalendar, int patient);
int daily_event_sub_PatientCount(DAILYCALENDAR *dailycalendar, int *patient);
int daily_event_sub_MedicalCapacityOverCheck(DAILYCALENDAR *dailycalendar, int patient);
int virus_injection(HUMANS *humans, int virus_type, int population);
int virus_injection2(HUMANS *humans, int virus_type, int population);
int virus_invasion(HUMANS *humans, int virus_type, int population);
int virus_invasion_sub_cons_hide(HUMANS *humans);
int daily_statistics_calc(DAILYCALENDAR *dailycalendar);
int infection_dice(HUMANS *humans );
int Simulation_Run( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line );
int Simulation_Run2( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line );
int Simulation_Run3( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line );
int Simulation_Run4( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line );
int actions_property_set_ex1( ACTIONS *actions );
int actions_property_set_ex2( ACTIONS *actions );
int actions_property_set_ex3( ACTIONS *actions );
int actions_property_set_ex4( ACTIONS *actions );
int Simulation_Run1_ex1(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run2_ex1(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run3_ex1(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run3_ex2(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run4_ex1(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run5_ex1(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run5_ex2(FILE *fp_log, FILE *fp_glp, int seed );
int Simulation_Run5_ex3(FILE *fp_log, FILE *fp_glp, int seed );
int main(int argc, char **argv){
//int len, res, ret;
//TCHAR str[100];
//time_t t;
FILE *fp_log, *fp_glp;
int seed;
//int n, nday,world_line, ret;//何度目の世界線をシミュレートしているのか
//ACTIONS actions;
setlocale( LC_ALL, "Japanese" );//ロケール（地域言語）を日本語でセット
if(argc&&argv){;}
_tprintf( _T("**********************************************************\n") );
_tprintf( _T("シミュレーション開始stt\n") );
fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
fp_glp = now_date_fopen( _T("複数回シミュレーショングラフ用データ"), _T(".csv") );
seed = time(NULL);
//seed = 1;
//Simulation_Run1_ex1( fp_log, fp_glp, seed );
//Simulation_Run2_ex1( fp_log, fp_glp, seed );
//Simulation_Run3_ex1( fp_log, fp_glp, seed );
//Simulation_Run3_ex2( fp_log, fp_glp, seed );
///Simulation_Run4_ex1( fp_log, fp_glp, seed );
//Simulation_Run5_ex1( fp_log, fp_glp, seed );//医療崩壊のチェック
//Simulation_Run5_ex2( fp_log, fp_glp, seed );//医療崩壊のチェック
Simulation_Run5_ex3( fp_log, fp_glp, seed );//医療崩壊のチェック。クラスター間感染率変化
_tprintf( _T("シミュレーション終了end\n") );
fclose( fp_log );
fclose( fp_glp );
//_fgetts(str, 10, stdin);
#ifdef ShutdownEnd
// これを有効にするとプログラム終了後にPCをシャットダウンする
myPCshutdown();
#endif
return 1;
}
/***************************************************************
*
* ========目的========
* HUMANS型のデータを生成し、初期化してポインタを返す
* 人口を指定して人間集団の構造体配列を確保、ポインタを返す
* エラーではNULLを返す
*
* ========引数========
* 人数n
*
* ========詳細========
*
*/
HUMANS * humans_make(int n){
HUMANS * humans=NULL;
int i, x;
if( n<1 ) return NULL;
humans = (HUMANS *)malloc( sizeof(HUMANS)*(n) );
if( humans==NULL ) return NULL;
for(i=0;i<n;i++){
humans[i].no = i;
humans[i].cons = -1;// 体質 Constitution 感染後の予後を決める
humans[i].count = -1;//感染後の経過日数。感染した段階で0になる。通常は-1
for(x=0;x<VIRUS_N;x++){
humans[i].sars_v[ x ].type = VIRUS_0;
humans[i].sars_v[ x ].Population = 0;
humans[i].sars_v[ x ].Rt = 0.0;
}
humans[i].virusn = 0;//何種類に感染しているか
humans[i].condition = CONDITION_ALIVE;// 生きて動いてる
humans[i].mask = MASK_OFF;// マスクをつけない
humans[i].self_restraint = SELF_RESTRAINT_ON;//自粛する
humans[i].personality = PERSONALITY_MIDDLE;
humans[i].immunity = IMMUNITY_OFF;// 免疫を持たない
}
return humans;
}
/***************************************************************
*
* ========目的========
* HUMANS型のデータの使用メモリを開放する
*
*/
int humans_free(HUMANS *humans){
if( humans==NULL ) return 0;
free( humans );
return 1;
}
/***************************************************************
*
* ========目的========
* FLOCK 型のデータを生成し、初期化してポインタを返す
*
*/
FLOCK * flock_make(int n){
FLOCK * flock=NULL;
if( n<1 ) return NULL;
flock = (FLOCK *)malloc( sizeof(FLOCK) );
if( flock==NULL ) return NULL;
flock->humans = humans_make( n );
if( flock->humans==NULL ) return NULL;
flock->n = n;
return flock;
}
/***************************************************************
*
* ========目的========
* FLOCK 型のデータの使用メモリを開放する
*
*/
int flock_free(FLOCK *flock){
if( flock==NULL ) return 0;
humans_free( flock->humans );
free( flock );
return 1;
}
/***************************************************************
*
* ========目的========
* DAILYCALENDAR 型のデータを生成し、初期化してポインタを返す
*
*/
DAILYCALENDAR * DailyCalendar_make( int n, ACTIONS *actions ){
DAILYCALENDAR *dailycalendar = NULL;
if( n<3 ){//3以下のクラスターなど計算に値しないのでデフォルト値を適用
n = MAN_N;
}
if( actions==NULL ) return NULL;
dailycalendar = (DAILYCALENDAR *)malloc( sizeof(DAILYCALENDAR) );
dailycalendar->today = flock_make( n );
dailycalendar->yesterday = flock_make( n );
dailycalendar->day = 0;
dailycalendar->emergency_flg = EMERGENCY_FLG_OFF;
dailycalendar->medicalover_flg = MEDICALCAPACITYOVER_FLG_OFF;
dailycalendar->actions = *actions;
dailycalendar->statistics.day = -1;
flock_property_set( dailycalendar->today, actions );
if( dailycalendar->today==NULL || dailycalendar->yesterday==NULL ){// mallocが失敗しているなら実行は不可能
return NULL;
}
return dailycalendar;
}
/***************************************************************
*
* ========目的========
* DAILYCALENDARの使用メモリを開放する
*
*/
int DailyCalendar_free(DAILYCALENDAR *dailycalendar ){
if( dailycalendar==NULL ){
return 0;
}
flock_free( dailycalendar->today );
flock_free( dailycalendar->yesterday );
free( dailycalendar );
return 1;
}
/***************************************************************
*
* ========目的========
* 一日進めるために、todayの内容をyesterdayにコピーする
* 内容がまったく同じになれば今後も変わることはないので0を返す
* エラーは-1を返す
*/
int dayAdvance(DAILYCALENDAR *dailycalendar){
size_t size;
if( dailycalendar==NULL ) return -1;
size = sizeof( dailycalendar->yesterday->humans[0] ) * (dailycalendar->yesterday->n);
if( memcmp(dailycalendar->yesterday->humans, dailycalendar->today->humans, size )==0 ){
//_tprintf( _T("全く同じ\n") );
return 0;
}
//_tprintf( _T("size=%d\n"), size );
memcpy(dailycalendar->yesterday->humans, dailycalendar->today->humans, size);
(dailycalendar->day)++;
return 1;
}
/***************************************************************
*
* ========目的========
* 行動や体質の設定条件を指定する構造体を適切に設定する
*
*
* ========引数========
*
*
* ========詳細========
*
*
*
*/
int actions_property_set( ACTIONS *actions,
int mask_m,
int cons_strong_m,
int cons_hide_m,
int cons_middle_m,
int cons_severe_m,
int cons_death_m,
int self_restraint_m,
int self_restraint_limit_m,
int personality_Influencer_m,
int personality_alone_m,
int personality_middle_m,
int personality_Influencer_weight_m,
int personality_alone_weight_m,
int immunity_m,
int modulo_M,
int mod_risk_wait_m ){
if( actions==NULL ) return 0;
if( mask_m>=0 && mask_m<=PER_DEF ){
actions->mask_m = mask_m;
}
else{
_tprintf( _T("maskの値が範囲外 0 <= %d <= %d\n"), mask_m, PER_DEF );
return 0;
}
if( (cons_strong_m + cons_hide_m + cons_middle_m + cons_severe_m + cons_death_m)==PER_DEF &&
(cons_strong_m>=0 && cons_hide_m>=0 && cons_middle_m>=0 && cons_severe_m>=0 && cons_death_m>=0 ) ){//100分率がきちんと指定されている
actions->cons_strong_m = cons_strong_m;
actions->cons_hide_m = cons_hide_m;
actions->cons_middle_m = cons_middle_m;
actions->cons_severe_m = cons_severe_m;
actions->cons_death_m = cons_death_m;
}
else{
_tprintf( _T("consの値が範囲外 合計 %d になるべきです %d \n"), PER_DEF, cons_strong_m + cons_hide_m + cons_middle_m + cons_severe_m + cons_death_m );
return 0;
}
if( self_restraint_m>=0 && self_restraint_m<=PER_DEF ){//100分率がきちんと指定されている
actions->self_restraint_m = self_restraint_m;
}
else{
_tprintf( _T("self_restraint_mの値が範囲外 0 <= %d <= %d\n"), self_restraint_m, PER_DEF );
return 0;
}
if( self_restraint_limit_m>=0 && self_restraint_limit_m<=PER_DEF ){//100分率がきちんと指定されている
actions->self_restraint_limit_m = self_restraint_limit_m;
}
else{
_tprintf( _T("self_restraint_limit_mの値が範囲外 0 <= %d <= %d\n"), self_restraint_limit_m, PER_DEF );
return 0;
}
if( (personality_Influencer_m + personality_alone_m + personality_middle_m)==PER_DEF &&
(personality_Influencer_m>=0 && personality_alone_m>=0 && personality_middle_m>=0) ){//100分率がきちんと指定されている
actions->personality_Influencer_m = personality_Influencer_m;
actions->personality_alone_m = personality_alone_m;
actions->personality_middle_m = personality_middle_m;
}
else{
_tprintf( _T("personality_の値が範囲外 合計 %d になるべきです %d \n"), PER_DEF, personality_Influencer_m + personality_alone_m + personality_middle_m );
return 0;
}
if( personality_Influencer_weight_m>=0 && personality_alone_weight_m>=0 ){
actions->personality_Influencer_weight_m = personality_Influencer_weight_m;
actions->personality_alone_weight_m = personality_alone_weight_m;
}
else{
_tprintf( _T("personality_.._weightの値が範囲外 合計 より大きい値0 になるべきです %d %d \n"), personality_Influencer_weight_m, personality_alone_weight_m );
return 0;
}
if( immunity_m>=0 && immunity_m<=PER_DEF ){
actions->immunity_m = immunity_m;
}
else{
_tprintf( _T("immunity_mの値が範囲外 0 <= %d <= %d\n"), immunity_m, PER_DEF );
return 0;
}
if( modulo_M>0 ){
actions->modulo_M = modulo_M;
}
else{
_tprintf( _T("modulo_Mの値が範囲外 0 < %d \n"), modulo_M );
return 0;
}
if( mod_risk_wait_m>=0 && mod_risk_wait_m<=PER_DEF ){//100分率がきちんと指定されている
actions->mod_risk_wait_m = mod_risk_wait_m;
}
else{
_tprintf( _T("mod_risk_wait_mの値が範囲外 0 <= %d <= %d\n"), mod_risk_wait_m, PER_DEF );
return 0;
}
return 1;
}
/***************************************************************
*
* ========目的========
*
* 人間群に対して、性格や行動特徴や遺伝的特徴を設定する
*
* ========引数========
*
*
* ========詳細========
*
*
*
*/
int flock_property_set( FLOCK *flock, ACTIONS *actions ){
int i;
int num;
if( flock==NULL || actions==NULL ) return 0;
// マスクの有無
num = flock->n * (actions->mask_m) / PER_DEF;
//_tprintf( _T("num= %d \n"), num );
for(i=0;i<num;i++){
flock->humans[i].mask = MASK_ON;
}
num = i + flock->n - num;
for(;i<num;i++){
flock->humans[i].mask = MASK_OFF;
}
// 体質
num = flock->n * (actions->cons_strong_m) / PER_DEF;
for(i=0;i<num;i++){
flock->humans[i].cons = CONS_STRONG;
}
num = i + flock->n * (actions->cons_hide_m) / PER_DEF;
for(;i<num;i++){
flock->humans[i].cons = CONS_HIDE;
}
num = i + flock->n * (actions->cons_middle_m) / PER_DEF;
for(;i<num;i++){
flock->humans[i].cons = CONS_MIDDLE;
}
num = i + flock->n * (actions->cons_severe_m) / PER_DEF;
for(;i<num;i++){
flock->humans[i].cons = CONS_SEVERE;
}
num = i + flock->n * (actions->cons_death_m) / PER_DEF;
//_tprintf( _T("num= %d \n"), num );
for(;i<num;i++){
flock->humans[i].cons = CONS_DEATH;
}
// 行動自粛
num = flock->n * (actions->self_restraint_m) / PER_DEF;
for(i=0;i<num;i++){
flock->humans[i].self_restraint = SELF_RESTRAINT_ON;
//_tprintf( _T("行動自粛= %d \n"), i );
}
num = i + flock->n - num;
for(;i<num;i++){
flock->humans[i].self_restraint = SELF_RESTRAINT_OFF;
}
// 行動特性(性格)
num = flock->n * (actions->personality_Influencer_m) / PER_DEF;
for(i=0;i<num;i++){
flock->humans[i].personality = PERSONALITY_INFLUENCER;
}
num = i + flock->n * (actions->personality_alone_m) / PER_DEF;
for(;i<num;i++){
flock->humans[i].personality = PERSONALITY_ALONE;
}
num = i + flock->n * (actions->personality_middle_m) / PER_DEF;
for(;i<num;i++){
flock->humans[i].personality = PERSONALITY_MIDDLE;
}
// 免疫の有無
num = flock->n * (actions->immunity_m) / PER_DEF;
for(i=0;i<num;i++){
flock->humans[i].immunity = IMMUNITY_ON;
}
num = i + flock->n - num;
for(;i<num;i++){
flock->humans[i].immunity = IMMUNITY_OFF;
}
flock_property_print(flock, stdout);
return 1;
}
/***************************************************************
*
* ========目的========
*
* 人間群に対して、設定した情報をカウントし、どの属性が何人に適用されているか表示する
* 本関数は元々デバッグ用の確認関数であったが、後で統計的な検討をする際にこの情報が必要であることに気付いたため
* 急遽flock_property_print( FLOCK *flock) から flock_property_print( FLOCK *flock, FILE *fp)
* に拡張し、標準出力ファイル出力兼用にした。同じ内容を毎回計算するのはスマートじゃないけど、
* カウント用関数を分離するほど重要でもないしボトルネックにもならないし、余力があったら手を入れようと思う。
* ========引数========
*
*
* ========詳細========
*
*/
int flock_property_print( FLOCK *flock, FILE *fp){
ACTIONS actions;
int i;
int mask, nomask;
int immunity, noimmunity;
if( flock==NULL ) return 0;
mask = 0;
nomask = 0;
for(i=0;i<flock->n;i++){
if( flock->humans[i].mask == MASK_ON ) mask++;
if( flock->humans[i].mask == MASK_OFF ) nomask++;
}
_ftprintf( fp, _T("マスクあり %d 人 マスク無し %d 人\r\n"), mask, nomask );
// 体質
actions.cons_strong_m = 0;
actions.cons_hide_m = 0;
actions.cons_middle_m = 0;
actions.cons_severe_m = 0;
actions.cons_death_m = 0;
for(i=0;i<flock->n;i++){
if( flock->humans[i].cons == CONS_STRONG ) (actions.cons_strong_m)++;
if( flock->humans[i].cons == CONS_HIDE ) (actions.cons_hide_m )++;
if( flock->humans[i].cons == CONS_MIDDLE ) (actions.cons_middle_m)++;
if( flock->humans[i].cons == CONS_SEVERE ) (actions.cons_severe_m)++;
if( flock->humans[i].cons == CONS_DEATH ) (actions.cons_death_m )++;
}
_ftprintf( fp, _T("無敵体質 %4d 人\r\n"), actions.cons_strong_m );
_ftprintf( fp, _T("忍者体質 %4d 人\r\n"), actions.cons_hide_m );
_ftprintf( fp, _T("普通体質 %4d 人\r\n"), actions.cons_middle_m );
_ftprintf( fp, _T("重症体質 %4d 人\r\n"), actions.cons_severe_m );
_ftprintf( fp, _T("死ぬ体質 %4d 人\r\n"), actions.cons_death_m );
// 行動自粛
actions.self_restraint_m = 0;
for(i=0;i<flock->n;i++){
if( flock->humans[i].self_restraint == SELF_RESTRAINT_ON ) (actions.self_restraint_m)++;
}
_ftprintf( fp, _T("自粛する人 %d 人\r\n"), actions.self_restraint_m );
// 行動特性(性格)
actions.personality_Influencer_m = 0;
actions.personality_alone_m = 0;
actions.personality_middle_m = 0;
for(i=0;i<flock->n;i++){
if( flock->humans[i].personality == PERSONALITY_INFLUENCER ) (actions.personality_Influencer_m)++;
if( flock->humans[i].personality == PERSONALITY_ALONE ) (actions.personality_alone_m )++;
if( flock->humans[i].personality == PERSONALITY_MIDDLE ) (actions.personality_middle_m )++;
}
_ftprintf( fp, _T("活発な性格 %4d 人\r\n"), actions.personality_Influencer_m );
_ftprintf( fp, _T("隠者な性格 %4d 人\r\n"), actions.personality_alone_m );
_ftprintf( fp, _T("普通な性格 %4d 人\r\n"), actions.personality_middle_m );
// 免疫の有無
immunity = 0;
noimmunity = 0;
for(i=0;i<flock->n;i++){
if( flock->humans[i].immunity == IMMUNITY_ON ) immunity++;
if( flock->humans[i].immunity == IMMUNITY_OFF ) noimmunity++;
}
_ftprintf( fp, _T("免疫あり %d 人 免疫無し %d 人\r\n"), immunity, noimmunity );
return 1;
}
/***************************************************************
*
* ========目的========
*
* 設定された行動情報を出力する
* ========引数========
*
*
* ========詳細========
*
*/
int actions_property_print(DAILYCALENDAR *dailycalendar, FILE *fp){
if( dailycalendar==NULL ) return 0;
_ftprintf( fp, _T("全人口は %d 個のクラスタに分けられる\r\n"), dailycalendar->actions.modulo_M );
_ftprintf( fp, _T("1クラスタ約 %d人\r\n"), (dailycalendar->today->n)/(dailycalendar->actions.modulo_M) );
_ftprintf( fp, _T("クラスタ間感染成功率は %d/%d =%f\r\n"), (dailycalendar->actions.mod_risk_wait_m), PER_DEF, (double)(dailycalendar->actions.mod_risk_wait_m)/PER_DEF );
return 1;
}
#ifndef MERSENNE_TWISTER
/***************************************************************
*
* ========目的========
* 0〜n-1までのn種類の整数の乱数を返す
*
*
* ========引数========
*
* ========詳細========
* ( (1.0 / (RAND_MAX + 1.0) ) * rand() ) で乱数を一旦0〜1.0 の範囲に限定し
* これにnを掛ける事でバラつきをよりランダムにする。
* 単純に剰余算を用いると、乱数のばらつきが小さくなる。
* 大数の法則が現れるような試行回数であれば、十分に「均等に」ランダム出現となる
* xorshiftやメルセンヌツイスタの方が乱数生成としては優れているが、
* 本シミュレーションでは乱数は毎回のダイス判定に一回ずつ使うのみであり
* 恐らくrandで十分である。むしろ、なれない方法を使ってうっかり乱数の不均等を作りこむ方が
* シミュレーションに悪影響を与える。
*/
int irand(int n){
int re;
re = ( (1.0 / (RAND_MAX + 1.0) ) * rand() ) * n;
return re;
}
/***************************************************************
*
* ========目的========
* ダイスの目の数とあたり判定の個数を指定し、当たったかどうかだけ返す
* 別の表現をすれば、分子/分母の確率で1を返し、(分母-分子)/分母の確率で0を返す
*
* ========引数========
*
* ========詳細========
* 例えば、分母denominator に6、分子numerator に2を指定すれば、
* 1/3の確率で1を返し、2/3の確率で0を返す。
* 乱数シミュレーションでこの値を間違えると重大かつ気付きにくいバグに繋がるので、
* 引数の矛盾についてはきっちり検証する
*
*/
int roll_the_dice(int denominator, int numerator){
int re;
if( (denominator < numerator) || (denominator < 1) || (numerator < 0) ){
_tprintf( _T("roll_the_dice(%d, %d)において引数が矛盾しています\n"), denominator, numerator );
}
re = ( (1.0 / (RAND_MAX + 1.0) ) * rand() ) * denominator;
if( re < numerator ) return 1;
return 0;
}
/***************************************************************
*
* ========目的========
* ダイスの目の数とあたり判定の個数を指定し、当たったかどうかだけ返す
* 別の表現をすれば、分子/分母の確率で1を返し、(分母-分子)/分母の確率で0を返す
*
* ========引数========
*
* ========詳細========
* 例えば、分母denominator に6、分子numerator に2を指定すれば、
* 1/3の確率で1を返し、2/3の確率で0を返す。
* 乱数シミュレーションでこの値を間違えると重大かつ気付きにくいバグに繋がるので、
* 引数の矛盾についてはきっちり検証する
*
*/
int roll_the_dice2(int denominator, int numerator){
int re;
if( (denominator < numerator) || (denominator < 1) || (numerator < 0) ){
_tprintf( _T("roll_the_dice(%d, %d)において引数が矛盾しています\n"), denominator, numerator );
}
re = ( (1.0 / (RAND_MAX + 1.0) ) * rand() ) * denominator;
if( re < numerator ) _tprintf( _T("roll_the_dice(%d分の%d) であたり\n"), denominator, numerator );
else _tprintf( _T("roll_the_dice(%d分の%d) ではずれ\n"), denominator, numerator );
if( re < numerator ) return 1;
return 0;
}
#endif
/***************************************************************
*
* ========目的========
* 感染者を見つけて、その個人IDを *i に設定して1を返す。
* 感染者がいなければ0を返す。
* iは0初期化して呼び出すこと
* iのインクリメントは「呼び出し側で行う」理由は以下詳細にて。
*
* ========引数========
*
* ========詳細========
* 本関数は、検索開始位置のインクリメントを呼び出し側に任せている。
* インデックスであるi変数を-1初期化として本関数で毎回インクリメントする仕様にすれば
* 呼び出し元の記述が1行減るのであるが、以下のヒューマンバグを予防するため意図的に面倒にしている。
* ・-1初期化を義務とした場合、うっかり0初期化しても先頭要素を無視するだけでバグが顕在化しにくい
* ・0初期化の方が感覚的にわかりやすい、やりやすい
* ・呼び出し元のインクリメント忘れは、テスト時に高確率でハングするので発見しやすい
* ミスは発見しやすい方がよいのだ。
*
*/
int find_infection(FLOCK *flock, int *i){
if( flock==NULL || i==NULL ) return 0;
if( (*i) < 0 || (*i)-1 > flock->n ){
_tprintf( _T("不正なインデックス。 %d 記述を見直せ\n"), (*i) );
}
for(;(*i)<flock->n;(*i)++){
//_tprintf( _T("i=%d\n"), *i );
if( flock->humans[(*i)].virusn>0 ){//感染発見
//_tprintf( _T("感染者発見\n") );
return 1;
}
}
return 0;
}
/***************************************************************
*
* ========目的========
* 特定の人間の中に、指定された種類のウイルス種が発生・感染しているか調べる
* ウイルスがいなければ0を返す。
* 定義されていないウイルス種を指定された場合はエラーで-1を返す
*
* ========引数========
*
* ========詳細========
*
*/
int find_Virus_type(HUMANS *humans, int type){
int x;
if( humans==NULL ) return -1;
if( type!=VIRUS_0 && type!=VIRUS_1 && type!=VIRUS_2 && type!=VIRUS_3 && type!=VIRUS_4 ){
_tprintf( _T("%d そんなものはない…なんだよ、%dって…\n"), type, type );
return -1;
}
for(x=0;x<VIRUS_N;x++){
if( humans->sars_v[x].type == type ) return 1;
}
return 0;
}
/***************************************************************
*
* ========目的========
* SARS_V sarsv 構造体のウイルスデータの中に実効再生産数Rtが残っているかどうか判定する
* 一個でもRtが0より大きい値で残っていれば1を返す。
*
* ========引数========
*
* ========詳細========
* 感染試行を行ってもよいかどうかの判定なので、一個でも実効再生産数Rtが残っているウイルスがいれば
* 1を返す。
*
*/
int is_Virus_Rt(SARS_V *sars_v){
int x;
if( sars_v==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
if( sars_v[x].Rt > 0 ) return 1;
}
return 0;
}
/***************************************************************
*
* ========目的========
* SARS_V sarsv 構造体のウイルスデータの中で、感染可能なだけウイルス個体数が増えているかどうか判定する
* 一種類でも感染可能な数を超えていれば1を返す。
*
* ========引数========
*
* ========詳細========
*
*
*/
int is_Virus_Population(SARS_V *sars_v){
int x;
if( sars_v==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
if( sars_v[x].Population >= INVASION_LIMIT ){
return 1;
}
}
return 0;
}
/***************************************************************
*
* ========目的========
* 一日分の病状の変化を行う
* 本人の体質と経過日数に応じて、「現状維持」「重症化」「死亡」の変化を起こす
* また、緊急事態宣言に応じた行動変容で外出自粛する場合は
*
* ========引数========
*
* ========詳細========
*
*/
int daily_condition(DAILYCALENDAR *dailycalendar){
int i;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
i = 0;
flock = dailycalendar->today;
while( find_infection( flock, &i ) ){//感染者を探す
if( flock->humans[i].count == TIME_LIMIT_INCUBATION ){// 発症タイミング
//_tprintf( _T("発症\n") );
flock->humans[i].condition = CONDITION_INCUB;
}
else if( flock->humans[i].count == TIME_LIMIT_SEVERE && //重症化リミット通過
flock->humans[i].condition == CONDITION_INCUB && //発症中
(flock->humans[i].cons==CONS_SEVERE || flock->humans[i].cons==CONS_DEATH) ){//重症化体質であるまたは死亡体質である
//_tprintf( _T("重症化\n") );
flock->humans[i].condition = CONDITION_SEVEE;
}
else if( flock->humans[i].count == TIME_LIMIT_DEATH && //死亡リミット通過
flock->humans[i].condition == CONDITION_SEVEE && //重症中
flock->humans[i].cons==CONS_DEATH ){// 死亡体質である
//_tprintf( _T("死亡\n") );
flock->humans[i].condition = CONDITION_DEATH;
}
else if( flock->humans[i].count == TIME_LIMIT_RECOVERY && flock->humans[i].condition == CONDITION_INCUB){//治癒
//_tprintf( _T("治癒%d\n"), i );
flock->humans[i].condition = CONDITION_ALIVE;
flock->humans[i].immunity = IMMUNITY_ON;
daily_condition_sub_human_recovery( &(flock->humans[i]) );
}
i++;
}
return 0;
}
/***************************************************************
*
* ========目的========
* daily_condition()の下請けとして、人間の感染状態を解き、治癒状態にする
*
*
* ========引数========
*
* ========詳細========
*
*/
int daily_condition_sub_human_recovery(HUMANS *humans){
int x;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
humans->sars_v[x].Population = 0;
//humans->sars_v[x].type = VIRUS_0;
}
humans->virusn = 0;
humans->condition = CONDITION_ALIVE;
humans->immunity = IMMUNITY_ON;
return 1;
}
/***************************************************************
*
* ========目的========
* :::::::::注意::::::::::
* 本関数は軽量なれど、時刻を秒単位でしか認識できない。マイクロ秒まで必要であれば
* FILE * now_date_micro_fopen(TCHAR *fname, TCHAR *ext)を使う
*
* 呼び出された瞬間の日時を付加した一意のファイル名を生成し、
* そのファイル名をもってファイルを書き込みでopenしファイルポインタを返す
* fnameには基本ファイル名を、extに拡張子指定する。拡張子にはドットを含むこと。
* fname、extは書き込み不可のリテラル文字列でかまわないし
* NULLでもよい。その場合は日時だけでファイル名が構成される。
*
* ========引数========
* TCHAR *fname 基本ファイル名。拡張子は含まない。NULLでもよいがあまり意味がない。
* TCHAR *ext 拡張子。NULLでもよい。
*
* ========詳細========
* 例えば、now_date_fopen("ログファイル", ".log"); として2021/11/12 13:31:22に呼び出せば、
* ログファイル(2021年11月12日13時31分22秒).log
* としてファイルが開かれる。
* now_date_fopen("", ""); として2021/11/12 13:31:22に呼び出せば、
* (2021年11月12日13時31分22秒)
* というファイルが開かれる。
* もしも、同名のファイルがあった場合にはエラーとしてNULLを返す
* 別の実装として、1秒待ってから再起呼び出しというのもありえるが、この関数が衝突するような実行環境下では
* そもそもかなり呼び出しが頻繁であるということであり、大量の待ち行列が発生することを考慮しなければならない。
* 大量の待ち行列では、スタックオーバーフローの懸念もある。特に、本関数はTCHAR fname_bufという
* 大きめの自動変数を使っているので、スタックの負荷はでかい。
* 秒以下刻み環境下では、他の関数を実装すべし
*
*/
FILE * now_date_fopen(TCHAR *fname, TCHAR *ext){
int year, month, day, hh, mm, ss;
TCHAR fname_buf[128*256];
FILE *fp=NULL;
get_yymmddhhmmss_int_now( &year, &month, &day, &hh, &mm, &ss );
_sntprintf(fname_buf, 128*256-2, _T("%s(%04d年%02d月%02d日%02d時%02d分%02d秒)%s"),
fname, year, month, day, hh, mm, ss, ext);
if( File_Existence(fname_buf)==0 ){
fp = _tfopen(fname_buf, _T("wb") );
}
else{
return NULL;
}
//else{
// Sleep(1001);
// fp = now_date_fopen(fname, ext);
//}
return fp;
}
/***************************************************************
*
* ========目的========
* now_date_fopen()の高機能版。秒単位まで同じであっても、何とか連番を入れて
* 一意なファイル名を生成する。
* 基本はnow_date_fopenに準じる。
*
* ========引数========
* TCHAR *fname 基本ファイル名。拡張子は含まない。NULLでもよいがあまり意味がない。
* TCHAR *ext 拡張子。NULLでもよい。
*
* ========詳細========
* 基本はnow_date_fopenに準じる。
* 例えば、now_date_fopen("ログファイル", ".log"); として2021/11/12 13:31:22に呼び出せば、
* ログファイル(2021年11月12日13時31分22秒).log
* としてファイルが開かれる。重複があった場合は、数字を増やしながら
* ログファイル(2021年11月12日13時31分22秒)0001.log
* ログファイル(2021年11月12日13時31分22秒)0002.log
* ...
* というファイルが開かれる。
*
*
*/
FILE * now_date_fopen_n(TCHAR *fname, TCHAR *ext){
int year, month, day, hh, mm, ss;
TCHAR fname_buf[128*256];
FILE *fp=NULL;
int n;
get_yymmddhhmmss_int_now( &year, &month, &day, &hh, &mm, &ss );
_sntprintf(fname_buf, 128*256-2, _T("%s(%04d年%02d月%02d日%02d時%02d分%02d秒)%s"),
fname, year, month, day, hh, mm, ss, ext);
n = 1;
while( File_Existence(fname_buf)==1 ){
_sntprintf(fname_buf, 128*256-2, _T("%s(%04d年%02d月%02d日%02d時%02d分%02d秒)%04d%s"),
fname, year, month, day, hh, mm, ss, n, ext);
n++;
}
fp = _tfopen(fname_buf, _T("wb") );
return fp;
}
/***************************************************************
*
* ========目的========
* 呼び出された瞬間の日付を格納する。
*
*/
int get_yymmddhhmmss_int_now( int *year, int *month, int *day, int *hh, int *mm, int *ss ){
time_t time_t_strct;
struct tm *tm_strct;
if( year==NULL || month==NULL || day==NULL || hh==NULL || mm==NULL || ss==NULL ) return 0;
time( &time_t_strct );
tm_strct = localtime( &time_t_strct );
*year = tm_strct->tm_year + 1900;
*month = tm_strct->tm_mon + 1;
*day = tm_strct->tm_mday;
*hh = tm_strct->tm_hour;
*mm = tm_strct->tm_min;
*ss = tm_strct->tm_sec;
return 1;
}
/***************************************************************
*
* ========目的========
*
* 呼び出された瞬間の日時を付加した一意のファイル名を生成し、
* そのファイル名をもってファイルを書き込みでopenしファイルポインタを返す
* fnameには基本ファイル名を、extに拡張子指定する。拡張子にはドットを含むこと。
* fname、extは書き込み不可のリテラル文字列でかまわないし
* NULLでもよい。その場合は日時だけでファイル名が構成される。
*
* ========引数========
* TCHAR *fname 基本ファイル名。拡張子は含まない。NULLでもよいがあまり意味がない。
* TCHAR *ext 拡張子。NULLでもよい。
*
* ========詳細========
* マイクロ秒まで対応した関数にしたかったが、どうも環境依存が過ぎるのでやめた
*
*/
FILE * now_date_micro_fopen(TCHAR *fname, TCHAR *ext){
int year, month, day, hh, mm, ss, ms;
TCHAR fname_buf[128*256];
FILE *fp=NULL;
get_yymmddhhmmssms_int_now( &year, &month, &day, &hh, &mm, &ss, &ms );
_sntprintf(fname_buf, 128*256-2, _T("%s(%04d年%02d月%02d日%02d時%02d分%02d秒.%06d)%s"),
fname, year, month, day, hh, mm, ss, ms, ext);
if( File_Existence(fname_buf)==0 ){
fp = _tfopen(fname_buf, _T("wb") );
}
else{
return NULL;
}
//else{
// Sleep(1001);
// fp = now_date_fopen(fname, ext);
//}
return fp;
}
/***************************************************************
*
* ========目的========
* 呼び出された瞬間の日付を格納する。
*
* ========詳細========
* 時刻の正確な取得(マイクロ秒まで)について、色々やってみたのであるが
* そもそもwin732bit環境でclock_getres()が実装されていない可能性に思い至る。
* 本関数は廃止。まったく。
*/
int get_yymmddhhmmssms_int_now( int *year, int *month, int *day, int *hh, int *mm, int *ss, int *ms ){
//struct timeval time_ms;
struct timespec ts;
struct tm *tm_strct;
time_t time_t_strct;
if( year==NULL || month==NULL || day==NULL || hh==NULL || mm==NULL || ss==NULL || ms==NULL ) return 0;
//gettimeofday(&time_ms, NULL);
clock_getres(CLOCK_REALTIME, &ts);
//tm_strct = localtime(&time_ms.tv_sec);//tv_sec に入っている値が、通常time( &time_t );として取得されるtime_tの値
//tm_strct = localtime( (time_t*)(&(ts.tv_sec)) );//tv_sec に入っている値が、通常time( &time_t );として取得されるtime_tの値
//localtime_r( &ts.tv_sec, &tm_strct);
time( &time_t_strct );
_tprintf( _T("time_t_strct=%ld\n"), time_t_strct );
_tprintf( _T("ts.tv_sec =%ld\n"), ts.tv_sec );
_tprintf( _T("ts.tv_nsec =%ld\n"), ts.tv_nsec );
time_t_strct = ts.tv_sec;
tm_strct = localtime( &time_t_strct );
_tprintf( _T("tm_strct->tm_year=%d\n"), tm_strct->tm_year );
*year = tm_strct->tm_year + 1900;
*month = tm_strct->tm_mon + 1;
*day = tm_strct->tm_mday;
*hh = tm_strct->tm_hour;
*mm = tm_strct->tm_min;
*ss = tm_strct->tm_sec;
*ms = ts.tv_nsec;//マイクロ秒
_tprintf( _T("%04d年%02d月%02d日%02d時%02d分%02d秒.%06d\n"), *year, *month, *day, *hh, *mm, *ss, *ms );
return 1;
}
/***************************************************************
*
* ========目的========
* 指定されたファイル名があるかどうか確認する。
* ファイルがあれば1、なければ0を返す。
*
*/
int File_Existence(TCHAR *fname){
int handle, re;
struct _tfinddatai64_t find;
if( fname==NULL ) return 0;
re = 0;
handle = _tfindfirsti64(fname, &(find));
if(handle>=0) re = 1;
_findclose(handle);
return re;
}
/***************************************************************
*
* ========目的========
* グループの状態を出力する。
*
*
*/
int fprintf_dailycalendar(FILE *fp, DAILYCALENDAR *dailycalendar){
if( fp==NULL || dailycalendar==NULL ) return 0;
_ftprintf(fp, _T("%d日目\t"), dailycalendar->day );
fprintf_dailycalendar_flock(fp, dailycalendar->today);
_ftprintf(fp, _T("\r\n") );
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* 全人口の状態を出力する。
*
*
*/
int fprintf_dailycalendar_flock(FILE *fp, FLOCK *flock){
int i;
if( fp==NULL || flock==NULL ) return 0;
_ftprintf(fp, _T("%d人\t"), flock->n );
for(i=0;i<flock->n;i++){
fprintf_dailycalendar_humans(fp, &(flock->humans[i]) );
_ftprintf(fp, _T("\t") );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* 全人口の状態を出力する。
*
*
*/
int fprintf_dailycalendar_humans(FILE *fp, HUMANS *humans){
int x;
if( fp==NULL || humans==NULL ) return 0;
_ftprintf(fp, _T("ID%04d"), humans->no );
if( humans->condition==CONDITION_DEATH ){//死亡
_ftprintf(fp, _T("死亡") );
}
else if( humans->virusn>0 ){//感染済み
for(x=0;x<VIRUS_N;x++){
if( humans->sars_v[x].type==VIRUS_1 ) _ftprintf(fp, _T("V1") );
if( humans->sars_v[x].type==VIRUS_2 ) _ftprintf(fp, _T("V2") );
if( humans->sars_v[x].type==VIRUS_3 ) _ftprintf(fp, _T("V3") );
if( humans->sars_v[x].type==VIRUS_4 ) _ftprintf(fp, _T("V4") );
}
_ftprintf(fp, _T("感染") );
}
else{
_ftprintf(fp, _T("健康") );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayの細かい統計情報を記録する
*
*
*/
int fprintf_daily_statistics(FILE *fp, DAILYCALENDAR *dailycalendar){
STATISTICS statistics;
int x;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("[感染統計情報]\r\n") );
_ftprintf(fp, _T("感染%04d日目\r\n"), statistics.day );
_ftprintf(fp, _T("%d人中\r\n"), statistics.n );
_ftprintf(fp, _T("健康者 \t%4d\t人\r\n"), statistics.living );
_ftprintf(fp, _T("患者 \t%4d\t人\r\n"), statistics.patient );
_ftprintf(fp, _T("死者 \t%4d\t人\r\n"), statistics.dead );
_ftprintf(fp, _T("新規感染\t%4d\t人\r\n"), statistics.infect_new );
_ftprintf(fp, _T("新規発症\t%4d\t人\r\n"), statistics.onset_new );
for(x=0;x<VIRUS_N;x++){
_ftprintf(fp, _T("ウイルスタイプ%02d\t%4d\t人感染中\r\n"), x+1, statistics.sars_v[x] );
}
_ftprintf(fp, _T("感染者総数 \t%4d\t人\r\n"), statistics.infect_total );
for(x=0;x<VIRUS_N;x++){
_ftprintf(fp, _T("ウイルスタイプ%02d\t%4d\t人感染済\r\n"), x+1, statistics.sars_v_total[x] );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayの感染者各個人の細かい統計情報を記録する
*
*
*/
int fprintf_humans_statistics(FILE *fp, DAILYCALENDAR *dailycalendar){
int i, x, count;
FLOCK *flock;
if( fp==NULL || dailycalendar==NULL ) return 0;
_ftprintf(fp, _T("\t[感染者詳細情報]\r\n") );
count = 0;
flock = dailycalendar->today;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
i = 0;
while( find_infection( flock, &i ) ){
count++;
_ftprintf(fp, _T("\t::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("\t感染者%d人目\t%04dさん"), count, flock->humans[i].no );
if( flock->humans[i].condition==CONDITION_INCUB )_ftprintf(fp, _T("\t発症中") );
if( flock->humans[i].condition==CONDITION_SEVEE )_ftprintf(fp, _T("\t重症") );
if( flock->humans[i].condition==CONDITION_DEATH )_ftprintf(fp, _T("\t死亡") );
_ftprintf(fp, _T("\r\n"), count, flock->humans[i].no );
_ftprintf(fp, _T("\t感染%02d日目\r\n"), flock->humans[i].count );
for(x=0;x<VIRUS_N;x++){
if( flock->humans[i].sars_v[x].Population > 0 ){
_ftprintf(fp, _T("\tウイルスタイプ%02d\t%lld粒 Rt=%f\r\n"), flock->humans[i].sars_v[x].type, flock->humans[i].sars_v[x].Population, flock->humans[i].sars_v[x].Rt );
}
}
i++;
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayの感染に関する統計情報を取得する
*
*
*/
int get_statistics_flock(DAILYCALENDAR *dailycalendar, STATISTICS *statistics){
int i, x, flg;
FLOCK *flock;
if( dailycalendar==NULL || statistics==NULL ) return 0;
statistics->day = dailycalendar->day;
statistics->n = dailycalendar->today->n;//グループの人数
statistics->living = 0;//健康な生者の人数
statistics->patient = 0;//患者の人数
statistics->infect_new = 0;
statistics->onset_new = 0;
statistics->dead = 0;//死者数
statistics->infect_total = 0;//これまで感染した人の総計
for(x=0;x<VIRUS_N;x++){
statistics->sars_v [ x ] = 0;
statistics->sars_v_total[ x ] = 0;
}
flock = dailycalendar->today;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
for(i=0;i<flock->n;i++){
if( flock->humans[i].condition==CONDITION_ALIVE ){//健康に動いて生きている
if( flock->humans[i].virusn>0 ){//未発症感染者
(statistics->patient)++;
}
else{//生きている健康者
(statistics->living)++;
}
}
else if( flock->humans[i].condition==CONDITION_INCUB ){//発症している
(statistics->patient)++;
if( flock->humans[i].count==TIME_LIMIT_INCUBATION ) (statistics->onset_new)++;//新規発症者数
}
else if( flock->humans[i].condition==CONDITION_SEVEE ){//重症
(statistics->patient)++;
}
else if( flock->humans[i].condition==CONDITION_DEATH ){//死んでいる
(statistics->dead)++;
}
if( flock->humans[i].count==0 ) (statistics->infect_new)++;//新規感染者数
if( flock->humans[i].virusn>0 ){//現時点でのウイルス感染数カウント
for(x=0;x<VIRUS_N;x++){
if( flock->humans[i].sars_v[x].Population > 0 ){
(statistics->sars_v[ flock->humans[i].sars_v[x].type-1 ])++;
}
}
}
flg = 0;
for(x=0;x<VIRUS_N;x++){// 感染履歴のカウント
if( flock->humans[i].sars_v[x].type != VIRUS_0 ){
(statistics->sars_v_total[ flock->humans[i].sars_v[x].type-1 ])++;
flg = 1;
}
}
if( flg == 1 ){// 一度でも感染したことがある
(statistics->infect_total)++;
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* todayのウイルスごとの度数分布を出力する
*
*
*/
int fprintf_daily_virus_Frequency_distribution(FILE *fp, DAILYCALENDAR *dailycalendar){
int i, x;
FLOCK *flock;
if( fp==NULL || dailycalendar==NULL ) return 0;
_ftprintf(fp, _T("%04d日目|"), dailycalendar->day );
flock = dailycalendar->today;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
i = 0;
while( find_infection( flock, &i ) ){
for(x=0;x<VIRUS_N;x++){
if( flock->humans[i].sars_v[x].Population > 0 ){
_ftprintf(fp, _T("%d"), flock->humans[i].sars_v[x].type );
}
}
i++;
}
_ftprintf(fp, _T("\r\n") );
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayの最終的な細かい統計情報を記録する
* ウイルスの感染履歴や免疫の有無も調べる
* 最終的な結果を算出するのが目的なので、fprintf_daily_statistics()と異なり、感染歴を勘定する
*/
int fprintf_daily_statistics2(FILE *fp, DAILYCALENDAR *dailycalendar){
STATISTICS statistics;
int x;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("[感染統計最終結果]\r\n") );
_ftprintf(fp, _T("感染%04d日目までシミュレーション\r\n"), statistics.day );
_ftprintf(fp, _T("%d人中\r\n"), statistics.n );
_ftprintf(fp, _T("生存者 \t%4d\t人\r\n"), statistics.living );
_ftprintf(fp, _T("患者 \t%4d\t人\r\n"), statistics.infect_total );
_ftprintf(fp, _T("死者 \t%4d\t人\r\n"), statistics.dead );
_ftprintf(fp, _T("新規感染\t%4d\t人\r\n"), statistics.infect_new );
_ftprintf(fp, _T("新規発症\t%4d\t人\r\n"), statistics.onset_new );
for(x=0;x<VIRUS_N;x++){
_ftprintf(fp, _T("ウイルスタイプ%02d\t%4d\t人感染\r\n"), x+1, statistics.sars_v_total[x] );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayの最終的な細かい統計情報を記録する
* データ検証のために単一のデータのみ出力する
*/
int fprintf_daily_statistics3(FILE *fp, DAILYCALENDAR *dailycalendar){
STATISTICS statistics;
//int x;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
//_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
//_ftprintf(fp, _T("[感染統計最終結果]\r\n") );
_ftprintf(fp, _T("感染%04d日目\t"), statistics.day );
//_ftprintf(fp, _T("%d人中\r\n"), statistics.n );
//_ftprintf(fp, _T("生存者\t%4d\t人\r\n"), statistics.living );
//_ftprintf(fp, _T("患者 \t%4d\t人\r\n"), statistics.patient );
//_ftprintf(fp, _T("死者 \t%4d\t人\r\n"), statistics.dead );
//_ftprintf(fp, _T("新規感染\t%4d\t人\r\n"), statistics.infect_new );
_ftprintf(fp, _T("新規発症\t%4d\t人\r\n"), statistics.onset_new );
//for(x=0;x<VIRUS_N;x++){
// _ftprintf(fp, _T("ウイルスタイプ%02d\t%4d\t人感染\r\n"), x+1, statistics.sars_v_total[x] );
//}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
*
* 度数分布を出力するために、各ウイルス数へのインデックスを生成する。
* ソートと名づけられているが、ソート順にアクセスするためのインデックスを生成する
* 間接ソート処理である。
* インデックス生成をするだけなので、元データは(当然)変化させない
* 本関数は、安定ソートである。
*
* ========引数========
* int *index ソート対象。この値でソート順アクセスする。
* int *sars_v 各ウイルス種ごとに感染者数を格納している。
* int len ソート対象の要素数。実質的にはVIRUS_Nだが、抽象化のためあえて引数化。
*
* ========詳細========
* statistics 構造体のsars_v メンバはsars_v[ VIRUS_N ]として定義されており、先頭から各ウイルスの感染者数を格納している。
* このデータは先頭からウイルス種ごとに並べられているという前提があるため、本関数で勝手に並べ替えることはできない。
* しかし、fprintf_daily_virus_Frequency_distribution_Sorting()は出現頻度の高いウイルスをまとめて、降順に
* 度数分布表示する関数である。
* 度数分布のデータを並べ替え表示するため、出現頻度の高いウイルスの感染者データから順にアクセスするための
* インデックスを生成するのが本関数の主旨である。
* 関数の大要はインデックスのインデックスを生成する処理となりちょっと想像しづらいので、具体例を示す。
* VIRUS_N = 4
* sars_v[0] ウイルス初期型 感染者 100人
* sars_v[1] ウイルス変異型1 感染者 50人
* sars_v[2] ウイルス変異型2 感染者 250人
* sars_v[3] ウイルス変異型3 感染者 50人
* このとき、降順に並べるのであれば、これらのデータにアクセスするためのインデックスは
* index[0] = 2
* index[1] = 0
* index[2] = 1
* index[3] = 3
* である。これを先頭から順にsars_vデータの添字として使えば
* index[0] = 2 → sars_v[ index[0] ] = 250
* index[1] = 0 → sars_v[ index[1] ] = 100
* index[2] = 1 → sars_v[ index[2] ] = 50
* index[3] = 3 → sars_v[ index[3] ] = 50
* となり、indexを0から順に添字として使うことで自動的に降順アクセスとなる。
* indexは、それぞれのデータの多い順ランキング(0が最頻)と言い換えることもできる。
*
*/
int bubbleSort(int *index, int *sars_v, int len){
int i, j, tmp;
if( index==NULL || sars_v==NULL ) return 0;
for(i=0;i<len;i++){
index[i] = i;
}
for(i=0; i<len; i++){
for(j=len-1; j>i; j--){
if( sars_v[ index[j] ] > sars_v[ index[j-1] ] ){
tmp = index[j];
index[j] = index[j-1];
index[j-1] = tmp;
}
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* todayのウイルスごとの度数分布を出力する
* 最頻値を左寄せ。
* 頻度が同じなら数値の若い順
*
*/
int fprintf_daily_virus_Frequency_distribution_Sorting(FILE *fp, DAILYCALENDAR *dailycalendar){
int i, x;
int *index;
STATISTICS statistics;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
_ftprintf(fp, _T("%04d日目|"), dailycalendar->day );
//i = 0;
//flock = dailycalendar->today;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
index = malloc( sizeof(int)*VIRUS_N );
bubbleSort( index, statistics.sars_v, VIRUS_N);
for(x=0;x<VIRUS_N;x++){
for(i=0;i<statistics.sars_v[ index[x] ];i++){
_ftprintf(fp, _T("%d"), index[x]+1 );
}
}
_ftprintf(fp, _T("\r\n") );
fflush( fp );
free( index );
return 1;
}
/***************************************************************
*
* ========目的========
* todayのグラフ用CSVデータを出力する
*
*
*/
int fprintf_daily_graph(FILE *fp, DAILYCALENDAR *dailycalendar, int world_line){
int x;
STATISTICS statistics;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
_ftprintf(fp, _T("感染者数\t%04d日目\t%d\t%d\r\n"), statistics.day, world_line, statistics.patient );
_ftprintf(fp, _T("死者数\t%04d日目\t%d\t%d\r\n"), statistics.day, world_line, statistics.dead );
for(x=0;x<VIRUS_N;x++){
_ftprintf(fp, _T("ウイルス(%02d)感染者数\t%04d日目\t%d\t%4d\r\n"), x+1, statistics.day, world_line, statistics.sars_v[x] );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* todayのグラフ用CSVデータを出力する
* デバッグ用に、感染者総数を出力する。
*
*/
int fprintf_daily_graph2(FILE *fp, DAILYCALENDAR *dailycalendar, int world_line){
int x;
STATISTICS statistics;
if( fp==NULL || dailycalendar==NULL ) return 0;
if( get_statistics_flock(dailycalendar, &statistics)==0 ){// 統計データ取得失敗
_tprintf( _T("統計データ取得失敗\n") );
return 0;
}
_ftprintf(fp, _T("感染者総数\t%04d日目\t%d\t%d\r\n"), statistics.day, world_line, statistics.infect_total );
for(x=0;x<VIRUS_N;x++){
_ftprintf(fp, _T("ウイルス(%02d)感染者総数\t%04d日目\t%d\t%4d\r\n"), x+1, statistics.day, world_line, statistics.sars_v_total[x] );
}
fflush( fp );
return 1;
}
/***************************************************************
*
* ========目的========
* DailyCalendarを作成する際に用いられた設定データを出力する
*
*/
int fprintf_daily_Calendar_setting(FILE *fp, int world_line, DAILYCALENDAR *dailycalendar){
ACTIONS *actions;
double risk;
if( fp==NULL ) return 0;
actions = &(dailycalendar->actions);
_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("この世界での設定\r\n::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("%03d回目の世界線\r\n"), world_line );
_ftprintf(fp, _T("全人口%d\r\n"), dailycalendar->today->n );
_ftprintf(fp, _T("基本的なウイルスのRt=%f\r\n"), REPTNUM_VIRUS_1 );
_ftprintf(fp, _T("タイプ４ウイルスのRt=%f\r\n"), REPTNUM_VIRUS_4 );
_ftprintf(fp, _T("変異率=1/%d\r\n"), MUTATION_PROBABILITY );
_ftprintf(fp, _T("変異が起こるタイムリミットは%d日\r\n"), MUTATION_DAY_LIMIT );
_ftprintf(fp, _T("マスク装着による感染元人間からのウイルス数は%f倍に減少\r\n"), MASK_DEFENSE );
_ftprintf(fp, _T("マスク装着による感染先人間へ のウイルス数は%f倍に減少\r\n"), MASK_WEAKNESS );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%以上が感染者となった時点で発令\r\n"), (double)(actions->self_restraint_limit_m)/PER_DEF*100 );
_ftprintf(fp, _T("緊急事態宣言は全人口の%d人以上が感染者となった時点で発令\r\n"), (dailycalendar->today->n)*(actions->self_restraint_limit_m)/PER_DEF );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%未満にまでが感染が減少した時点で解除\r\n"), (double)(actions->self_restraint_limit_m)/PER_DEF*100/2 );
_ftprintf(fp, _T("緊急事態宣言は全人口の%d人未満にまでが感染が減少した時点で解除\r\n"), (dailycalendar->today->n)*(actions->self_restraint_limit_m)/PER_DEF/2 );
_ftprintf(fp, _T("緊急事態宣言が発令されると全人口の%f%%が行動を自粛\r\n"), (double)(actions->self_restraint_m)/PER_DEF*100 );
_ftprintf(fp, _T("行動を自粛した場合、感染確率は%f倍に減少\r\n"), EMERGENCY_RATIO );
_ftprintf(fp, _T("クラスター数%d\r\n"), actions->modulo_M );
_ftprintf(fp, _T("クラスターサイズ%d\r\n"), (dailycalendar->today->n)/(actions->modulo_M) );
_ftprintf(fp, _T("クラスター間感染率%f%%\r\n"), (double)(actions->mod_risk_wait_m)/PER_DEF*100 );
risk = (dailycalendar->today->n)/(actions->modulo_M)*(dailycalendar->today->n)/(actions->modulo_M)*(double)(actions->mod_risk_wait_m)/PER_DEF;
_ftprintf(fp, _T("2クラスター間感染リスク%fp\r\n"), risk );
_ftprintf(fp, _T("医療崩壊ラインは%d人を超えたとき\r\n"), (dailycalendar->today->n)*MEDICALCAPACITYOVER/PER_DEF );
return 1;
}
/***************************************************************
*
* ========目的========
* DailyCalendarを作成する際に用いられた設定データを出力する
* ACTIONSしか引数にできないとき用
*/
int fprintf_daily_Calendar_setting_act(FILE *fp, int world_line, ACTIONS *actions, int n){
double risk;
if( fp==NULL ) return 0;
_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("この世界での設定\r\n::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("%03d回目の世界線\r\n"), world_line );
_ftprintf(fp, _T("全人口%d\r\n"), n );
_ftprintf(fp, _T("基本的なウイルスのRt=%f\r\n"), REPTNUM_VIRUS_1 );
_ftprintf(fp, _T("タイプ４ウイルスのRt=%f\r\n"), REPTNUM_VIRUS_4 );
_ftprintf(fp, _T("変異率=1/%d\r\n"), MUTATION_PROBABILITY );
_ftprintf(fp, _T("変異が起こるタイムリミットは%d日\r\n"), MUTATION_DAY_LIMIT );
_ftprintf(fp, _T("マスク装着による感染元人間からのウイルス数は%f倍に減少\r\n"), MASK_DEFENSE );
_ftprintf(fp, _T("マスク装着による感染先人間へ のウイルス数は%f倍に減少\r\n"), MASK_WEAKNESS );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%以上が感染者となった時点で発令\r\n"), (double)(actions->self_restraint_limit_m)/PER_DEF*100 );
_ftprintf(fp, _T("緊急事態宣言は全人口の%d人以上が感染者となった時点で発令\r\n"), n*(actions->self_restraint_limit_m)/PER_DEF );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%未満にまでが感染が減少した時点で解除\r\n"), (double)(actions->self_restraint_limit_m)/PER_DEF*100/2 );
_ftprintf(fp, _T("緊急事態宣言は全人口の%d人未満にまでが感染が減少した時点で解除\r\n"), n*(actions->self_restraint_limit_m)/PER_DEF/2 );
_ftprintf(fp, _T("緊急事態宣言が発令されると全人口の%f%%が行動を自粛\r\n"), (double)(actions->self_restraint_m)/PER_DEF*100 );
_ftprintf(fp, _T("行動を自粛した場合、感染確率は%f倍に減少\r\n"), EMERGENCY_RATIO );
_ftprintf(fp, _T("クラスター数%d\r\n"), actions->modulo_M );
_ftprintf(fp, _T("クラスターサイズ%d\r\n"), n/(actions->modulo_M) );
_ftprintf(fp, _T("クラスター間感染率%f%%\r\n"), (double)(actions->mod_risk_wait_m)/PER_DEF*100 );
risk = (n)/(actions->modulo_M)*(n)/(actions->modulo_M)*(double)(actions->mod_risk_wait_m)/PER_DEF;
_ftprintf(fp, _T("2クラスター間感染リスク%fp\r\n"), risk );
_ftprintf(fp, _T("医療崩壊ラインは%d人を超えたとき\r\n"), (n*MEDICALCAPACITYOVER/PER_DEF) );
return 1;
}
/***************************************************************
*
* ========目的========
* DailyCalendarを作成する際に用いられた設定データを出力する
*
*/
int fprintf_daily_Calendar_setting2(FILE *fp, int world_line){
if( fp==NULL ) return 0;
_ftprintf(fp, _T("::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("この世界での設定\r\n::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::\r\n") );
_ftprintf(fp, _T("%03d回目の世界線\r\n"), world_line );
_ftprintf(fp, _T("基本的なウイルスのRt=%f\r\n"), REPTNUM_VIRUS_1 );
_ftprintf(fp, _T("タイプ４ウイルスのRt=%f\r\n"), REPTNUM_VIRUS_4 );
_ftprintf(fp, _T("変異率=1/%d\r\n"), MUTATION_PROBABILITY );
_ftprintf(fp, _T("変異が起こるタイムリミットは%d日\r\n"), MUTATION_DAY_LIMIT );
_ftprintf(fp, _T("マスク装着による感染元人間からのウイルス数は%f倍に減少\r\n"), MASK_DEFENSE );
_ftprintf(fp, _T("マスク装着による感染先人間へ のウイルス数は%f倍に減少\r\n"), MASK_WEAKNESS );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%以上が感染者となった時点で発令\r\n"), (double)(ACT_SELF_RESTRAINT_LIMIT_M)/PER_DEF*100 );
_ftprintf(fp, _T("緊急事態宣言は全人口の%f%%未満にまでが感染が減少した時点で解除\r\n"), (double)(ACT_SELF_RESTRAINT_LIMIT_M)/PER_DEF*100/2 );
_ftprintf(fp, _T("緊急事態宣言が発令されると全人口の%f%%が行動を自粛\r\n"), (double)ACT_SELF_RESTRAINT_M/PER_DEF*100 );
_ftprintf(fp, _T("行動を自粛した場合、感染確率は%f倍に減少\r\n"), EMERGENCY_RATIO );
_ftprintf(fp, _T("さらにマスクをつける素質のある人はマスクを付けるようになる\r\n") );
return 1;
}
/***************************************************************
*
* ========目的========
* シミュレーションに時間がかかりすぎる場合、この関数を有効にして
* 寝ている間に演算実行、
* 終了後一定時間でPCをシャットダウンするのもアリ。
*
*/
int myPCshutdown(void){
int res;
TCHAR timer[3];
res = _tsystem( _T("C:\\WINDOWS\\system32\\shutdown.exe -s -t 240 -c \"windowsを終了します\"") );
if(res){;}
_tprintf( _T("シャットダウンをやめるなら1を入力\n") );
while(1){
_fgetts(timer, 2, stdin);
if(timer[0]=='1'){
res = _tsystem( _T("C:\\WINDOWS\\system32\\shutdown.exe -a") );
break;
}
Sleep(100);
}
if(res){;}
return 1;
}
/***************************************************************
*
* ========目的========
* 感染を起こす。感染において、再生産数は少数だが、感染の有無は1か0の真偽値である。
* 単純な再生産数ではもちろん少数を決定的に反映できないので、
* 毎回再生産数を元にした乱数抜き出しを感染の是非に用いる。
*
*
* ========引数========
*
* ========詳細========
* 本関数は、感染における確率浮動を適切に反映するために実装されている。
* 例えば、再生産数が2であれば、感染は2人に対して起こる。
* しかし、実際には感染は確率的であり、「大数の法則が現れる」ような大量感染では平均2人に感染
* するといえるが、全ての事例で確定的に2人に感染するわけではない。そもそもそのような確定的過ぎる
* 実装は感染シミュレーションとしてはいささか問題である(そのような実装は方程式を解くのと変わらない)。
* また、再生産数は整数ではなく、本質的には小数・分数である。再生産数が1.4の場合、例えば100人の感染源
* から140人に感染するような確率であるが、全体を通してこれを調整するのはちょっとめんどうである。
* この問題を解決するため、本関数では以下の実装を採用している。
* ・感染機会の細粒化
* ・感染毎での確率浮動
* 具体的には、次のように計算する。
* 実効再生産数Rt0が1.4の場合、感染機会(期待値)を0.5回ずつに細分する。
* 1回目、0.5回分の感染シミュレート。乱数(さいころ)で5/10の確率で感染を実行。これは、期待値を0.5消費してくじを引くようなもの。残り0.9
* 2回目、0.5回分の感染シミュレート。乱数(さいころ)で5/10の確率で感染を実行。残り0.4
* 3回目、0.4回分の感染シミュレート。乱数(さいころ)で4/10の確率で感染を実行。残り0.0
* 言い換えると、実効再生産数を期待値の引換券と捉え、さいころを振るごとに消費していく感覚である。
* このように、消費する期待値に応じた乱数さいころを振り、その確率的な結果を元に感染を実行するのである。
* 毎回このように乱数で決めていけば、全体としては大数の法則に従うが、毎回かなりのばらつきが出るであろう。
* この方法においては、細粒度の細かさも大切である。あまり細かくすると、一回の感染につき大数の法則が
* 適用されすぎてしまい、かえって確定的になりかねない。また、計算負荷も大きくなりすぎる。
* ほどほどにバラける細度がベストである。今回は、1人への感染を「確率50%を2回」の期待値に分解することにする。
* もちろん、端数のときは端数に対応した確率をかける。
* 期待値0.4なら、10個中4個が当たり(感染成立)、6個がはずれ(感染不成立)の目をもつさいころを振るのと同義である。
*
*/
int virus_infection(HUMANS *humans, SARS_V *sars_v){
if( humans && sars_v ){;}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分の感染を一通り起こす
* 本関数は、実効再生産数を感染可能になった初日に全て使う。
* デバッグや感染状況の解析を行いやすい反面、感染のバラつきに関しては不十分なので
* 本格的なシミュレーション時にはdaily_infection2()が望ましい。
*
* ========引数========
*
* ========詳細========
* 本関数は感染可能初日に一気に次の感染者を発生させるので、感染拡大が5日周期という周期性を持つようになる。
* これは感染拡大の要因推測や数値解析しやすい反面、ランダム性に欠け、一部の感染要因を的確に反映しないかもしれない。
* 特に、感染源となる人間の体内では、感染初期ほど変異ウイルス数の絶対量が少なく、極めて不利である。
* そのため、本関数では変異ウイルスが感染元から感染先へ継代される可能性が極めて低い。
* また、この仕様は、Rt数が正確に反映されているかどうか試算しやすい反面、ばらつきに欠ける。
* 即ち、実効再生産数Rtが2.4であれば、感染可能になった初日に(初期設定では)感染確率0.5の感染試行を4回、
* 感染確率0.4の感染試行を1回行う。
*
* irand()は擬似乱数でしかも周期が短い乱数とされているが、本実装では毎回感染者選定などのダイス転がし
* に毎回一回きり使うだけなので、周期が短くても十分なバラつきが得られると仮定している。
* targetの選定は、感染源となる人間の添字 i からどれだけ離れているかを乱数で規定し、剰余で後方参照しつつ
* 自身はtargetにはならない。例えば、
* iが2で全住民数が5のとき、targetは
* 2 + 1 の3
* 2 + 2 の4
* 2 + 3 の0
* 2 + 4 の1
* の4パターンである。それぞれ0～4までの乱数に+1したものをiにプラスし、住民数5で剰余を得ている
*
*
*
*/
int daily_infection(DAILYCALENDAR *dailycalendar){
int i, target;
FLOCK *flock;
SARS_V sars_v[ VIRUS_N ];
if( dailycalendar==NULL ) return 0;
//_tprintf( _T("\r%d日目の感染"), dailycalendar->day );
flock = dailycalendar->yesterday;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
i = 0;
while( find_infection( flock, &i ) ){//感染源となる人間を探す
//_tprintf( _T("感染者発見id=%d\n"), i );
daily_infection_sub_viruscopy( flock->humans[i].sars_v, sars_v, flock->humans[i].mask );
if( dailycalendar->emergency_flg == EMERGENCY_FLG_ON &&
flock->humans[i].self_restraint==SELF_RESTRAINT_ON ){// 緊急事態宣言下で、かつ、感染元に自粛する意思がある
daily_infection_sub_Rt_n( &(dailycalendar->today->humans[i]), EMERGENCY_RATIO );
}
if( is_Virus_Rt(sars_v)==1 && is_Virus_Population(sars_v)==1 && flock->humans[i].count >= TIME_LIMIT_CAN_INFECT ){
//daily_infection_sub_viruscopy( flock->humans[i].sars_v, sars_v, flock->humans[i].mask );
//if( is_Virus_Rt(sars_v)==1 ) daily_infection_sub_Rt_0( &(dailycalendar->today->humans[i]) );
while( is_Virus_Rt(sars_v)==1 ){// ウイルス群集のデータどれか一種類でも、実効再生産数が残っている
target = (i + irand( flock->n-1 ) + 1) % ( flock->n );//この計算は、感染源からの配列距離を算出するもの
daily_infection_sub_virus_invasion( &(dailycalendar->today->humans[ target ]), sars_v );
}
daily_infection_sub_Rt_0( &(dailycalendar->today->humans[i]) );
//_tprintf( _T("感染元%dのRtを消去\n"), i );
}
i++;
}
return 1;
}
// daily_infection()の下請け。この一日分の感染を計算するために、実効再生産数を含めたウイルス情報をコピーする
// マスクをつけていると排出ウイルス量が0.8倍される
int daily_infection_sub_viruscopy(SARS_V *sars_v1, SARS_V *sars_v2, int mask){
int x;
for(x=0;x<VIRUS_N;x++){
sars_v2[x].type = sars_v1[x].type;
sars_v2[x].Rt = sars_v1[x].Rt;
if( mask==MASK_ON ) sars_v2[x].Population = sars_v1[x].Population * MASK_WEAKNESS;
else sars_v2[x].Population = sars_v1[x].Population;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分の感染を一通り起こす
* 本関数は、実効再生産数を感染可能になった日から一単位ずつ使って感染させていく。
* 感染元が実効再生産数 Rt=2.4のウイルスに感染している場合、Rtを回数券のように毎日一枚ずつ使っていく。
* 感染拡大によりよくランダム性が反映されるので、実際のシミュレーション時にはこちらがオススメ
*
* ========引数========
*
* ========詳細========
* 本関数は毎日少しずつ感染者を増やしていくのでランダム性がよい。
* 特に、変異ウイルスの動向をよりよく反映すると思われる。
* 感染源となる人間の体内では、感染初期ほど変異ウイルス数の絶対量が少なく、極めて不利である。
* 変異ウイルスが発生した人間、または変異ウイルスに感染した人間は、感染初期であれば変異ウイルスの
* 個体数が少ないのでなかなか感染にまで至らず、変異ウイルスの人間から人間への継代に失敗する可能性が高い。
* 感染の態様としては、
* 実効再生産数Rtが2.4であれば、感染可能になった日から(初期設定では)感染確率0.5の感染試行を1回ずつ行っていく
* 実装における初期設定では人人感染を起こすときの確率試行の細かさPROBABILITY_PARTICLEが0.5として定義されている。
* この場合、感染から5日目に感染可能になったとすると
* 5日目 感染確率0.5の感染試行1回(Rtの残数1.9)
* 6日目 感染確率0.5の感染試行1回(Rtの残数1.4)
* 7日目 感染確率0.5の感染試行1回(Rtの残数0.9)
* 8日目 感染確率0.5の感染試行1回(Rtの残数0.4)
* 9日目 感染確率0.4の感染試行1回(Rtの残数0.0)
* という形で5日間感染試行を続ける。
*
* rand()は擬似乱数でしかも周期が短い乱数とされているが、本実装では毎回感染者選定などのダイス転がし
* に毎回一回きり使うだけなので、周期が短くても十分なバラつきが得られると仮定している。
* targetの選定は、感染源となる人間の添字 i からどれだけ離れているかを乱数で規定し、剰余で後方参照しつつ
* 自身はtargetにはならない。例えば、
* iが2で全住民数が5のとき、targetは
* 2 + 1 の3
* 2 + 2 の4
* 2 + 3 の0
* 2 + 4 の1
* の4パターンである。それぞれ0～4までの乱数に+1したものをiにプラスし、住民数5で剰余を得ている
*
*
*
*/
int daily_infection2(DAILYCALENDAR *dailycalendar){
int i, target, mod, mod_wait;
FLOCK *flock;
SARS_V sars_v[ VIRUS_N ];
if( dailycalendar==NULL ) return 0;
//_tprintf( _T("\r%d日目の感染"), dailycalendar->day );
flock = dailycalendar->yesterday;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
mod = dailycalendar->actions.modulo_M;
mod_wait = dailycalendar->actions.mod_risk_wait_m;
i = 0;
while( find_infection( flock, &i ) ){//感染源となる人間を探す
//_tprintf( _T("感染者発見id=%d\n"), i );
daily_infection2_sub_viruscopy( flock->humans[i].sars_v, sars_v, flock->humans[i].mask );
if( dailycalendar->emergency_flg == EMERGENCY_FLG_ON &&
flock->humans[i].self_restraint==SELF_RESTRAINT_ON ){// 緊急事態宣言下で、かつ、感染元に自粛する意思がある
daily_infection_sub_Rt_n( &(dailycalendar->today->humans[i]), EMERGENCY_RATIO );
}
if( is_Virus_Rt(sars_v)==1 /*&& is_Virus_Population(sars_v)==1*/ && flock->humans[i].count >= TIME_LIMIT_CAN_INFECT ){
//daily_infection_sub_viruscopy( flock->humans[i].sars_v, sars_v, flock->humans[i].mask );
//if( is_Virus_Rt(sars_v)==1 ) daily_infection_sub_Rt_0( &(dailycalendar->today->humans[i]) );
while( is_Virus_Rt(sars_v)==1 ){// ウイルス群集のデータどれか一種類でも、実効再生産数が残っている
target = (i + irand( flock->n-1 ) + 1) % ( flock->n );//この計算は、感染源からの配列距離を算出するもの
if( target%mod == i%mod ){//法が同じ
daily_infection_sub_virus_invasion( &(dailycalendar->today->humans[ target ]), sars_v );
}
else if( roll_the_dice(PER_DEF, mod_wait)==1 ){//法が異なるのでダイスロール
daily_infection_sub_virus_invasion( &(dailycalendar->today->humans[ target ]), sars_v );
}
else{//法が異なることにより感染に失敗したのでRtだけ減らす
daily_infection2_sub_Rt_Decrease_sars_v_n( sars_v );
}
}
daily_infection2_sub_Rt_Decrease_n( &(dailycalendar->today->humans[i]) );
//_tprintf( _T("感染元%dのRtを消去\n"), i );
}
i++;
}
return 1;
}
// daily_infection2()の下請け。この一日分の感染を計算するために、実効再生産数を含めたウイルス情報をコピーする
// マスクをつけていると排出ウイルス量が0.8倍される
// Rtは4.5を超えると使い切れない(感染後、約9日間ほどの感染期間があるので、0.5ずつだと使い切れない)
// oneday_ticketは、実効再生産数の最大値からコピー量を9で割って算出した値。
int daily_infection2_sub_viruscopy(SARS_V *sars_v1, SARS_V *sars_v2, int mask){
int x;
//_tprintf( _T("ONEDAY_TICKET=%f\n"), ONEDAY_TICKET );
for(x=0;x<VIRUS_N;x++){
sars_v2[x].type = sars_v1[x].type;
if( sars_v1[x].Rt >= ONEDAY_TICKET ){ sars_v2[x].Rt = ONEDAY_TICKET; }
else { sars_v2[x].Rt = sars_v1[x].Rt; }
if( mask==MASK_ON ){ sars_v2[x].Population = sars_v1[x].Population * MASK_WEAKNESS; }
else { sars_v2[x].Population = sars_v1[x].Population; }
//_tprintf( _T("sars_v2[x].Population=%lld\n"), sars_v2[x].Population );
}
return 1;
}
// daily_infection()の下請け。Rtを0.5引く
//
int daily_infection2_sub_Rt_Decrease_n(HUMANS *humans){
int x;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
if( humans->sars_v[x].Rt >= PROBABILITY_PARTICLE ) humans->sars_v[x].Rt = humans->sars_v[x].Rt - PROBABILITY_PARTICLE;
else humans->sars_v[x].Rt = 0;
}
return 1;
}
// daily_infection()の下請け。Rtを0.5引く
//
int daily_infection2_sub_Rt_Decrease_sars_v_n(SARS_V *sars_v){
int x;
if( sars_v==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
if( sars_v[x].Rt >= PROBABILITY_PARTICLE ) sars_v[x].Rt = sars_v[x].Rt - PROBABILITY_PARTICLE;
else sars_v[x].Rt = 0;
}
return 1;
}
/***************************************************************
*
* ========目的========
* daily_infection()の下請け。
* SARS_V で指定されたウイルスを感染させる。SARS_VのRtは書き換えられる。
* 各ウイルスのRtが0になるまで感染のチャレンジが行われる
* "*sars_v" は構造体配列の先頭ポインタではなく、構造体本体のポインタ。
*
* ========引数========
*
* ========詳細========
* 本関数で大切なのは、「変異ウイルスの放出過程を再現すること」である。
* ある人間の体内で、変異型2が発生したとする。
* 確率的に、変異は原種ウイルス数が増えた後に起こる可能性が高い。それだけチャンスが増えるからである。
* この点は変異ウイルスが外部に出る際に不利に働く。感染者の臨床において後半に出てきた変異ウイルスは
* 新参者であり、どうがんばっても体内での多数派にはなれない。
* これはつまり体外に出て新たな人間に感染するチャンスも低いことを意味する。
* つまり、「変異ウイルスは、変異した人間の体から出るときが最も競争率が高い」ということである。
* ハーディ・ワインベルグの法則のとおり、患者の飛沫に新参の変異ウイルスが入る可能性は低い。
* 一方で、感染先にもしも入ることができれば、原株も変異株も身一つで入植した新参であり、両方とも
* 個体数が拮抗し、次の感染時にはほぼ平等となるだろう。
* 「二人目からは平等」を実現するために、
* 1. ウイルス数に応じた感染チャンスのくじ引き
* 2. 感染が決まれば、どの株も同数が感染する
* として実装している。特に、2をせず、元ウイルス数に応じた感染数にすると、
* 感染が進んでも元のウイルス個体数が常に「ハーディワインベルグ比」として維持されることになり、
* 変異株が日の目を見る可能性は実質的にはゼロである。
*
*
*/
int daily_infection_sub_virus_invasion(HUMANS *humans, SARS_V *sars_v){
int x;
double Rt, dice;
long long max;
int denominator, numerator;
if( humans==NULL || sars_v==NULL ) return 0;
max = INVASION_LIMIT;
denominator = max/INVASION_LIMIT_FOOT_CUT;
for(x=0;x<VIRUS_N;x++){
//if( find_Virus_type( humans, sars_v[x].type )==0 ){//感染対象のhumansにtypeウイルスが未感染
if( sars_v[x].Rt > 0 ){//実効再生産数の処理が必要
//_tprintf( _T("sars_v[%d].Rt = %f\n"), x, sars_v[x].Rt );
if( sars_v[x].Rt < PROBABILITY_PARTICLE ){
Rt = sars_v[x].Rt;
sars_v[x].Rt = 0;
}
else{
Rt = PROBABILITY_PARTICLE;
sars_v[x].Rt = sars_v[x].Rt - PROBABILITY_PARTICLE;
}
//denominator = max/INVASION_LIMIT_FOOT_CUT;
numerator = sars_v[x].Population/INVASION_LIMIT_FOOT_CUT;
//_tprintf( _T("numerator=%d max=%lld denominator=%d sars_v[x].Population=%lld\n"), numerator, max, denominator, sars_v[x].Population );
if( numerator <= 0 ) continue;
if( sars_v[x].Population >= max || //感染に必要なウイルス数を満たしている
roll_the_dice( denominator, numerator) ==1 ){//ウイルス数は少ないけどダイスロールであたりを引いた
dice = (double)irand( 10 )/10;//0.0～0.9の値を返す。
if( dice < Rt ){
virus_invasion(humans, sars_v[x].type, Default_VIRUS_INVASION_N );
//_tprintf( _T("dice=%f 感染成立\n"), dice );
}
else{
//_tprintf( _T("dice=%f 感染不成立\n"), dice );
}
}
}
}
return 1;
}
// daily_infection()の下請け。Rtを0にする
//
int daily_infection_sub_Rt_0(HUMANS *humans){
int x;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
humans->sars_v[x].Rt = 0;
}
return 1;
}
// daily_infection()の下請け。Rtをn倍にする
//
int daily_infection_sub_Rt_n(HUMANS *humans, double n){
int x;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){
humans->sars_v[x].Rt = (humans->sars_v[x].Rt) * n;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分の突然変異を一通り起こす。
* 本関数では、既に感染が長期間に及んでいる人間の体内では変異が発生しないという実装である。
* これは必ずしも実態に沿った実装ではないが、日本のように重症者をほぼ隔離できる医療システムを維持している、
* あるいは長期感染者が、人ごみをうろうろするような世紀末状態でない国であれば概ね仮定しても良いと思われる。
*
* ========引数========
*
* ========詳細========
*
*/
int daily_mutation(DAILYCALENDAR *dailycalendar){
int i, mut, virs;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
// 突然変異処理は、「daily_infection()感染」の直前に行うので、yesterdayに対して行う
flock = dailycalendar->yesterday;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
i = 0;
while( find_infection( flock, &i ) ){//感染者を探す
if( flock->humans[i].count > MUTATION_DAY_LIMIT ) goto LABEL;
if( find_Virus_type( &(flock->humans[i]), VIRUS_2 ) && find_Virus_type( &(flock->humans[i]), VIRUS_3 ) ){
//_tprintf( _T("Warning!!!!!!!!!!!!W\n") );
//_tprintf( _T("%d日目にタイプ2とタイプ3の変異型ウイルスに同時感染し危険な変異型4発生\n"), dailycalendar->day );
virus_injection2( &(dailycalendar->today->humans[ i ]), VIRUS_4, 1);
}
mut = irand( MUTATION_PROBABILITY );
if( mut==MUTATION_PROBABILITY-1 ){//変異率の逆数の乱数(0～MUTATION_PROBABILITY-1)を取得し、変異判定
virs = irand( 2 );//変異率に対して、変異型1 OR 2のどちらになるかはさらに二択。
//_tprintf( _T("%d日目に変異型ウイルス%d発生\n"), dailycalendar->day, virs+2 );
virus_injection2( &(dailycalendar->today->humans[ i ]), virs+2, 1);
}
LABEL:
i++;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分の突然変異を一通り起こす
* 本関数は、重症者のような長期間感染し続けている患者の体内で変異が発生し、
* さらに市中に広がるような実装である。
* ウイルスが集団内でどのように振舞うかを検討するという意味では不適切な実装であるが、
* 重症者集団が一種の遺伝子バンクとして機能してしまう、医療崩壊した地域ではありうる挙動ではある
*
* ========引数========
*
* ========詳細========
*
*/
int daily_mutation2(DAILYCALENDAR *dailycalendar){
int i, mut, virs;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
// 突然変異処理は、「daily_infection()感染」の直前に行うので、yesterdayに対して行う
flock = dailycalendar->yesterday;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
i = 0;
while( find_infection( flock, &i ) ){//感染者を探す
if( find_Virus_type( &(flock->humans[i]), VIRUS_2 ) && find_Virus_type( &(flock->humans[i]), VIRUS_3 ) ){
//_tprintf( _T("Warning!!!!!!!!!!!!W\n") );
//_tprintf( _T("%d日目にタイプ2とタイプ3の変異型ウイルスに同時感染し危険な変異型4発生\n"), dailycalendar->day );
virus_injection( &(dailycalendar->today->humans[ i ]), VIRUS_4, 1);
}
mut = irand( MUTATION_PROBABILITY );
if( mut==MUTATION_PROBABILITY-1 ){//変異率の逆数の乱数(0～MUTATION_PROBABILITY-1)を取得し、変異判定
virs = irand( 2 );//変異率に対して、変異型1 OR 2のどちらになるかはさらに二択。
//_tprintf( _T("%d日目に変異型ウイルス%d発生\n"), dailycalendar->day, virs+2 );
virus_injection( &(dailycalendar->today->humans[ i ]), virs+2, 1);
}
i++;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分のウイルス増殖を一通り起こす
* この関数は、各ウイルス個体数合計がVIRUS_MAXまで増殖する。
* つまり、人間の体内ではVIRUS_MAXまでしかウイルスが増えることはできない。
* ウイルスの増殖に使える資源は当然一人の人間の持つ体という限界点がある。
* 実稼動では本関数を使うこと。
* ========引数========
*
* ========詳細========
*
*/
int daily_proliferation(DAILYCALENDAR *dailycalendar){
int i, x;
long long tmp, sum1, sum2;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
flock = dailycalendar->today;
for(i=0;i<flock->n;i++){
if( flock->humans[i].virusn>0 ){//感染発見
sum1 = SUM_POPULATION( flock->humans[i] );
if( sum1 >= VIRUS_MAX ){
if( flock->humans[i].count<50 ){//継続感染はカウントする意味がない
(flock->humans[i].count)++;
}
continue;
}
for(x=0;x<VIRUS_N;x++){
if( flock->humans[i].sars_v[x].type!=VIRUS_0 ){
tmp = flock->humans[i].sars_v[x].Population;
//_tprintf( _T("増殖予定ウイルス%d %lld → %lld\n"), x+1, flock->humans[i].sars_v[x].Population, (flock->humans[i].sars_v[x].Population) * PROLIFERATION );
flock->humans[i].sars_v[x].Population = (flock->humans[i].sars_v[x].Population) * PROLIFERATION;
sum2 = SUM_POPULATION( flock->humans[i] );
if( sum2 > VIRUS_MAX ){
//_tprintf( _T("sum=%lld tmp= %lld VIRUS_MAX=%lld\n"), sum2, tmp, VIRUS_MAX );
//_tprintf( _T("限界値なので修正%lld → %lld\n"), flock->humans[i].sars_v[x].Population, tmp + VIRUS_MAX - sum1 );
flock->humans[i].sars_v[x].Population = tmp + VIRUS_MAX - sum1;
}
//else if( flock->humans[i].sars_v[x].Population < tmp ){ //事前にバックアップしておいた値より減っていたら限界とする
// flock->humans[i].sars_v[x].Population = tmp;
// //_tprintf( _T("限界値なので修正%lld → %lld\n"), flock->humans[i].sars_v[x].Population, tmp );
//}
}
}
(flock->humans[i].count)++;
//_tprintf( _T("(flock->humans[i].count) %d\n"), (flock->humans[i].count) );
//_tprintf( _T("v1=%lld\n"), flock->humans[i].sars_v[0].Population );
//_tprintf( _T("v2=%lld\n"), flock->humans[i].sars_v[1].Population );
//_tprintf( _T("v3=%lld\n"), flock->humans[i].sars_v[2].Population );
//_tprintf( _T("v4=%lld\n\n"), flock->humans[i].sars_v[3].Population );
}
//else if( flock->humans[i].condition == CONDITION_INCUB){
// _tprintf( _T("flock->humans[i].virusn=%d\n\n"), flock->humans[i].virusn );
// _tprintf( _T("flock->humans[i].condition=%d\n\n"), flock->humans[i].condition );
//}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 1日分のウイルス増殖を一通り起こす
* この関数は単純で、各ウイルス個体数がそれぞれVIRUS_MAXまで増殖する。
* つまり、4種類のウイルスに感染した人間の体内では4×VIRUS_MAXまで増える。
* これは、一人の人間の体内に存在できる最大数がVIRUS_MAXであるという設計理念と反する。
* 実稼動ではdaily_proliferation()を使うこと。
* ========引数========
*
* ========詳細========
*
*/
int daily_proliferation2(DAILYCALENDAR *dailycalendar){
int i, x;
long long tmp;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
flock = dailycalendar->today;
for(i=0;i<flock->n;i++){
if( flock->humans[i].virusn>0 ){//感染発見
for(x=0;x<VIRUS_N;x++){
if( flock->humans[i].sars_v[x].type!=VIRUS_0 ){
tmp = flock->humans[i].sars_v[x].Population;
//_tprintf( _T("増殖予定%lld → %lld\n"), flock->humans[i].sars_v[x].Population, (flock->humans[i].sars_v[x].Population) * PROLIFERATION );
flock->humans[i].sars_v[x].Population = (flock->humans[i].sars_v[x].Population) * PROLIFERATION;
if( flock->humans[i].sars_v[x].Population > VIRUS_MAX ){
//_tprintf( _T("限界値なので修正%lld → %lld\n"), flock->humans[i].sars_v[x].Population, VIRUS_MAX );
flock->humans[i].sars_v[x].Population = VIRUS_MAX;
}
else if( flock->humans[i].sars_v[x].Population < tmp ){ //事前にバックアップしておいた値より減っていたら限界とする
flock->humans[i].sars_v[x].Population = tmp;
//_tprintf( _T("限界値なので修正%lld → %lld\n"), flock->humans[i].sars_v[x].Population, tmp );
}
}
}
(flock->humans[i].count)++;
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* todayのイベント発生条件を調べ、該当する場合に当該イベントを発生させる
* さらに、緊急事態宣言の発令と解除で違う値を返す
*
* ========引数========
*
* ========詳細========
*
* チェック条件例
* ・緊急事態宣言State of Emergency の条件チェック
* ・emergency_flg と 現在の感染状況の比較チェック
*
*
*
* イベント例
*
*
*
*
*/
int daily_event( DAILYCALENDAR *dailycalendar ){
int ret = 1;
int patient;
if( dailycalendar==NULL ) return 0;
if( daily_event_sub_PatientCount( dailycalendar, &patient)==0 ) return -1;
if( daily_event_sub_EmergencyOnCheck(dailycalendar, patient)==1 &&
dailycalendar->emergency_flg != EMERGENCY_FLG_ON ){//緊急事態宣言に値するのにまだ宣言されていない
dailycalendar->emergency_flg = EMERGENCY_FLG_ON;
_tprintf( _T("緊急事態宣言 発令\n") );
ret = 2;
}
else if( daily_event_sub_EmergencyOffCheck(dailycalendar, patient)==1 &&
dailycalendar->emergency_flg != EMERGENCY_FLG_OFF ){//緊急事態宣言解除に値するのにまだ解除されていない
dailycalendar->emergency_flg = EMERGENCY_FLG_OFF;
_tprintf( _T("緊急事態宣言 解除\n") );
ret = 3;
}
if( daily_event_sub_MedicalCapacityOverCheck(dailycalendar, patient)==1 ){
dailycalendar->medicalover_flg = MEDICALCAPACITYOVER_FLG_ON;
//_tprintf( _T("医療崩壊しました\n") );
}
return ret;
}
/***************************************************************
*
* ========目的========
* daily_event()の下請け。緊急事態宣言発令に値するか調べ、
* 発令するべき感染者数なら1を、発令すべきではなければ0を返す
*
* ========引数========
*
* ========詳細========
*
*/
int daily_event_sub_EmergencyOnCheck(DAILYCALENDAR *dailycalendar, int patient){
int patient_limit;
if( dailycalendar==NULL ) return -1;
patient_limit = (dailycalendar->statistics.n - dailycalendar->statistics.dead) * (dailycalendar->actions.self_restraint_limit_m)/PER_DEF;
if( patient >= patient_limit ) return 1;
return 0;
}
/***************************************************************
*
* ========目的========
* daily_event()の下請け。緊急事態宣言解除に値するか調べ、
* 解除するべき感染者数なら1を、解除すべきではなければ0を返す
*
* ========引数========
*
* ========詳細========
*
*/
int daily_event_sub_EmergencyOffCheck(DAILYCALENDAR *dailycalendar, int patient){
int patient_limit;
if( dailycalendar==NULL ) return -1;
patient_limit = (dailycalendar->statistics.n - dailycalendar->statistics.dead) * ((dailycalendar->actions.self_restraint_limit_m)/2)/PER_DEF;
if( patient < patient_limit ) return 1;
return 0;
}
/***************************************************************
*
* ========目的========
* daily_event()の下請け。緊急事態宣言に関わる、現在の患者数をカウントしpatientに格納する
* 関数が無事終了すれば1を返す。エラーは0を返す
*
* ========引数========
*
* ========詳細========
*
*/
int daily_event_sub_PatientCount(DAILYCALENDAR *dailycalendar, int *patient){
int i;
FLOCK *flock;
if( dailycalendar==NULL ) return 0;
(*patient) = 0;
flock = dailycalendar->today;//何度もアクセスするポインタは次数を減らした方が高速化する可能性がある
for(i=0;i<flock->n;i++){
if( flock->humans[i].condition==CONDITION_ALIVE ){//健康に動いて生きている
if( flock->humans[i].virusn>0 ){//未発症感染者
(*patient)++;
}
}
else if( flock->humans[i].condition==CONDITION_INCUB ){//発症している
(*patient)++;
}
else if( flock->humans[i].condition==CONDITION_SEVEE ){//重症
(*patient)++;
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* daily_event()の下請け。医療のキャパシティを超えていないかチェックする。
* ここでいうキャパシティーオーバーとは、俗に言う医療崩壊を指す。
* キャパシティーを超えていたら1を、キャパシティーを超えてなければ0を返す
*
* ========引数========
*
* ========詳細========
*
*/
int daily_event_sub_MedicalCapacityOverCheck(DAILYCALENDAR *dailycalendar, int patient){
if( dailycalendar==NULL ) return -1;
if( (dailycalendar->today->n)*(MEDICALCAPACITYOVER)/PER_DEF < patient ) return 1;
return 0;
}
/***************************************************************
*
* ========目的========
* 実際のウイルス注入(injection)処理。指定された人間に、指定タイプのウイルスを指定個数侵入させる最も基本的な関数
* 体質等の要素は一切考慮しない、処理としてはウイルス注入に近い。
* 突然変異による変異型発生も、実際の処理としては本関数のように注入することである→ただしRtの関係で派生関数があるvirus_injection2
*
* ========引数========
*
* ========詳細========
*
*/
int virus_injection(HUMANS *humans, int virus_type, int population){
int x;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){//まずは既に感染していないかチェック
if( humans->sars_v[x].type==virus_type ){//既に感染しているので個体数だけ増加
humans->sars_v[x].Population = humans->sars_v[x].Population + population;
//_tprintf( _T("既に感染済み\n") );
return 1;
}
}
for(x=0;x<VIRUS_N;x++){//空き領域を探す
//_tprintf( _T("x=%d humans->sars_v[x].type=%d\n"), x, humans->sars_v[x].type );
if( humans->sars_v[x].type==VIRUS_0 ){//空き発見
humans->sars_v[x].type = virus_type;
humans->sars_v[x].Population = population;
if ( virus_type==VIRUS_1 ){humans->sars_v[x].Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_2 ){humans->sars_v[x].Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_3 ){humans->sars_v[x].Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_4 ){humans->sars_v[x].Rt = REPTNUM_VIRUS_4;}
//_tprintf( _T("humans->sars_v[%d].type=%d\n"), x, humans->sars_v[x].type );
//_tprintf( _T("humans->sars_v[%d].Population=%d\n"), x, humans->sars_v[x].Population );
//_tprintf( _T("humans->sars_v[%d].Rt=%f\n"), x, humans->sars_v[x].Rt );
break;
}
}
(humans->virusn)++;
if( humans->count<0 ) humans->count = 0;
return 1;
}
/***************************************************************
*
* ========目的========
* 実際のウイルス注入(injection)処理。
* 変異によるウイルス発生は本関数で再現する。
* virus_injection()との違いは、実効再生産数Rtを既にいるウイルスと同じにする点である。
* 突然変異による変異型発生も、実際の処理としては本関数で注入することである
*
* ========引数========
*
* ========詳細========
* 本シミュレーションでは、ウイルス1からウイルス2やウイルス3のタイプが変異し、
* 同一の人間の体内にウイルス2とウイルス3が両方存在したとき、2,3の変異遺伝子(RNA)を
* 取り込んだハイブリッドとしてウイルス4が生まれる。
* その変異ウイルスだが、基本的には1に感染後ある程度時間がたってから変異することになるはずである。
* しかし、virus_injection()では実効再生産数Rtをデフォルト値設定する。これは実態に合わない。
* なぜなら、実効再生産数Rtは一人の人間に対して算出される値であり、本プログラムでは感染機会チケットの扱いである。
* 例えば、タイプ1がRt=2.4で0.9まで使い切ったところで変異が起こってタイプ2が生じた場合、
* if( virus_type==VIRUS_2 ){humans->sars_v[x].Rt = REPTNUM_VIRUS_1;}のようにRtをデフォルト値で初期化してはならない。
* たとえウイルスの株の種類が増えても、感染機会まで増えてはならないからである。(REPTNUM_VIRUS=2.4)
* もしも許せば、この人間はタイプ1の0.9枚のチケットを使い切ったあとにさらに差額の1.5回分タイプ2の感染機会を持つことになる。
* 実効再生産数は一人の人間がどれだけ新規の患者を獲得するかの大まかな指標なので、
* 変異するたびに増えるのはおかしいのである。最悪の場合、チケットが切れる前に新たな変異が起こればさらに
* 感染機会は増えるのだ。これを防止するために、本関数ではRtを既にあるウイルスからコピーして帳尻を合わせる。
*
*/
int virus_injection2(HUMANS *humans, int virus_type, int population){
int x, Rt;
if( humans==NULL ) return 0;
for(x=0;x<VIRUS_N;x++){//まずは既に感染していないかチェック
if( humans->sars_v[x].type==virus_type ){//既に感染しているので個体数だけ増加
humans->sars_v[x].Population = humans->sars_v[x].Population + population;
//_tprintf( _T("既に感染済み\n") );
return 1;
}
}
Rt = 0;
if ( virus_type==VIRUS_1 ){ Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_2 ){ Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_3 ){ Rt = REPTNUM_VIRUS_1;}
else if( virus_type==VIRUS_4 ){ Rt = REPTNUM_VIRUS_4;}
for(x=0;x<VIRUS_N;x++){//Rtを探す
if( humans->sars_v[x].type!=VIRUS_0 ){//空き発見
Rt = humans->sars_v[x].Rt;
break;
}
}
for(x=0;x<VIRUS_N;x++){//空き領域を探す
//_tprintf( _T("x=%d humans->sars_v[x].type=%d\n"), x, humans->sars_v[x].type );
if( humans->sars_v[x].type==VIRUS_0 ){//空き発見
humans->sars_v[x].type = virus_type;
humans->sars_v[x].Population = population;
humans->sars_v[x].Rt = Rt;
//_tprintf( _T("humans->sars_v[x].Rt=%f\n"), humans->sars_v[x].Rt );
break;
}
}
(humans->virusn)++;
if( humans->count<0 ) humans->count = 0;
return 1;
}
/***************************************************************
*
* ========目的========
* 実際のウイルス侵入処理。指定された人間に、指定タイプのウイルスを指定個数侵入させる
* Population オーダーは×100。つまり、Populationが1なら100個体として扱い、感染成立。
*
* ========引数========
*
* ========詳細========
*
*/
int virus_invasion(HUMANS *humans, int virus_type, int population){
if( humans==NULL ) return 0;
if( humans->immunity == IMMUNITY_ON ){// 免疫があるので感染不成立
//_tprintf( _T("免疫があるので感染しませんでした\n") );
return 1;
}
if( humans->mask == MASK_ON ){
population = population * MASK_DEFENSE;
if( population<1 ) return 1;
}
virus_injection(humans, virus_type, population);
if( humans->cons == CONS_HIDE ){virus_invasion_sub_cons_hide( humans );}
return 1;
}
//virus_invasion()の下請け。体質が不顕性感染タイプの場合、Rtを半分に減ずる(感染後半にもう一度感染をおこすため)
int virus_invasion_sub_cons_hide(HUMANS *humans){
int x;
if(humans==NULL) return 0;
for(x=0;x<VIRUS_N;x++){
if( humans->sars_v[x].Rt > 0 ) humans->sars_v[x].Rt = (humans->sars_v[x].Rt)/2;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 当日の感染状況の統計データを計算して格納する
*
* ========引数========
*
* ========詳細========
* 本来は、計算した日付を保持することで無駄な計算を省けるはずだが、他の出力関数がいろいろなタイミングで
* 呼び出す関係上、最新データであることを保証するのが難しいため常に値を計算する仕様。
*
*/
int daily_statistics_calc(DAILYCALENDAR *dailycalendar){
if( dailycalendar==NULL ) return 0;
//if( dailycalendar->day != dailycalendar->statistics.day ){//カレンダーの日付と統計データの日付が違う
get_statistics_flock(dailycalendar, &(dailycalendar->statistics) );
//}
return 1;
}
// 感染を実行する関数。
// 感染元のHUMAN構造体のメンバ構造体SARS_Vの実行再生産数を元に感染を試行する。
// 最も単純な実装では、1回の感染で実行再生産数を1減ずる。例えば、実行再生産数が2.4ならば、
// 1回目 感染で2.4 → 1.4
// 2回目 感染で1.4 → 0.4
// 3回目 感染で0.4 → 0.0 確率は乱数で1～10まで生成し、4以下で感染成立、5以上で感染不成立
// となる。
// しかし、実際にはある程度のランダム性があるべきなので、例えば自由度を上げて
// 次のように実装すると、実行再生産数を確率的に反映した上でランダム性が増す。
// 自由度50%として1回の感染で実行再生産数を0.5減ずる。例えば、実行再生産数が2.4ならば、
// 1回目 感染で2.4 → 1.9 確率は乱数で1～10まで生成し、5以下で感染成立、6以上で感染不成立
// 1回目 感染で1.9 → 1.4 確率は乱数で1～10まで生成し、5以下で感染成立、6以上で感染不成立
// 1回目 感染で1.4 → 0.9 確率は乱数で1～10まで生成し、5以下で感染成立、6以上で感染不成立
// 1回目 感染で0.9 → 0.4 確率は乱数で1～10まで生成し、5以下で感染成立、6以上で感染不成立
// 1回目 感染で0.4 → 0.0 確率は乱数で1～10まで生成し、4以下で感染成立、5以上で感染不成立
// となる。
/***************************************************************
*
* ========目的========
*
* ========引数========
*
* ========詳細========
*
*/
int infection_dice(HUMANS *humans ){
if( humans==NULL ) return 0;
return 1;
}
/***************************************************************
*
* ========目的========
* 感染をシミュレートする。
*
*
* ========引数========
* int seed;//rand関数の初期化に使った値。後でシミュレーションを再現するときに使う用
*
* ========詳細========
* 乱数生成器を初期化する関数srand( seed );の種となった値seedは記録しておく(厳密には擬似乱数)。
* もしも、シミュレーションに奇妙な動きがあった場合、再現するためには乱数生成パターンも再現する必要がある。
* Cのrand()で一般的に用いられる乱数生成アルゴリズムは線形合同法といわれており、初期値はsrand()を通して設定できる。
* このsrand()に与える初期値があれば、同じ擬似乱数列を得られるため、全く同じ感染パターンとなる。
*
*
*/
int Simulation_Run( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line ){
#ifdef Details
FILE *fp1, *fp3;
#endif
FILE *fp2;
int day, ret;
DAILYCALENDAR *dailycalendar;
if( fp_log==NULL || fp_glp==NULL ) return 0;
//fp1 = now_date_fopen( _T("感染ログ"), _T(".log") );
//fp2 = now_date_fopen( _T("統計情報"), _T(".log") );
fp2 = now_date_fopen_n( _T("統計情報"), _T(".log") );
#ifdef Details
fp1 = now_date_fopen_n( _T("感染ログ"), _T(".log") );
fp3 = now_date_fopen_n( _T("度数分布"), _T(".log") );
#endif
if( n < 100 ) n = 1000;
if( nday< 1 ) nday = 20;
//srand( seed );
dailycalendar = DailyCalendar_make( n, actions );
fprintf_daily_Calendar_setting(fp2, world_line, dailycalendar);
flock_property_print(dailycalendar->today, fp2);
virus_injection( &(dailycalendar->today->humans[1]), VIRUS_1, 1);//感染源注入
virus_injection( &(dailycalendar->today->humans[2]), VIRUS_1, 1);//感染源注入
// virus_injection( &(dailycalendar->today->humans[1]), VIRUS_2, 1);//感染源注入
#ifdef Details
_ftprintf(fp2, _T("乱数生成の種は%d\r\n"), seed );
//fprintf_dailycalendar(fp1, dailycalendar);//感染者の詳細まで全て記録する。めちゃ重い上データもでかい
_ftprintf(fp3, _T("人口%d人\r\n"), dailycalendar->today->n );
#endif
for(day=0;day<nday;day++){
dayAdvance( dailycalendar );//カレンダーの日めくり(yesterdayをtodayへコピー)
daily_mutation( dailycalendar );//変異(yesterdayから算出したウイルス量をtodayへ書き込み)
daily_proliferation( dailycalendar );//増殖(todayのウイルス量へPROLIFERATIONを掛け合わせる)
daily_condition( dailycalendar );//治癒や重症化など体調変化(todayの内容を参照し、必要に応じてtodayのパラメータを書き換え)
daily_infection2( dailycalendar );//感染(yesterdayから感染元となりうる条件を満たした人間をみつけ、todayへ感染させる)
ret = daily_event( dailycalendar );
if ( ret==2 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言発令\r\n"), day+1 );
else if( ret==3 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言解除\r\n"), day+1 );
fprintf_daily_graph(fp_glp, dailycalendar, world_line);
//_ftprintf(fp_log, _T("第%2d世界\t"), world_line );
//fprintf_daily_statistics3(fp_log, dailycalendar);
fprintf_daily_statistics(fp2, dailycalendar);//todayの統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
fprintf_humans_statistics(fp2, dailycalendar);//todayの感染者各個人の統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
fflush( fp2 );
#ifdef Details
//fprintf_dailycalendar(fp1, dailycalendar);//感染者の詳細まで全て記録する。めちゃ重い上データもでかい
//fprintf_daily_virus_Frequency_distribution(fp3, dailycalendar);//ウイルス毎の感染者数の度数分布表出力
fprintf_daily_virus_Frequency_distribution_Sorting(fp3, dailycalendar);
fflush( fp3 );
#endif
}
daily_statistics_calc( dailycalendar );
//if( dailycalendar.statistics.sars_v_total[VIRUS_4-1] > 0 ){
_ftprintf(fp_log, _T("ウイルスタイプ4\t%d\r\n"), dailycalendar->statistics.sars_v_total[VIRUS_4-1] );
//}
//fprintf_daily_statistics2(fp_log, dailycalendar);//世界線記録用
//fprintf_daily_statistics3(fp_log, dailycalendar);//世界線記録用
fflush( fp_log );
DailyCalendar_free( dailycalendar ) ;
fclose(fp2);
#ifdef Details
fclose(fp1);
fclose(fp3);
#endif
return 1;
}
/***************************************************************
*
* ========目的========
* 感染をシミュレートする。
*
*
* ========引数========
* int seed;//rand関数の初期化に使った値。後でシミュレーションを再現するときに使う用
*
* ========詳細========
*
*
*
*/
int Simulation_Run2( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line ){
#ifdef Details
FILE *fp1, *fp3;
#endif
FILE *fp2;
int day;
DAILYCALENDAR *dailycalendar;
if( fp_log==NULL || fp_glp==NULL ) return 0;
fp2 = now_date_fopen( _T("統計情報"), _T(".log") );
#ifdef Details
fp1 = now_date_fopen_n( _T("感染ログ"), _T(".log") );
fp3 = now_date_fopen_n( _T("度数分布"), _T(".log") );
#endif
if( n < 100 ) n = 1000;
if( nday< 1 ) nday = 20;
srand( seed );
dailycalendar = DailyCalendar_make( n, actions );
fprintf_daily_Calendar_setting(fp2, world_line, dailycalendar);
flock_property_print(dailycalendar->today, fp2);
virus_injection( &(dailycalendar->today->humans[1]), VIRUS_1, 1);//感染源注入
virus_injection( &(dailycalendar->today->humans[2]), VIRUS_1, 1);//感染源注入
#ifdef Details
_ftprintf(fp2, _T("乱数生成の種は%d\r\n"), seed );
fprintf_dailycalendar(fp1, dailycalendar);
_ftprintf(fp3, _T("人口%d人\r\n"), dailycalendar->today->n );
#endif
for(day=0;day<nday;day++){
dayAdvance( dailycalendar );//カレンダーの日めくり
daily_mutation( dailycalendar );//変異
daily_proliferation( dailycalendar );//増殖
daily_condition( dailycalendar );//治癒や重症化など体調変化
daily_infection2( dailycalendar );//感染
daily_event( dailycalendar );
#ifdef Details
//fprintf_dailycalendar(fp1, dailycalendar);
fprintf_daily_statistics(fp2, dailycalendar);//todayの統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
fprintf_humans_statistics(fp2, dailycalendar);//todayの感染者各個人の統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
//fprintf_daily_virus_Frequency_distribution(fp3, dailycalendar);//ウイルス毎の感染者数の度数分布表出力
fprintf_daily_virus_Frequency_distribution_Sorting(fp3, dailycalendar);
fflush( fp2 );
fflush( fp3 );
#endif
}
daily_statistics_calc( dailycalendar );
_ftprintf(fp_log, _T("ウイルスタイプ4\t%d\r\n"), dailycalendar->statistics.sars_v_total[VIRUS_4-1] );
//fprintf_daily_statistics2(fp_log, dailycalendar);//世界線記録用
fflush( fp_log );
DailyCalendar_free( dailycalendar ) ;
fclose(fp2);
#ifdef Details
fclose(fp1);
fclose(fp3);
#endif
return 1;
}
/***************************************************************
*
* ========目的========
* 感染をシミュレートする。
* とにかく軽量で、出力等も最小限
*
* ========引数========
* int seed;//rand関数の初期化に使った値。後でシミュレーションを再現するときに使う用
*
* ========詳細========
*
*
*
*/
int Simulation_Run3( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line ){
#ifdef Details
FILE *fp1, *fp3;
FILE *fp2;
#endif
int day;
DAILYCALENDAR *dailycalendar;
if( fp_log==NULL || fp_glp==NULL ) return 0;
#ifdef Details
fp1 = now_date_fopen_n( _T("感染ログ"), _T(".log") );
fp2 = now_date_fopen( _T("統計情報"), _T(".log") );
fp3 = now_date_fopen_n( _T("度数分布"), _T(".log") );
#endif
if( n < 100 ) n = 1000;
if( nday< 1 ) nday = 20;
//srand( seed );
dailycalendar = DailyCalendar_make( n, actions );
flock_property_print(dailycalendar->today, fp_glp);
actions_property_print(dailycalendar, fp_glp);
//_ftprintf(fp_glp, _T("seed=%d\r\n"), seed );
_ftprintf(fp_glp, _T("パラメータ\t経過日数\t世界線\t数\r\n") );
virus_injection( &(dailycalendar->today->humans[1]), VIRUS_1, 1);//感染源注入
virus_injection( &(dailycalendar->today->humans[2]), VIRUS_1, 1);//感染源注入
#ifdef Details
fprintf_daily_Calendar_setting(fp2, world_line, dailycalendar);
flock_property_print(dailycalendar->today, fp2);
_ftprintf(fp2, _T("乱数生成の種は%d\r\n"), seed );
fprintf_dailycalendar(fp1, dailycalendar);
_ftprintf(fp3, _T("人口%d人\r\n"), dailycalendar->today->n );
#endif
for(day=0;day<nday;day++){
dayAdvance( dailycalendar );//カレンダーの日めくり
daily_mutation( dailycalendar );//変異
daily_proliferation( dailycalendar );//増殖
daily_condition( dailycalendar );//治癒や重症化など体調変化
daily_infection2( dailycalendar );//感染
daily_statistics_calc( dailycalendar );
daily_event( dailycalendar );
//fprintf_daily_graph(fp_glp, dailycalendar, world_line);
#ifdef Details
//fprintf_dailycalendar(fp1, dailycalendar);
fprintf_daily_statistics(fp2, dailycalendar);//todayの統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
fprintf_humans_statistics(fp2, dailycalendar);//todayの感染者各個人の統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
//fprintf_daily_virus_Frequency_distribution(fp3, dailycalendar);//ウイルス毎の感染者数の度数分布表出力
fprintf_daily_virus_Frequency_distribution_Sorting(fp3, dailycalendar);
fflush( fp2 );
fflush( fp3 );
#endif
}
daily_statistics_calc( dailycalendar );
//_ftprintf(fp_log, _T("ウイルスタイプ4\t%d\r\n"), dailycalendar->statistics.sars_v_total[VIRUS_4-1] );
_ftprintf(fp_log, _T("%d\t"), dailycalendar->statistics.sars_v_total[VIRUS_4-1] );//ウイルス種4の感染者数
//_ftprintf(fp_log, _T("%d\t"), dailycalendar->statistics.sars_v_total[VIRUS_1-1] );//ウイルス種1の感染者数
//fprintf_daily_statistics2(fp_log, dailycalendar);//世界線記録用
fflush( fp_log );
DailyCalendar_free( dailycalendar ) ;
#ifdef Details
fclose(fp1);
fclose(fp2);
fclose(fp3);
#endif
return 1;
}
/***************************************************************
*
* ========目的========
* 感染をシミュレートする。
* 度数分布は出力するが、統計情報は出力しない。統計情報ファイルには、実行時のパラメータだけ記録される
* 本関数は、主に特定のパラメータのときの感染動向などを外部に見せるための情報を出力する
*
* ========引数========
* int seed;//rand関数の初期化に使った値。後でシミュレーションを再現するときに使う用
*
* ========詳細========
*
*
*
*/
int Simulation_Run4( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line ){
#ifdef Details
FILE *fp1;
#endif
FILE *fp3;
FILE *fp2;
int day, re;
DAILYCALENDAR *dailycalendar;
TCHAR fname[ FNAME_MAX+1 ];
if( fp_log==NULL || fp_glp==NULL ) return 0;
_sntprintf(fname, FNAME_MAX, _T("統計情報_世界%03d線%05d人%03dクラスタ"), world_line, n, actions->modulo_M );
//fp2 = now_date_fopen_n( _T("統計情報"), _T(".log") );
fp2 = now_date_fopen_n( fname, _T(".log") );
#ifdef Details
fp1 = now_date_fopen_n( _T("感染ログ"), _T(".log") );
#endif
_sntprintf(fname, FNAME_MAX, _T("度数分布_世界%03d線%05d人%03dクラスタ"), world_line, n, actions->modulo_M );
//fp3 = now_date_fopen_n( _T("度数分布"), _T(".log") );
fp3 = now_date_fopen_n( fname, _T(".log") );
dailycalendar = DailyCalendar_make( n, actions );
fprintf_daily_Calendar_setting(fp2, world_line, dailycalendar);
flock_property_print(dailycalendar->today, fp2);
virus_injection( &(dailycalendar->today->humans[1]), VIRUS_1, 1);//感染源注入
virus_injection( &(dailycalendar->today->humans[2]), VIRUS_1, 1);//感染源注入
_ftprintf(fp2, _T("乱数生成の種は%d\r\n"), seed );
#ifdef Details
fprintf_dailycalendar(fp1, dailycalendar);
#endif
_ftprintf(fp3, _T("人口%d人\r\n"), dailycalendar->today->n );
for(day=0;day<nday;day++){
dayAdvance( dailycalendar );//カレンダーの日めくり
if( day==50) virus_injection( &(dailycalendar->today->humans[99]), VIRUS_4, 100);//感染源注入
daily_mutation( dailycalendar );//変異
daily_proliferation( dailycalendar );//増殖
daily_condition( dailycalendar );//治癒や重症化など体調変化
daily_infection2( dailycalendar );//感染
daily_statistics_calc( dailycalendar );
re = daily_event( dailycalendar );
if ( re==2 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言発令\r\n"), day+1 );
else if( re==3 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言解除\r\n"), day+1 );
fprintf_daily_graph(fp_glp, dailycalendar, world_line);
fprintf_daily_virus_Frequency_distribution_Sorting(fp3, dailycalendar);
fflush( fp3 );fflush( fp2 );fflush( fp_glp );
#ifdef Details
//fprintf_dailycalendar(fp1, dailycalendar);
fprintf_daily_statistics(fp2, dailycalendar);//todayの統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
fprintf_humans_statistics(fp2, dailycalendar);//todayの感染者各個人の統計情報を記録する
_ftprintf(fp2, _T("\r\n") );
//fprintf_daily_virus_Frequency_distribution(fp3, dailycalendar);//ウイルス毎の感染者数の度数分布表出力
fflush( fp2 );
#endif
}
daily_statistics_calc( dailycalendar );
_ftprintf(fp_log, _T("ウイルスタイプ4\t%d\r\n"), dailycalendar->statistics.sars_v_total[VIRUS_4-1] );
//fprintf_daily_statistics2(fp_log, dailycalendar);//世界線記録用
fflush( fp_log );
DailyCalendar_free( dailycalendar ) ;
fclose(fp2);
fclose(fp3);
#ifdef Details
fclose(fp1);
#endif
return 1;
}
/***************************************************************
*
* ========目的========
* 感染をシミュレートする。
* 医療崩壊の条件を探るため、医療崩壊の回数をカウントする
* 度数分布は出力するが、統計情報は出力しない。統計情報ファイルには、実行時のパラメータだけ記録される
* 本関数は、主に特定のパラメータのときの感染動向などを外部に見せるための情報を出力する
*
* ========引数========
* int seed;//rand関数の初期化に使った値。後でシミュレーションを再現するときに使う用
*
* ========詳細========
*
*
*
*/
int Simulation_Run5( int n, int nday, FILE *fp_log, FILE *fp_glp, int seed, ACTIONS *actions, int world_line ){
FILE *fp3;
FILE *fp2;
int day, re, ret;
DAILYCALENDAR *dailycalendar;
TCHAR fname[ FNAME_MAX+1 ];
if( fp_log==NULL || fp_glp==NULL ) return 0;
_sntprintf(fname, FNAME_MAX, _T("統計情報_世界%03d線%05d人%03dクラスタ"), world_line, n, actions->modulo_M );
//fp2 = now_date_fopen_n( _T("統計情報"), _T(".log") );
fp2 = now_date_fopen_n( fname, _T(".log") );
_sntprintf(fname, FNAME_MAX, _T("度数分布_世界%03d線%05d人%03dクラスタ"), world_line, n, actions->modulo_M );
//fp3 = now_date_fopen_n( _T("度数分布"), _T(".log") );
fp3 = now_date_fopen_n( fname, _T(".log") );
dailycalendar = DailyCalendar_make( n, actions );
fprintf_daily_Calendar_setting(fp2, world_line, dailycalendar);
flock_property_print(dailycalendar->today, fp2);
virus_injection( &(dailycalendar->today->humans[1]), VIRUS_1, POPULATION_ORDER);//感染源注入
virus_injection( &(dailycalendar->today->humans[2]), VIRUS_1, POPULATION_ORDER);//感染源注入
virus_injection( &(dailycalendar->today->humans[3]), VIRUS_1, POPULATION_ORDER);//感染源注入
virus_injection( &(dailycalendar->today->humans[4]), VIRUS_1, POPULATION_ORDER);//感染源注入
_ftprintf(fp2, _T("乱数生成の種は%d\r\n"), seed );
_ftprintf(fp3, _T("人口%d人\r\n"), dailycalendar->today->n );//度数分布出力
for(day=0;day<nday;day++){
ret = dayAdvance( dailycalendar );//カレンダーの日めくり
if( day==50){
virus_injection( &(dailycalendar->today->humans[100]), VIRUS_4, 100);//感染源注入
virus_injection( &(dailycalendar->today->humans[101]), VIRUS_4, 100);//感染源注入
virus_injection( &(dailycalendar->today->humans[102]), VIRUS_4, 100);//感染源注入
virus_injection( &(dailycalendar->today->humans[103]), VIRUS_4, 100);//感染源注入
}
daily_mutation( dailycalendar );//変異
daily_proliferation( dailycalendar );//増殖
daily_condition( dailycalendar );//治癒や重症化など体調変化
daily_infection2( dailycalendar );//感染
daily_statistics_calc( dailycalendar );
re = daily_event( dailycalendar );
if ( re==2 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言発令\r\n"), day+1 );
else if( re==3 ) _ftprintf(fp2, _T("%04d日から 緊急事態宣言解除\r\n"), day+1 );
//fprintf_daily_graph(fp_glp, dailycalendar, world_line);
fprintf_daily_virus_Frequency_distribution_Sorting(fp3, dailycalendar);//度数分布出力
if( ret!=1 && day>50 ){//これ以上続けても変化しない恒常状態になったので終わり
break;
}
fflush( fp3 );fflush( fp2 );//fflush( fp_glp );
}
daily_statistics_calc( dailycalendar );
if(dailycalendar->medicalover_flg==MEDICALCAPACITYOVER_FLG_ON){//医療崩壊が起きた
_tprintf( _T("医療崩壊しました\n") );
_ftprintf(fp_log, _T("1\t") );
}
else{
_ftprintf(fp_log, _T("0\t") );
}
//fprintf_daily_statistics2(fp_log, dailycalendar);//世界線記録用
if(dailycalendar->medicalover_flg==MEDICALCAPACITYOVER_FLG_ON){//医療崩壊が起きた
_ftprintf(fp2, _T("医療崩壊した\r\n") );
}
fflush( fp_log );fflush( fp_glp );
DailyCalendar_free( dailycalendar ) ;
fclose(fp2);
fclose(fp3);
return 1;
}
/***************************************************************
*
* ========目的========
* actions_property_set() の設定のサンプル。main関数側で直接指定しても良いが、
* 少し長くなるのと、パラメーターが多いので指定方法のサンプルとしていくつか定義する
* 非常に無防備な集団の設定。
*
*/
int actions_property_set_ex1( ACTIONS *actions ){
int ret;
if( actions==NULL ) return 0;
ret = actions_property_set( actions,
3000, //int mask,マスクをつける確率
0, //int cons_strong,無敵個体の、全人口に対する割合
0, //int cons_hide,隠れ感染者の、全人口に対する割合
9800, //int cons_middle,通常感染者の、全人口に対する割合
100, //int cons_severe,重症化する人間の、全人口に対する割合
100, //int cons_death,感染後14日目に死ぬ人間の、全人口に対する割合
ACT_SELF_RESTRAINT_M, //int self_restraint,行動を自粛する割合。
ACT_SELF_RESTRAINT_LIMIT_M, //int self_restraint_limit,行動を自粛する閾値(緊急事態宣言の条件)
0, //int personality_Influencer,活動的な人物の割合
0, //int personality_alone,インドアな人物の割合
10000, //int personality_middle,平均的行動傾向人物の割合
30000, //int personality_Influencer_weight,インフルエンサーの能力強化値
3000, //int personality_alone_weight インドアな人物の能力弱体化値
0, //int immunity,免疫を持つ人の割合
100, //int modulo_M,クラスタ分けのための法
10000 //int mod_risk_wait;//異なるクラスタ間で感染するための通過率。万分率で指定。10000で通常感染、0で完全遮蔽。
);
return ret;
}
/***************************************************************
*
* ========目的========
* actions_property_set() の設定のサンプル。main関数側で直接指定しても良いが、
* 少し長くなるのと、パラメーターが多いので指定方法のサンプルとしていくつか定義する
* 非常に無防備な集団の設定。
*
*/
int actions_property_set_ex2( ACTIONS *actions ){
int ret;
if( actions==NULL ) return 0;
ret = actions_property_set( actions,
5000, //int mask,マスクをつける確率
0, //int cons_strong,無敵個体の、全人口に対する割合
0, //int cons_hide,隠れ感染者の、全人口に対する割合
9800, //int cons_middle,通常感染者の、全人口に対する割合
100, //int cons_severe,重症化する人間の、全人口に対する割合
100, //int cons_death,感染後14日目に死ぬ人間の、全人口に対する割合
ACT_SELF_RESTRAINT_M, //int self_restraint,行動を自粛する割合。
ACT_SELF_RESTRAINT_LIMIT_M, //int self_restraint_limit,行動を自粛する閾値(緊急事態宣言の条件)
0, //int personality_Influencer,活動的な人物の割合
0, //int personality_alone,インドアな人物の割合
10000, //int personality_middle,平均的行動傾向人物の割合
30000, //int personality_Influencer_weight,インフルエンサーの能力強化値
3000, //int personality_alone_weight インドアな人物の能力弱体化値
0, //int immunity,免疫を持つ人の割合
100, //int modulo_M,クラスタ分けのための法
10000 //int mod_risk_wait;//異なるクラスタ間で感染するための通過率。100分率で指定。100で通常感染、0で完全遮蔽。
);
return ret;
}
/***************************************************************
*
* ========目的========
* actions_property_set() の設定のサンプル。main関数側で直接指定しても良いが、
* 少し長くなるのと、パラメーターが多いので指定方法のサンプルとしていくつか定義する
* 非常に無防備な集団の設定。
*
*/
int actions_property_set_ex3( ACTIONS *actions ){
int ret;
if( actions==NULL ) return 0;
ret = actions_property_set( actions,
5000, //int mask,マスクをつける確率
0, //int cons_strong,無敵個体の、全人口に対する割合
0, //int cons_hide,隠れ感染者の、全人口に対する割合
9800, //int cons_middle,通常感染者の、全人口に対する割合
100, //int cons_severe,重症化する人間の、全人口に対する割合
100, //int cons_death,感染後14日目に死ぬ人間の、全人口に対する割合
ACT_SELF_RESTRAINT_M, //int self_restraint,行動を自粛する割合。
ACT_SELF_RESTRAINT_LIMIT_M, //int self_restraint_limit,行動を自粛する閾値(緊急事態宣言の条件)
0, //int personality_Influencer,活動的な人物の割合
0, //int personality_alone,インドアな人物の割合
10000, //int personality_middle,平均的行動傾向人物の割合
30000, //int personality_Influencer_weight,インフルエンサーの能力強化値
3000, //int personality_alone_weight インドアな人物の能力弱体化値
0, //int immunity,免疫を持つ人の割合
200, //int modulo_M,クラスタ分けのための法
500 //int mod_risk_wait;//異なるクラスタ間で感染するための通過率。100分率で指定。100で通常感染、0で完全遮蔽。
);
return ret;
}
/***************************************************************
*
* ========目的========
* actions_property_set() の設定のサンプル。main関数側で直接指定しても良いが、
* 少し長くなるのと、パラメーターが多いので指定方法のサンプルとしていくつか定義する
* 非常に無防備な集団の設定。
*
*/
int actions_property_set_ex4( ACTIONS *actions ){
int ret;
if( actions==NULL ) return 0;
ret = actions_property_set( actions,
6000, //int mask,マスクをつける確率
0, //int cons_strong,無敵個体の、全人口に対する割合
0, //int cons_hide,隠れ感染者の、全人口に対する割合
9800, //int cons_middle,通常感染者の、全人口に対する割合
100, //int cons_severe,重症化する人間の、全人口に対する割合
100, //int cons_death,感染後14日目に死ぬ人間の、全人口に対する割合
8000, //int self_restraint,行動を自粛する割合。デフォルトでは8000（80%）
ACT_SELF_RESTRAINT_LIMIT_M, //int self_restraint_limit,行動を自粛する閾値(緊急事態宣言の条件)
0, //int personality_Influencer,活動的な人物の割合
0, //int personality_alone,インドアな人物の割合
10000, //int personality_middle,平均的行動傾向人物の割合
30000, //int personality_Influencer_weight,インフルエンサーの能力強化値
3000, //int personality_alone_weight インドアな人物の能力弱体化値
0, //int immunity,免疫を持つ人の割合
50, //int modulo_M,クラスタ分けのための法
500 //int mod_risk_wait;//異なるクラスタ間で感染するための通過率。100分率で指定。100で通常感染、0で完全遮蔽。
);
return ret;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するデバッグ用シミュレーションSimulation_Run
* のお試し1
*
*/
int Simulation_Run1_ex1(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, nday;
ACTIONS actions;
_ftprintf(fp_glp, _T("パラメータ\t経過日数\t世界線\t数\r\n") );
nday = 365*4;
n = 3000;
ret = actions_property_set_ex1( &actions );
if( ret==0 ) return 0;
for(world_line=0;world_line<20;world_line++){
//for(world_line=0;world_line<1;world_line++){
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run( n, nday, fp_log, fp_glp, seed, &actions, world_line );
seed++;
}
return 1;
}
/***************************************************************
*
* ========目的========
* 条件を変えながら世界線をやり直すシミュレーションSimulation_Run2
* のお試し1
*
*/
int Simulation_Run2_ex1(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, nday;
ACTIONS actions;
_ftprintf(fp_glp, _T("パラメータ\t経過日数\t世界線\t数\r\n") );
nday = 365*4;
ret = actions_property_set_ex2( &actions );
if( ret == 0 ) return 0;
for(n=200;n<=10000;n=n+200){
for(world_line=0;world_line<10;world_line++){
Simulation_Run( n, nday, fp_log, fp_glp, seed, &actions, world_line );
seed++;
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するデバッグ用シミュレーションSimulation_Run3
* のお試し1
*
*/
int Simulation_Run3_ex1(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, nday;
ACTIONS actions;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex3( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t流行数\t平均患者数\r\n") );
for(n=500;n<=10000;n=n+500){
_ftprintf(fp_log, _T("ウイルスタイプ4感染数\t") );
_ftprintf(fp_log, _T("人口%d\t"), n );
if(n==500){
_ftprintf(fp_log, _T("=COUNTIF(E2:GV2,\">5\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(E2:GV2)\t") );
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<200;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run3( n, nday, fp_log, fp_glp, seed, &actions, world_line );
//seed=seed+3;
}
_ftprintf(fp_log, _T("\r\n"), n );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するデバッグ用シミュレーションSimulation_Run3
* のお試し2
* クラスタ分けを細かくしていく
*
*/
int Simulation_Run3_ex2(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, nday;
ACTIONS actions;
int clstn;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex3( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t流行数\t平均患者数\r\n") );
for(clstn=10;clstn<=50;clstn=clstn+10){
actions.modulo_M = clstn;
for(n=500;n<=10000;n=n+500){
_ftprintf(fp_log, _T("V4感染数(%dクラスタ)\t"), clstn );
_ftprintf(fp_log, _T("人口%d\t"), n );
if(n==500){
_ftprintf(fp_log, _T("=COUNTIF(E2:GV2,\">5\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(E2:GV2)\t") );
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<200;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run3( n, nday, fp_log, fp_glp, seed, &actions, world_line );
//seed=seed+3;
}
_ftprintf(fp_log, _T("\r\n"), n );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するシミュレーションSimulation_Run4
* のお試し1
*
*
*/
int Simulation_Run4_ex1(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, nday;
ACTIONS actions;
int clstn;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex4( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t流行数\t平均患者数\r\n") );
for(clstn=10;clstn<=10;clstn=clstn+10){
//actions.modulo_M = clstn;
for(n=10000;n<=10000;n=n+500){
_ftprintf(fp_log, _T("V4感染数(%dクラスタ)\t"), clstn );
_ftprintf(fp_log, _T("人口%d\t"), n );
if(n==500){
_ftprintf(fp_log, _T("=COUNTIF(E2:GV2,\">5\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(E2:GV2)\t") );
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<10;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run4( n, nday, fp_log, fp_glp, seed, &actions, world_line );
//seed=seed+3;
}
_ftprintf(fp_log, _T("\r\n"), n );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するシミュレーションSimulation_Run5
* のお試し1
*
*
*/
int Simulation_Run5_ex1(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, i, nday;
ACTIONS actions;
int clstn;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex4( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t流行数\t平均患者数\r\n") );
for(clstn=100;clstn<=200;clstn=clstn+100){
actions.modulo_M = clstn;
//for(n=1000;n<=10000;n=n+1000){
for(i=5;i<=100;i=i+5){
n = i*100;
_ftprintf(fp_log, _T("医療崩壊の有無(%dクラスタ)\t"), clstn );
_ftprintf(fp_log, _T("人口%d\t"), n );
if(n==1000){
_ftprintf(fp_log, _T("=COUNTIF(E2:GV2,\">5\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(E2:GV2)\t") );
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<50;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run5( n, nday, fp_log, fp_glp, seed, &actions, world_line );
}
_ftprintf(fp_log, _T("\r\n"), n );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するシミュレーションSimulation_Run5
* のお試し2
* 医療崩壊リスクを調べるため
* クラスターのサイズを変化させていく
*
*/
int Simulation_Run5_ex2(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, i, nday;
ACTIONS actions;
int clstsize;
int FstFlg=1;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex4( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t\t流行数\t平均患者数\r\n") );
for(clstsize=1000;clstsize<=10000;clstsize=clstsize+500){
//for(n=1000;n<=10000;n=n+1000){
for(i=100;i<=100;i=i+10){
n = i*1000;
actions.modulo_M = n/clstsize;
_ftprintf(fp_log, _T("医療崩壊の有無(%dクラスタ)\t"), n/clstsize );
_ftprintf(fp_log, _T("人口%d\tクラスタサイズ%d\t"), n, clstsize );
if(FstFlg==1){
_ftprintf(fp_log, _T("=COUNTIF(F2:GV2,\"=1\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(F2:GV2)\t") );
FstFlg = 0;
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<20;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run5( n, nday, fp_log, fp_glp, seed, &actions, world_line );
}
_ftprintf(fp_log, _T("\r\n"), n );
fflush( fp_log );fflush( fp_glp );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
}
return 1;
}
/***************************************************************
*
* ========目的========
* 一定の条件下で計算の妥当性を検証するシミュレーションSimulation_Run5
* のお試し1
* 医療崩壊リスクを調べるため
* クラスター間の感染率を変化させていく
*
*/
int Simulation_Run5_ex3(FILE *fp_log, FILE *fp_glp, int seed ){
int ret;
int world_line;//何度目の世界線をシミュレートしているのか
int n, i, nday;
ACTIONS actions;
int mod_risk_wait;
int FstFlg=1;
nday = 365*4;
//n = 3000;
ret = actions_property_set_ex4( &actions );
if( ret==0 ) return 0;
srand( seed );
_ftprintf(fp_log, _T("\t\t\t流行数\t平均患者数\r\n") );
for(mod_risk_wait=700;mod_risk_wait<=1200;mod_risk_wait=mod_risk_wait+50){
//for(n=1000;n<=10000;n=n+1000){
for(i=100;i<=100;i=i+10){
n = i*1000;
actions.mod_risk_wait_m = mod_risk_wait;
_ftprintf(fp_log, _T("医療崩壊の有無\t") );
_ftprintf(fp_log, _T("人口%d(クラスタ数%d)\t"), n, actions.modulo_M );
_ftprintf(fp_log, _T("クラスタ間感染率%0.3f%\t"), (double)mod_risk_wait*100/PER_DEF );
if(FstFlg==1){
_ftprintf(fp_log, _T("=COUNTIF(F2:GV2,\"=1\")\t") );
_ftprintf(fp_log, _T("=AVERAGE(F2:GV2)\t") );
FstFlg = 0;
}
else{
_ftprintf(fp_log, _T("\t\t") );
}
for(world_line=0;world_line<50;world_line++){
//for(world_line=0;world_line<1;world_line++){
fprintf_daily_Calendar_setting_act(fp_glp, world_line, &actions, n);
_tprintf( _T("\n%d回目の世界\n"), world_line );
Simulation_Run5( n, nday, fp_log, fp_glp, seed, &actions, world_line );
}
_ftprintf(fp_log, _T("\r\n"), n );
fflush( fp_log );fflush( fp_glp );
//fclose( fp_log );
//fp_log = now_date_fopen( _T("複数回シミュレーション結果"), _T(".log") );
}
}
return 1;
}
以上です。