NSSOLAdvent Calendar 2023

日鉄ソリューションズ株式会社

いい意味でこんなはずじゃなかった、マイクラで物理シミュレーション

Last updated at 2023-12-25Posted at 2023-12-09

いい意味で、こんなはずじゃなかった

この記事を読むことで、マイクラ世界で一風変わった物理シミュレーションをする方法と、そこで得られたエンジニアリングに関する学びについて知ることがができます。

この記事はマイクラ内の物理シミュレーションMODを使って、ビュフォンの針を用いた円周率計算を行う話です。こう聞くと簡単な話に思えてしまいますよね。

私も当初は、一様乱数を作ればちょちょいと再現できるだろうと軽く考えていたのですが、実際に始めてみたら、次のような想定外の技術課題が次々と現れてきて、結果的にエンジニアリングの楽しさを味わえるお題になりました。

今回、進んではハマりつつ得られた学び：

健気な針たちの様子。マイクラでもこんなものが作れます：

以下、実験方法と学びについて説明していきます。

実験方法

ビュフォンの針

まず、今回、円周率の近似値計算に用いるビュフォンの針とは、地面に多数の平行線を引き、そこに針をランダムに落すと、どれかの線と針が交差する確率から円周率の近似値が求まるという数学の理論です。地面の平行線の間隔を $D$ に、針の長さを $L=D/2$ としたとき、交差する確率 $p$ は次の式で与えられます。

$$p = \frac{1}{\pi}$$

従って、落とした針の本数を $T$ 本、交差した本数を $C$ 本とすると円周率の近似値は次の式で与えられます。

$$円周率 = \frac{1}{p} \fallingdotseq \frac{T}{C} $$

（理論的には）落とした針の本数が多ければ多いほど近似が改善していきます。
以下の実験では平行線の間隔 $D=12$、針の長さ $L=6$ として計算しています。

マインクラフトの世界でビュフォンの針を実現する方法

実験環境準備

マインクラフトの世界内で次の図1のようなビュフォンの針を作り、ランダムにばら撒くのですが、当然ながらデフォルトのマインクラフトではその四角四面な座標系からはみ出すことはできません。

図1：巨人でも刺さらなさそうな直径1〜3m、長さ7mの針の勇姿。なお、コンピュータとエンジン（スラスター）も搭載しています。

そこで、マインクラフトの座標系からオブジェクトを解き放ち、実験を自動化するためにマインクラフトMODをいくつか利用します。今回、実験に用いた環境は次のとおりです。
※ここで、Java版マインクラフトとMODとは何なのかやインストール方法、注意点、トラブル解決法については知っていることを前提とします。ググれば大抵答えが得られますから。

MacBook Pro 2021
Oracle JDK 17.0.6
Java版マインクラフト 1.18.2
- 今回はPrism Launcherを用いて構成管理しました。
マインクラフトMOD
- Minecraft Forge 40.2.0: MODの基盤となるシステムMOD
- Valkyrien Skies 2.1.0-beta.12(Forge): 物理シミュレーションMOD
- Valkyrien Tournament 1.0.0-beta1(Forge): スラスターなどを追加するValkyrien Skiesのアドオン
- Eureka! 1.1.0-beta.9(Forge): 船（今回は針）を作るための舵などを追加するValkyrien Skiesのアドオン
- Kotlin for Forge 3.4.0: Valkyrien Skiesの依存ライブラリ
- CC:Tweaked 1.101.3(Forge): マイクラ内にLua言語でプログラミングできるコンピュータを追加するMOD
- Valkyrien Computers 1.5.2(Forge): Valkyrien Skiesの船（今回は針）の位置と角度を測定するセンサーなどを追加するCC:Tweakedのアドオン
ゲーム生成時の設定
- ゲームモード：クリエイティブ
- 構造物の生成：オフ
- ワールドタイプ：スーパーフラット
ゲーム生成後のコマンドによるゲーム設定変更
- 次のコマンドを実行して、時間と天気が変わらないようにする。
```
/gamerule doDaylighteCycle false
/time set noon
/gamerule doWeatherCycle false
/weather clear
```
- 次のコマンドを実行して、モブが湧かさないようにする。
```
/gamerule doMobSpawning false
```

