"""
ランダムな経済取引による富の分布シミュレーション。
概要:
多数の参加者間で行われるランダムな金銭交換が、どのように富の不均衡（指数分布）
を生み出すかをシミュレーションし、リアルタイムで可視化します。
詳細説明:
- 統計力学におけるエネルギー等分配や、経済物理学における所得分布のモデルです。
- 全員が同じ初期資金からスタートし、ランダムに選ばれた相手に一定額（vtrade）を支払います。
- 支払う資金がない場合は、持っている全額を支払います（負債はなし）。
- 多数の取引サイクルを経て、以下の3つの観点から解析を行います：
1. 個人の資産推移（カオスな状態）
2. ソートされた資産分布（格差の可視化）
3. 資産の頻度分布（理論的な指数分布 $f(x) = A e^{-x/\langle x \rangle}$ への収束）
関連リンク: :doc:`randomtrade_usage`
"""
import sys
import random
import numpy as np
from math import exp
import matplotlib.pyplot as plt
#=============================
# デフォルトパラメータ
#=============================
NP_DEF = 50 # 参加人数
VAL_DEF = 50 # 初期資金（平均額）
VTRADE_DEF = 1 # 1回の取引額
NITER_DEF = 500 # 最大反復サイクル数
NPLOT_DEF = 100 # グラフ更新間隔
NFX_DEF = 21 # 頻度分布のメッシュ数
[ドキュメント]
def main():
"""
シミュレーションのメインループ。
取引の実行、統計計算、グラフの更新を制御します。
"""
# 引数処理
argv = sys.argv
npersons = int(argv[1]) if len(argv) >= 2 else NP_DEF
avg_value = int(argv[2]) if len(argv) >= 3 else VAL_DEF
vtrade = int(argv[3]) if len(argv) >= 4 else VTRADE_DEF
niter = int(argv[4]) if len(argv) >= 5 else NITER_DEF
nplotinterval = int(argv[5]) if len(argv) >= 6 else NPLOT_DEF
nfx = int(argv[6]) if len(argv) >= 7 else NFX_DEF
# 初期化：全員同じ金額からスタート
x_index = np.arange(npersons)
v = np.full(npersons, float(avg_value))
# 理論的な指数分布（ボルツマン分布相当）の準備
# f(x) = (N/avg) * exp(-x/avg)
beta = 1.0 / avg_value
coeff_a = npersons * beta
ftx = np.linspace(0, avg_value * 10, 101)
fty = [coeff_a * exp(-beta * x) for x in ftx]
print(f"\nSimulation settings: n={npersons}, avg={avg_value}, trade={vtrade}")
print(f"Theoretical eq: f(x) = {coeff_a:.4g} * exp(-{beta:.4g} * x)")
# 描画準備
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 10))
plt.subplots_adjust(hspace=0.4)
# 取引シミュレーション開始
for i in range(niter):
# 1サイクル：全員がランダムに1人と取引
for ip in range(npersons):
itarget = random.randrange(npersons)
if itarget == ip: continue # 自分自身とは取引しない
# 支払額の決定（所持金が足りない場合は全額）
payment = vtrade if v[ip] >= vtrade else v[ip]
v[ip] -= payment
v[itarget] += payment
# 表示更新タイミングの判定
if i == 0 or (i + 1) % nplotinterval == 0 or i == niter - 1:
v_sorted = np.sort(v)
v_max = np.max(v)
ax1.cla()
ax2.cla()
ax3.cla()
# グラフ1: 生の資産分布
ax1.bar(x_index, v, color='skyblue', alpha=0.7)
ax1.axhline(avg_value, color='red', linestyle='--', linewidth=0.8, label='Initial/Avg')
ax1.set_title(f'Individual Wealth (Step {i+1}/{niter})')
ax1.set_ylabel('Wealth')
ax1.set_xlim(-1, npersons)
# グラフ2: 格差の可視化（ソート済み）
ax2.bar(x_index, v_sorted, color='steelblue', alpha=0.7)
ax2.axhline(avg_value, color='red', linestyle='--', linewidth=0.8)
ax2.set_title('Wealth Distribution (Sorted)')
ax2.set_ylabel('Wealth (Sorted)')
ax2.set_xlim(-1, npersons)
# グラフ3: 頻度分布と理論曲線
# ヒストグラムの計算
counts, bins = np.histogram(v, bins=nfx, range=(0, max(v_max, 1)), density=False)
bin_width = bins[1] - bins[0]
bin_centers = (bins[:-1] + bins[1:]) / 2.0
y_dist = counts / bin_width # 密度に変換
ax3.step(bin_centers, y_dist, where='mid', label='Simulation', color='black')
ax3.plot(ftx, fty, label='Exponential Theory', color='red', linestyle='dashed', alpha=0.7)
ax3.set_title('Probability Density (Wealth Function)')
ax3.set_xlabel('Wealth')
ax3.set_ylabel('Frequency Density')
ax3.set_xlim(0, max(v_max, avg_value * 3))
ax3.legend(fontsize=9)
plt.pause(0.05)
print("\nSimulation completed.")
plt.show(block=False)
input("Press ENTER to exit>>")
if __name__ == '__main__':
main()