using PyPlot

# plotcomparison(mcint) は1変数函数 f(x) を積分する函数 mcint(f) を与えると
#
# f(x) = exp(-x^2/a)
#
# の積分を mcint(f) で計算した結果とexactな値 sqrt(a*π) を比較するための
# グラフを表示してくれる.
#
function plotcomparison(mcint)
    as = logspace(-1, 1, 10^6)
    f(x) = exp(-x^2/a)
    @time s_mcint = [mcint(x->exp(-x^2/a)) for a in as]
    s_exact = [sqrt(a*π) for a in as]

    sleep(0.1)
    figure(figsize=(5,3.5))
    plot(as, s_mcint, label="MC integration")
    plot(as, s_exact, label="exact value", ls="--")
    xscale("log")
    grid()
    legend()
end

# 標準正規分布の sample を標準正規分布の確率密度函数の値の逆数に変換する函数
#
sample2weight(sample) = 1 ./ (x -> exp(-x^2/2)/sqrt(2π)).(sample)

# 大域変数達

srand(2018)
N = 2^6
sample = randn(N)
weight = sample2weight(sample);

# 大域変数を含む函数は遅い.

mcint_naive(f) = mean(weight .* f.(sample))
plotcomparison(mcint_naive);

# 大域変数を定数で置き換えると速くなる.　
# ただし, 定数を書き換えても函数は変化しなくなる.

const sample_const = sample
const weight_const = weight
mcint_const(f) = mean(weight_const .* f.(sample_const))
plotcomparison(mcint_const);

# function-like object を使っても速くなる.

mutable struct MCInt
    sample::Array{Float64,1}
    weight::Array{Float64,1}
end

(mcint::MCInt)(f) = mean(mcint.weight .* f.(mcint.sample))
mcint_struct = MCInt(sample, weight)
plotcomparison(mcint_struct);

# mutable function-like object であれば中身を書き換えることができる.

srand(5)
mcint_struct.sample = randn(N)
mcint_struct.weight = sample2weight(mcint_struct.sample)
plotcomparison(mcint_struct);

# keyword arguments は次のような使い方をできない.

mcint_kwargs(f; sample=sample, weight=weight) = mean(weight .* f.(sample))
plotcomparison(mcint_kwargs);

# 定数を使った keyword arguments を使うと速い.

mcint_kwargs(f; sample=sample_const, weight=weight_const) = mean(weight .* f.(sample))
plotcomparison(mcint_kwargs);

# しかし, 定数を変化させても mcint_kargs 函数は変化しない.

srand(5)
sample_const = randn(N)
weight_const = sample2weight(sample_const)
plotcomparison(mcint_kwargs);

# mcint_kwargsを再定義すると定数変更が有効になる.

mcint_kwargs(f; sample=sample_const, weight=weight_const) = mean(weight .* f.(sample))
plotcomparison(mcint_kwargs);