針の構築方法

まずブロックを組み合わせて針（マイクラの座標系からはみ出せるValkyrien Skiesの船）を作ります。
次のムービーのように、鉄ブロックとスラスター（Thruster）、コマンドコンピューター、船の舵（Oak Ship Helm）、船リーダー（Ship Reader）を組み合わせて針を作ります。

注意点：
- Ship Readerを、隣にある似た名前のShip Raderと間違えないように。
- スラスターの向きはお互い逆向きになっていることに注意してください。これは針を回転させるために使うからです。
- コマンドコンピューターの上に舵（Oak Ship Helm）を置くには、しゃがみながら行う必要があります。しゃがまないとコンピュータの操作画面に入ってしまいます。コンピュータの操作画面からはESCキーを押すことで抜けられます。
- ムービーだと分かりにくいですが、最後の方でコマンドコンピューターの上に舵を置いた後にしていることは、次の通りです。
  1. 手を空にします。手に何かを持っていると舵を操作できないためです。
  2. しゃがんだまま舵を右クリックすることで、船の組み立てメニューを表示します。
  3. Assembleボタンをクリックすることで、船を組み立てます。これでマイクラの座標系から船がはみ出せるようになります。ESCキーを押して組み立てメニューから抜けます。
  4. 舵を壊します。壊さないと期待通りには針が落下しません。

このままだと太すぎるので、計算上、ビュフォンの針とみなすのは赤い線の部分（長さ6ブロック分）で、地面の線は青の線の部分です。この赤の線と青の線が交差したら、針が地面の並行線と交差したとみなします。

次に、針（船）とコンピュータに特定に名前を付けることと、コンピュータのプログラム格納用フォルダが作られる操作を行います。
次のムービーを参考に、以下のステップを実行してください。

ステップ：
1. 船のリネームコマンド /vs ship AAAA-AAAA-AAAA rename BBBB を実行して、針に名前を付けます。
  ここでAAAA-AAAA-AAAAは自動的に付けられた名前、BBBBは新たな名前です。ムービー内ではneedle14という新たな名前を付けています。針の名前は、針をテレポートさせるときに必要となります。
2. コマンドコンピュータ（の頭）を右クリックして起動し、label set CCCC コマンドを実行してコンピュータに名前を付けます。
3. 何かファイルを作って保存することで、マイクラを動かしている方のコンピュータ（以下、ホスト）にプログラム格納用フォルダを作ります。ムービー内では edit needle_comp4.txt（タイポ恥ずかしい・・・）コマンドを使い、エディタを開き中身が空のファイルを保存しています。（エディタの下に書いてありますが）エディタの保存や終了のメニューはエディタ内でCtrlキーを押すと表示されます。なお、コンピュータから抜けるにはESCキーを押します。

プログラムの準備

以上により、ホスト側のマイクラデータ格納フォルダにプログラム格納用フォルダが作られていますので、それを探します。
Prism Launcherを使っている場合は次の手順で探せます。

今回のマインクラフトインスタンスを右クリックして「編集...」を選びます。
ウィンドウの左の「ワールド」をクリックします。
ウィンドウの右下の「フォルダーを表示」をクリックします。
マイクラのセーブデータフォルダ（saves）が開くので、次のパスを辿ると目的のプログラム格納用フォルダに到達できます。
saves/{実験に用いているワールド名}/computercraft/computer/{先ほど作ったファイル（needle_comp4.txt）が格納された1からの連番のフォルダ}

そのディレクトリ直下に次の2つのファイルを格納します。

drop_randomly.lua: 針を上空からランダムに落下させ、その位置と角度情報をログファイルに書き出すプログラム（Lua言語）
env.lua: drop_randomly.lua の設定ファイル（Lua言語）

設定ファイル env.lua の内容を適宜変更します。

env.lua

constants = {
	SHIP_NAME = 先ほど設定した針（船）の名前（例:"needle"）,
	ORIGIN = { x=線のx座標（例:0.5）, y=地面の一つ上のz座標（例:-58.5）, z=線のz座標（例:0.5） },
	LOG_FILENAME = ログファイル名（例:"/buffons_needle.log"）,
	L = 針の長さ（例：6）,
	D = 線と線の間の距離（例:12）,
	DROP_W = ランダムに上空でブロックを落とす幅の半分（例:12 / 2 + 1.5）,
	MIN_DROP_HEIGHT = 針を落とす地面からの高さ（例:20）,
	DROP_HEIGHT_RANGE = ランダムに針を落とす地面からの高さの範囲（例:10）,
	KEEP_THRUSTER_RANGE = 地面に近づいた後にランダムにスラスターの噴射を維持する範囲（例:3）,
	NEAR_GROUND_HEIGHT = 地面に近づいたとみなす地面からの距離（例:3）
}

drop_randomly.lua はビュフォンの針の実験を行うプログラムです。
針のコンピュータを右クリックして、プロンプトで drop_randomly コマンドを実行すれば針が宙空にテレポートしてランダムな位置に降ってくるようになります。頭上に注意しつつ、速やかにその場から離れてください。

このプログラムは次のような処理を行います。
ざっくり言うと、何度も何度も、落ちた針を拾っては上空のランダムな位置から針を落として、ランダムに回転させることでビュフォンの針に求められる位置と角度に関するランダムな分布を得るという仕組みです。

47〜51行目: ログファイルがなければ、ログファイルを作成します。
61〜64行目: 上空のランダムな位置から針を落とします。落とす幅はx軸について±DROP_W。高さはMIN_DROP_HEIGHTからMIN_DROP_HEIGHT+DROP_HEIGHT_RANGEの範囲。
66行目: スラスターの噴射を開始し、針を回転させます。
68〜74行目: 高さがNEAR_GROUND_HEIGHTになるまで待ちます。
76〜79行目: ランダムに最大KEEP_THRUSTER_RANGE秒待って、スラスターを停止します。
81〜89行目: 針の速度が0.001以下になるまで待ちます。
91〜107行目: 簡易的に針一本を用いた円周率の近似計算をしていますが、本番の円周率計算には用いないため気にしなくて大丈夫です。
109〜117行目: ログファイルに停止した針の中央位置（x, y, z）や角度（yaw, roll, pitch）などを出力します。

drop_randomly.lua

require("env")

local SHIP_NAME = constants.SHIP_NAME
local ORIGIN = constants.ORIGIN
local LOG_FILENAME = constants.LOG_FILENAME
local L = constants.L
local D = constants.D
local DROP_W = constants.DROP_W
local MIN_DROP_HEIGHT = constants.MIN_DROP_HEIGHT
local DROP_HEIGHT_RANGE = constants.DROP_HEIGHT_RANGE
local KEEP_THRUSTER_RANGE = constants.KEEP_THRUSTER_RANGE
local NEAR_GROUND_HEIGHT = constants.NEAR_GROUND_HEIGHT


local reader = peripheral.wrap("back")


function start_thruster()
	redstone.setOutput("left", true)
	redstone.setOutput("right", true)
end

function stop_thruster()
	redstone.setOutput("left", false)
	redstone.setOutput("right", false)
end

function getSpeed(reader)
	local v = reader.getVelocity()
	return math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
end

function range(n)
	local i = 0
	return function()
		if i >= n then
			return nil
		end

		local r = i
		i = i + 1
		return r 
	end
end


if not fs.exists(LOG_FILENAME) then
	local log_file = fs.open(LOG_FILENAME, "a")
	log_file.write("date,origin_x,origin_y,origin_z,teleport_x,teleport_y,ship_x,ship_y,ship_z,ship_yaw,ship_roll,ship_pitch,pi\n")
	log_file.close()
end


local cross_count = 0
local total_count = 0

stop_thruster()

-- for i in range(1000) do
while true do
	local teleport_x = math.random() * DROP_W * 2 - DROP_W + (ORIGIN.x - 0.5)
	local teleport_y = math.random() * DROP_HEIGHT_RANGE + MIN_DROP_HEIGHT + (ORIGIN.y - 0.5)
	commands.exec("vs teleport " .. SHIP_NAME .. " " .. 
		teleport_x .. " " .. teleport_y .. " " .. (ORIGIN.z - 0.5))

	start_thruster()

	while true do
		local ship_pos = reader.getWorldspacePosition()
		if ship_pos.y <= ORIGIN.y + NEAR_GROUND_HEIGHT then
			break
		end
		sleep(0.1)
	end

	local keep_thruster_sec = math.random() * KEEP_THRUSTER_RANGE
	sleep(keep_thruster_sec)

	stop_thruster()

	while true do
		local speed = getSpeed(reader)
		if speed <= 0.001 then
			break
		end

		sleep(0.1)
	end
	sleep(0.5)

	local ship_pos = reader.getWorldspacePosition()
	local ship_rotation = reader.getRotation()
	local sin_roll = math.sin(ship_rotation.roll)
	local left_x = ship_pos.x + L / 2 * sin_roll
	local right_x = ship_pos.x - L / 2 * sin_roll

	if (ship_pos.x - ORIGIN.x) >= -D / 2 and (ship_pos.x - ORIGIN.x) < D / 2 then
		total_count = total_count + 1
		if (left_x - ORIGIN.x) * (right_x - ORIGIN.x) < 0 then
			cross_count = cross_count + 1
		end
	end
	
	local pi_approx = (total_count / cross_count)
	if cross_count == 0 then
		pi_approx = 0.0
	end

	local log_file = fs.open(LOG_FILENAME, "a")
	log_file.write(os.date("%F %T") .. 
		"," .. ORIGIN.x .. "," .. ORIGIN.y .. "," .. ORIGIN.z ..
		"," .. teleport_x .. "," .. teleport_y ..
		"," .. ship_pos.x .. "," .. ship_pos.y .."," .. ship_pos.z .. 
		"," .. ship_rotation.yaw .. "," .. ship_rotation.roll .. "," .. ship_rotation.pitch .. 
		"," .. pi_approx ..
		"\n")
	log_file.close()
end

-- print("total: " .. total_count .. ", crossed: " .. cross_count)
-- print("pi = " .. (total_count / cross_count))

得られたエンジニアリングの学び

上記のプログラムは最終的な完成形なんですけれども、ここに至るまでにも色々ありましてね。
そう来たかという学びを時系列順に共有させていただきます。何かの気付きになったら幸いです。

1. 仮想空間で計算機リソース不足問題に出くわすなんて・・・。

休日に数時間試行錯誤して針を数本完成させた後、「さーて、ブロックもプログラムもできたから、一晩、針に頑張ってもらうとするよ！」と家族に呟いてMacBookを閉じたのですが、今思えばフラグ発言だったか・・・

で、一晩明けてみて「さて、どれだけデータ量が増えているかな」と確認してみたら、「おや、ログが増えない？？」と。
どの針のログファイルのサイズも950KB近辺にあって、「これはもしや」とコンピュータMOD（CC:Tweaked）の仕様を確認してみたら、デフォルトのドライブサイズの上限は1MBでした（え、いまどき1MB？？1GBじゃなく？？）。ディスクフルになっていたということですね、ださい・・・。

仮想空間のコンピュータなら、少なくともホストの容量が許すくらいはドライブが使えると思い込んでおりました。
思い込みはだめで、ちゃんと計算リソースのキャパシティは確認しましょうということですね、はい。

さて、解決策ですが、設定ファイルでドライブの最大容量を定義できますので、そこを緩和すれば完了です。
具体的には先ほども出てきたマイクラのセーブデータフォルダ（saves）配下にある saves/{実験に用いているワールド名}/serverconfig/computercraft-server.toml でドライブの最大容量を設定できます。ドライブの最大容量は computer_space_limit という項目です。

デフォルトでは次のように1MBになっております。

これを次のように三桁足せばドライブの容量を1GBにすることができます。これで一年以上動かし続けても大丈夫になります。

2. 乱数が乱数じゃない問題（その1）：一様乱数になるはずが、あれ、なんか釣鐘型なんですけど・・・。

ディスク容量も大丈夫になったことだしとデータ量が増えてきたので円周率の近似値計算に入りました。
しかし、ここで、またひと問題が発覚。

当初は、上空から針を落とすときに地面の線からみて垂直方向に-6から6の範囲でランダムに落とすようにしておりました（最終版のプログラムでいうと env.lua　で DROP_W = 6 という設定にしたものと同じ条件で実験していました）。
理想的には、落とした針は線からの距離（diff_x = ship_x - origin_x）は-6から6の範囲で一様分布となってほしいのですが、実際の距離はChatGPT先生のお力をお借りしたところ、下図のようになりました。

お、おう、なんだか釣鐘型ですね。値の範囲も理想では-6〜6になってほしいところが、マイナス方向には2くらいはみ出し、逆にプラス方向は2近く足りないと・・・。全然ダメでした。

一様乱数性に関するKS検定のp値も、社会人になって（いや、この人生で）見た一番小さい値 $2.96\times 10^{-100}$ を指し示していました（p値が0.05を下回ったら高い確率で分布が一様ではないと言える（帰無仮説を棄却できる）ということ）。まあp値をみるまでもないですけどね。

上空で針を落とした時点ではx座標は一様分布になっているのですが、針の様子をよく観察すると地面で針がくるくる回る時にスライドする現象が起きているためでした。他には落下中にもスライドしている可能性があります。まあ物理シミュレーション自体にランダム性はありますからね。

解決策ですが、多少の無駄は発生しますが、上空から落とす時は-7.5〜7.5のより広い範囲で落とし、落ちて回転が停止した後の針の位置が線から-6〜6にあるもののみデータとして採用することしました。

それにしても、こういう、ちょっと軽く分析してみようというときにChatGPT先生の存在はありがたいですね。チャチャっと試せて。

3. 乱数が乱数じゃない問題（その2）：え？仮想空間内の針なのに不良品が紛れ込むの？？

さーて、これで乱数の問題は解決と思いきや、まだ終わってはいませんでした。
数日シミュレーションを回した後、13本の針を放り投げた結果データをマージしたものが一様乱数になっているかKS検定してみたところ、何だか針は線からの距離（diff_x = ship_x - origin_x。ただし、-6〜6にあるレコードのみ採用）と針の角度（ship_roll）の一様乱数性の検定結果のp値が0.1〜0.2と低めであることが気になりました。

そこで、マージ前の各針についても一様乱数性をKS検定してみたところ、次のようにp値が低い針がちらほらとあることが分かりました。例えば針 #1 は diff_x のp値が0.0598であり、高い確率で一様分布から外れているといえます。
正直、ヒストグラムを見ても一様分布に近いのか遠いのか分かりませんでした。検定は大事ですね。

OK/NG	針番号	diff_xのp値	ship_rollのp値
NG	#1	0.0598	0.533
NG	#2	0.429	0.456
NG	#3	0.859	0.4
NG	#4	0.0975	0.234
OK	#5	0.976	0.548
OK	#6	0.7	0.904
OK	#7	0.613	0.921
OK	#8	0.803	0.558
OK	#9	0.673	0.515
NG	#10	0.398	0.682
NG	#11	0.741	0.324
NG	#12	0.52	0.294
OK	#13	0.597	0.785
...	...	...	...

一番左の列 OK/NG は、やや保守的にp値が仮に両方とも0.5以上であればOK（一様分布ではないとは言い切れない）、未満はNGと判定したものです。この基準では約半数の針がNG（不良品）ということです。

実験結果の章で出てきますが、このように不良品の針を除外することによりマージ後のp値が0.6〜0.8と許容できる数値になりました。

しかし、仮想世界で作った針に不良品が混入するなんて・・・。
不良品が生まれる可能性として一つ考えられるのは、よくある話ですが、乱数生成アルゴリズムの問題ですね。

マイクラ内のコンピュータのプログラミング言語であるLuaの乱数アルゴリズムが何か調査してみました。
Lua 5.4から用いられている乱数生成アルゴリズムはxoshiro256**であり乱数性に問題ありませんが、マイクラ内のコンピュータが採用しているLuaのバージョンは5.2/3であり、そこではstdlibのrandom関数（線形合同法の類の乱数性がダメなもの）を利用しています。
乱数アルゴリズムが原因である可能性があるということです。

今回は締切の関係で試せませんでしたが、Luaでxoshiro256**系アルゴリズムを実装して一様乱数性が改善するか実験してみるといいですね。

4. そして、運用でカバー（苦笑）

これで、もう、問題はないだろうと思っていたら、次がありました。
多数の針（少なくとも10本以上）を同時に動かしていると、時々、物理シミュレーションがフリーズする（時間が止まったかのように針が空中で止まっている。人は動ける。）という現象が発生するという問題です。
間隔が短い時には3時間くらいで発生し、長いときは12時間くらいで発生したりします。タイミングがランダムに思えます。

ご参考程度にフリーズしたときのログ

おそらくスラスターの噴射処理あたりで、ConcurrentModificationExceptionが発生しているようです。

[08Dec2023 17:28:36.539] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: Error in physics pipeline background task
java.util.ConcurrentModificationException: null
	at java.util.ArrayList$Itr.checkForComodification(ArrayList.java:1013) ~[?:?]
	at java.util.ArrayList$Itr.next(ArrayList.java:967) ~[?:?]
	at org.valkyrienskies.tournament.ship.tournamentShipControl.applyForces(tournamentShipControl.java:163) ~[tournament-1.0.0-beta2+39af2d30ae-forge.jar%2368!/:?]
	at org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineStage.tickPhysics(VSPhysicsPipelineStage.kt:138) ~[VS2_Nightly_Forge.jar%2370!/:?]
	at org.valkyrienskies.core.impl.pipelines.VSPipelineImpl.tickPhysics(VSPipelineImpl.kt:109) ~[VS2_Nightly_Forge.jar%2370!/:?]
	at org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask.run(VSPhysicsPipelineBackgroundTask.kt:75) [VS2_Nightly_Forge.jar%2370!/:?]
	at org.valkyrienskies.core.impl.pipelines.VSPipelineImpl$physicsThread$1.invoke(VSPipelineImpl.kt:65) [VS2_Nightly_Forge.jar%2370!/:?]
	at org.valkyrienskies.core.impl.pipelines.VSPipelineImpl$physicsThread$1.invoke(VSPipelineImpl.kt:64) [VS2_Nightly_Forge.jar%2370!/:?]
	at kotlin.concurrent.ThreadsKt$thread$thread$1.run(Thread.kt:30) [kotlinforforge-3.4.0-obf.jar%2365!/:3.4.0]
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/ERROR] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: !!!!!!! VS PHYSICS THREAD CRASHED !!!!!!!
[08Dec2023 17:28:36.540] [Physics thread/WARN] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineBackgroundTask/]: Physics pipeline ending
[08Dec2023 17:28:41.538] [Server thread/WARN] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineStage/]: Too many game frames in the game frame queue. Is the physics stage broken?
[08Dec2023 17:28:42.541] [Server thread/WARN] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineStage/]: Too many game frames in the game frame queue. Is the physics stage broken?
[08Dec2023 17:28:43.553] [Server thread/WARN] [org.valkyrienskies.core.impl.pipelines.VSPhysicsPipelineStage/]: Too many game frames in the game frame queue. Is the physics stage broken?

で、締切も迫っていたし、原因は根深そうな気配を感じたので、よくある再起動ソリューションで解決です（苦笑）。
マイクラを一度終了して、再度起動、また drop_randomly プログラムを起動して回るという暫定対策です。

これの原因調査と根本対策は残された宿題です。重そう・・・。
せめてフリーズを検知して、パトランプを光らせるくらいの仕組みを導入したいところ。

実験結果

さて、実験結果なのですが、前章の課題が解決した今、簡単です。

針は線からの距離（diff_x = ship_x - origin_x。ただし、-6〜6にあるレコードのみ採用）と針の角度（ship_roll）の一様乱数性の検定結果のp値が両方とも0.5以上のレコードのみを残して、マージした結果を用いて次の計算式を用いて円周率の近似値を求めるだけです。
$$円周率 \fallingdotseq \frac{T}{C}$$ ただし、$T$ は落とした針の本数、$C$ は線と交差した針の本数です。

以下、次のNotebookから結果のみを抜粋します。冗長ですので詳細はNotebookをご確認ください。

円周率の近似値計算の詳細: https://github.com/sonoisa/misc/blob/main/minecraft/buffons_needle/filter_data_and_estimate_pi.ipynb

不良品のフィルタリング

プログラムを一時停止したり、針の物理シミュレーションがフリーズするたびにログファイルを分割していたため見かけ上、針は56本あります（ログファイルが56ファイルあるということです）。

やや保守的に、diff_xとship_rollの一様乱数性のKS検定のp値が両方とも0.5以上であればOK、それ以外は怪しい不良品とみなしたところ56本中、（たった）12本が残りました。下表はその結果です。p値が0.1を切るものも結構見当たります。
実は#14以降は最終版の針の形状とは少し違う形状だったので、もしかしたらそれがNGを多く生む原因になったのかもしれません。

OK/NG	針番号	diff_xのp値	ship_rollのp値
NG	#1	0.0598	0.533
NG	#2	0.429	0.456
NG	#3	0.859	0.4
NG	#4	0.0975	0.234
OK	#5	0.976	0.548
OK	#6	0.7	0.904
OK	#7	0.613	0.921
OK	#8	0.803	0.558
OK	#9	0.673	0.515
NG	#10	0.398	0.682
NG	#11	0.741	0.324
NG	#12	0.52	0.294
OK	#13	0.597	0.785
NG	#14	0.0202	0.855
NG	#15	0.731	0.0592
NG	#16	0.00576	0.0654
NG	#17	0.0583	0.41
NG	#18	0.687	0.0103
NG	#19	0.142	0.379
NG	#20	0.244	0.944
NG	#21	0.774	0.213
OK	#22	0.792	0.738
NG	#23	0.747	0.0207
NG	#24	0.313	0.798
NG	#25	0.0337	0.636
NG	#26	0.0506	0.715
NG	#27	0.324	0.855
NG	#28	0.716	0.146
NG	#29	0.523	0.0467
NG	#30	0.315	0.702
NG	#31	0.301	0.0271
NG	#32	0.553	0.205
NG	#33	0.617	0.253
NG	#34	0.0303	0.249
NG	#35	0.891	0.395
NG	#36	0.337	0.526
NG	#37	0.0668	0.48
NG	#38	0.301	0.00298
NG	#39	0.0241	0.0527
NG	#40	0.177	0.558
OK	#41	0.933	0.65
NG	#42	0.305	0.165
NG	#43	0.246	0.0402
NG	#44	0.825	0.457
OK	#45	0.952	0.566
NG	#46	0.718	0.111
NG	#47	0.26	0.853
NG	#48	0.337	0.709
NG	#49	0.158	0.72
OK	#50	0.852	0.776
OK	#51	0.764	0.619
NG	#52	0.112	0.829
NG	#53	0.164	0.984
NG	#54	0.322	0.24
OK	#55	0.752	0.788
NG	#56	0.423	0.0207

これは針を落とした回数でいうと 37716 回です。

マージしたデータの一様乱数性

品質の大丈夫そうな針のデータをマージしたもののdiff_xとship_rollの一様乱数性のKS検定のp値はそれぞれ0.801、0.596です。大丈夫ですね。
また、ヒストグラムは下図のとおりです。

ちなみに、不良品を除外しない場合の diff_x のp値は 0.00998 という全然ダメな数値になります。ship_roll のp値は 0.800 というこちらは悪くない値でした。

円周率の近似値

線から針の中点までの相対位置 diff_x = ship_x - origin_x と平行線と針のなす角度 ship_roll から針の末端の座標は次のように計算できます。前述したように、針の末端は下図の赤い線の末端です。
$$\text{diff_x}\ + \frac{L}{2} * \sin(\text{ship_roll})　と　\text{diff_x}\ - \frac{L}{2} * \sin(\text{ship_roll})$$

針との交差を判定する対象の線は diff_x が0の位置にあるので、末端の符号が逆、もしくは少なくとも片方が0であれば交差していると判定できます。

今回でいうと交差した本数は 11995 本、落とした全数は 37716 本です。従って、
$$円周率 \fallingdotseq \frac{37716}{11995} \fallingdotseq 3.1443$$
という結果が得られます。おお、悪くないですね。

ちなみに不良品を除外しない場合は $円周率 \fallingdotseq 165410/52919 \fallingdotseq 3.1257$ というちょっと悪い結果になりました。このように不良品が混ざると低めの値になる傾向がありました。

収束状況

円周率の近似値の収束状況も確認してみます。
下図は、針を一本一本落とすたびの円周率の近似値の変化です。赤の点線は正解である円周率です。

ただし、レコードの順序をシャッフルするとこのチャートは大きく上下に変化しますので、これはあくまでも一例と思ってください。

また、各桁の収束状況もみてみます。円周率の近似値の整数部、小数一桁目、小数二桁目の数字が長く変わらなくなったら収束したといえるでしょうから。

まず整数部です。下図のとおり、400〜500回も針を投げれば整数部は3に収束することが分かります。

ティアキンのゾナウギアを使って針を作っても、比較的現実的な時間で、円周率の近似値として 3 くらいは得られるということですね。

次に小数一桁目です。20,000回ほど投げれば1に収束するようです。

最後に小数二桁目です。下図のようにまだまだ数字が変わっている状況にありますので、収束しているとはいえませんね。

以上の分析から、現時点の円周率の近似値は 3.1 と結論付けるのが良さそうです。

感想

当初はビュフォンの針のシミュレーションなんて簡単なお題と思っていのですが、こんなに色々と、これは仕事か？と思うようなエンジニアリング的課題が転がっているとは想像しておりませんでした。
頭の中で想像してできるだろうからOKとするのではなく、何かに挑めば色々と知識が得られるというお話ですね。

本当はもう少し時間があれば、針を上空に移動させる方法としてテレポートを用いるのではなく、針にスラスターを付けまくって全方向に推進できるドローンにしてPID制御を用いて所定の上空位置に移動させるつもりでした。それはまたの機会に。

謝辞

今年もアドカレで何か書こうと思っていたものの、いいアイディアが思い浮かばずにおりましたところ、アドカレ初日のyamamorisobaさんの記事「Minecraft Education Editionでコードビルダーを使って円周率を計算してみた」を見かけて、最近、子供とはまっているマイクラ用MODを複数組み合わせればビュフォンの針を用いた円周率計算ができるだろうと閃きました。
本記事作成のきっかけを与えてくれたyamamorisobaさんとMOD開発者さんたちに感謝いたします。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up