text:
# SPXL長期模擬データから，運用ケース別の成績を計算する
using CSV, DataFrames, Dates, Plots, StatsBase
pyplot(fmt=:svg, fontfamily="meiryo") # plotの設定

# SPXL読み込み
dfmt = dateformat"yyyy-mm-dd"
spxl = DataFrame(CSV.File("./spxl_pred_1954.csv", dateformat=dfmt, header=true))

# -------------------------------------
# (1)まずはSP500とSPXLで，30年の運用データの中央値，上位下位5%の価格推移を計算する
# -------------------------------------
# SPXLを30年データに分割
# 1954/1/4~1984/1/4, 1954/1/5~1984/1/5...1981/4/28~2021/4/28というイメージ
st_date = spxl.timestamp[1]
data_len = length(spxl.timestamp[spxl.timestamp .<= st_date + Dates.Year(30)])
en_date = spxl.timestamp[data_len]
# データ数data_len(7541日)分のデータを，data_num(9407)個用意
data_num = length(spxl.spxl) - data_len + 1
spxl_30y = ones(data_len, data_num)
for i in 1:data_num
    spxl_30y[:,i] = spxl.spxl[i:i + data_len - 1]
end

# データの正規化(1行目の値を1とした比率に変換)
for c = 1:data_num
    spxl_30y[:,c] ./= spxl_30y[1,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_spxl = zeros(data_len, 3)
for r in 1:data_len
    quan_spxl[r,:] = quantile(spxl_30y[r,:], [0.05, 0.5, 0.95])
end

labels = ["下位5%", "中央値", "上位5%", "投資額"]
plot( (1:data_len) .* 30 / data_len, quan_spxl[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="価格[\$]", yscale=:log10, ylims=(10^-0.5, 10^2.5))
plot!( (1:data_len) .* 30 / data_len, quan_spxl[:,2], label=labels[2], linecolor=:green, ms=1, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_spxl[:,1], label=labels[1], linecolor=:red, yscale=:log10)
savefig("SPXL_30y_5parcentile.png")

# sp500を30年データに分割--------------------
sp500_30y = ones(data_len, data_num)
for i in 1:data_num
    sp500_30y[:,i] = spxl.sp500_yidx[i:i + data_len - 1]
end

# データの正規化(1行目の値を1とした比率に変換)
for c = 1:data_num
    sp500_30y[:,c] ./= sp500_30y[1,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_sp500 = zeros(data_len, 3)
for r in 1:data_len
    quan_sp500[r,:] = quantile(sp500_30y[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_sp500[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="価格[\$]", yscale=:log10, ylims=(10^-0.5, 10^2.5))
plot!( (1:data_len) .* 30 / data_len, quan_sp500[:,2], label=labels[2], linecolor=:green, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_sp500[:,1], label=labels[1], linecolor=:red, yscale=:log10)
savefig("SP500_30y_5parcentile.png")

# -------------------------------------
# (2)上記の価格推移で，500万+毎月10万を積立投資した場合を計算する
# -------------------------------------

# 年月日から購入金額を計算
purchases = zeros(data_len)
purchases[1] = 500 # 初回500万円
mid_flag = false
for i = 2:data_len
    if (Dates.Month(spxl.timestamp[i - 1]) != Dates.Month(spxl.timestamp[i]))
        purchases[i] = 10
        mid_flag = false
    elseif ( (Dates.Month(spxl.timestamp[i]) == Dates.Month(6) || Dates.Month(spxl.timestamp[i]) == Dates.Month(12)) &&
        Dates.Day(spxl.timestamp[i]) > Dates.Day(15) && mid_flag == false)
        purchases[i] = 20
        mid_flag = true
    end
end

# 購入金額を積算して，総投資額も計算しておく
total_purchases = accumulate(+, purchases)

# 購入株式数
quantity_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    quantity_spxl[:,n] = purchases ./ spxl_30y[:,n]
end

# 購入株式数を積算して保有株式数とする
stocks_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    stocks_spxl[:,n] = accumulate(+, quantity_spxl[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_spxl = copy(stocks_spxl)
for c = 1:data_num
    assets_spxl[:,c] .*= spxl_30y[:,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_spxl_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_spxl_assets[r,:] = quantile(assets_spxl[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10, ylims=(10^2.5, 10^6.0))
plot!( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,2], label=labels[2], linecolor=:green, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,1], label=labels[1], linecolor=:red, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black, yscale=:log10)
savefig("SPXL_30y_5parcentile_assets.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string(quan_spxl_assets[end, i] * (1 - 0.20315)) * "万円\n")
end
quan_spxl_assets[2513, 2]

print("下位5%の最大下落率[%]:" * string(100 - minimum(quan_spxl_assets[:,1] ./ total_purchases) * 100)[1:4])

# SP500も計算-------------------------
# 購入株式数
quantity_sp500 = zeros(size(sp500_30y))
for n in 1:data_num
    quantity_sp500[:,n] = purchases ./ sp500_30y[:,n]
end

# 購入株式数を積算して保有株式数とする
stocks_sp500 = zeros(size(sp500_30y))
for n in 1:data_num
    stocks_sp500[:,n] = accumulate(+, quantity_sp500[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_sp500 = copy(stocks_sp500)
for c = 1:data_num
    assets_sp500[:,c] .*= sp500_30y[:,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_sp500_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_sp500_assets[r,:] = quantile(assets_sp500[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10, ylims=(10^2.5, 10^6.0))
plot!( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,2], label=labels[2], linecolor=:green, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,1], label=labels[1], linecolor=:red, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black, yscale=:log10)
savefig("SP500_30y_5parcentile_assets.png")
# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string(quan_sp500_assets[end, i] * (1 - 0.20315))[1:6] * "万円\n")
end
quan_sp500_assets[2513, 2] # 10年時点の中央値

print("下位5%の最大下落率[%]:" * string(100 - minimum(quan_sp500_assets[:,1] ./ total_purchases) * 100)[1:4])

# -------------------------------------
# (3)運用パターン変更(完全積立投資)
# -------------------------------------
# 5300万を全て積立投資する．30*12=360ヶ月あるので，14.7万円ずつ．

# 年月日から購入金額を計算
purchases = zeros(data_len)
amount = (5300 / 360)
purchases[1] = amount # 初回
for i = 2:data_len
    if (Dates.Month(spxl.timestamp[i - 1]) != Dates.Month(spxl.timestamp[i])) # 月が変わったタイミングで投資
        purchases[i] = amount
    end
end

# 購入金額を積算して，総投資額も計算しておく
total_purchases = accumulate(+, purchases)

# 購入株式数
quantity_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    quantity_spxl[:,n] = purchases ./ spxl_30y[:,n]
end

# 購入株式数を積算して保有株式数とする
stocks_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    stocks_spxl[:,n] = accumulate(+, quantity_spxl[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_spxl = copy(stocks_spxl)
for c = 1:data_num
    assets_spxl[:,c] .*= spxl_30y[:,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_spxl_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_spxl_assets[r,:] = quantile(assets_spxl[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,2], label=labels[2], linecolor=:green, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_assets[:,1], label=labels[1], linecolor=:red, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black, yscale=:log10)
savefig("SPXL_30y_5parcentile_assets_積立時.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string(quan_spxl_assets[end, i] * (1 - 0.20315))[1:7] * "万円\n")
end
quan_spxl_assets[2513, 2] # 10年時点の中央値

print("下位5%の最大下落率[%]:" * string(100 - minimum(quan_spxl_assets[:,1] ./ total_purchases) * 100)[1:4])

# SP500も計算-------------------------
# 購入株式数
quantity_sp500 = zeros(size(sp500_30y))
for n in 1:data_num
    quantity_sp500[:,n] = purchases ./ sp500_30y[:,n]
end

# 購入株式数を積算して保有株式数とする
stocks_sp500 = zeros(size(sp500_30y))
for n in 1:data_num
    stocks_sp500[:,n] = accumulate(+, quantity_sp500[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_sp500 = copy(stocks_sp500)
for c = 1:data_num
    assets_sp500[:,c] .*= sp500_30y[:,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_sp500_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_sp500_assets[r,:] = quantile(assets_sp500[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,3], label=labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,2], label=labels[2], linecolor=:green, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_sp500_assets[:,1], label=labels[1], linecolor=:red, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black, yscale=:log10)
savefig("SP500_30y_5parcentile_assets_積立時.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string(quan_sp500_assets[end, i] * (1 - 0.20315))[1:7] * "万円\n")
end
quan_sp500_assets[2513, 2] # 10年時点の中央値

print("下位5%の最大下落率[%]:" * string(100 - minimum(quan_sp500_assets[:,1] ./ total_purchases) * 100)[1:4])

# -------------------------------------
# (4)現金-SPXLリバランス
# -------------------------------------
# 最初に500万円投資，30年間現金は保有，50:50 or 70:30で毎月頭にリバランス
# ただし，平均購入金額を下回ったら売却はしない．

# 年月日から投資への現金追加額を計算
purchases = zeros(data_len)
purchases[1] = 500 # 初回500万円
mid_flag = false
for i = 2:data_len
    if (Dates.Month(spxl.timestamp[i - 1]) != Dates.Month(spxl.timestamp[i]))
        purchases[i] = 10
        mid_flag = false
    elseif ( (Dates.Month(spxl.timestamp[i]) == Dates.Month(6) || Dates.Month(spxl.timestamp[i]) == Dates.Month(12)) &&
        Dates.Day(spxl.timestamp[i]) > Dates.Day(15) && mid_flag == false)
        purchases[i] = 20
        mid_flag = true
    end
end
# 購入金額を積算して，総投資額も計算しておく
total_purchases = accumulate(+, purchases)

# リバランス計算を関数化
# price: 株価推移
# ratio: spxl対現金をいくつにするか．0.7なら70:30
# purchases: 追加投資現金
function rebalance(prices, ratio, purchases, timestamp, fee)
    # for debug
    # prices = spxl_30y[:,1]
    # ratio = 0.7
    # timestamp = spxl.timestamp
    # fee = 0.45/100
    # 変数定義と初期化
    tax_rate = 20.315 / 100 # 税率
    ave_unit_price = prices[1] # 平均取得単価
    hold_stock = zeros(length(prices)) # 保有株式数推移
    hold_stock[1] = purchases[1] * ratio / prices[1]
    hold_cash = zeros(length(prices)) # 保有現金推移
    hold_cash[1] = purchases[1] * (1 - ratio)
    ratio_temp = zeros(length(prices)) # リバランス時の比率
    ratio_temp[1] = ratio
    # 2日目以降を計算
    for i in 2:length(prices)
        # 現金を追加
        hold_cash[i] = purchases[i]
        # 月初めにリバランス
        if (Dates.Month(timestamp[i - 1]) != Dates.Month(timestamp[i]))
            # 現金と株式の比率を計算
            hold_stock_price = sum(hold_stock) * prices[i] # 現在の株式評価額
            hold_cash_price = sum(hold_cash) # 手持ちの現金
            # ratio_temp[i] = hold_stock_price / (hold_cash_price + hold_stock_price) # 現状の株式：現金比率
            diff = ratio * (hold_cash_price + hold_stock_price) - hold_stock_price # 理想の株式価格までの差額
            # 差額が正なら追加投資
            if (diff > 0)
                hold_stock[i] += diff / prices[i] # 差額分の数量の株を追加購入
                hold_cash[i] -= diff * (1 + fee) # 購入+手数料分の現金を引く
                ave_unit_price = ( (sum(hold_stock[1:i - 1]) * ave_unit_price) + diff ) / sum(hold_stock) # 平均取得単価の更新
                # print(ave_unit_price)
                # print("\n")
            # 差額が負なら売却
            else
                # 株価が平均取得単価以上のときだけ売却．下回っていたらホールド．税金がかかるので，その分考慮して多めに売却する
                if (prices[i] >= ave_unit_price)
                    sale = -diff # 売却額
                    profit = sale - sale / prices[i] * ave_unit_price
                    tax = profit * tax_rate # 税金
                    hold_stock[i] -= sale / prices[i] # 差額分の数量の株を売却
                    hold_cash[i] += sale * (1 - fee) - tax # 売却分(税額・手数料考慮)の現金を追加
                end
            end
            # 現金と株式の比率を計算
            hold_stock_price = sum(hold_stock) * prices[i] # 現在の株式評価額
            hold_cash_price = sum(hold_cash) # 手持ちの現金
            ratio_temp[i] = hold_stock_price / (hold_cash_price + hold_stock_price) # 現状の株式：現金比率
        end
    end

    return (hold_stock, hold_cash)
end

hold_stocks = zeros(size(spxl_30y))
hold_cash = zeros(size(spxl_30y))
for n in 1:data_num
    prices = spxl_30y[:,n]
    ratio = 0.7
    timestamp = spxl.timestamp
    fee = 0.45 / 100
    (stock, cash) = rebalance(prices, ratio, purchases, timestamp, fee)
    hold_stocks[:,n] = stock
    hold_cash[:,n] = cash
end

# 購入株式数を積算して保有株式数とする
stocks_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    stocks_spxl[:,n] = accumulate(+, hold_stocks[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_spxl = copy(stocks_spxl)
for c = 1:data_num
    assets_spxl[:,c] .*= spxl_30y[:,c]
end

# 現金も積算する
total_cash = zeros(size(spxl_30y))
for n in 1:data_num
    total_cash[:,n] = accumulate(+, hold_cash[:,n])
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_assets[r,:] = quantile(assets_spxl[r,:], [0.05, 0.5, 0.95])
end
quan_cash = zeros(data_len, 3)
for r in 1:data_len
    quan_cash[r,:] = quantile(total_cash[r,:], [0.05, 0.5, 0.95])
end

plot( (1:data_len) .* 30 / data_len, quan_assets[:,3], label="金融資産額" * labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,2], label="金融資産額" * labels[2], linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,1], label="金融資産額" * labels[1], linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,3], label="現金" * labels[3], linestyle=:dot, linecolor=:blue)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,2], label="現金" * labels[2], linestyle=:dot, linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,1], label="現金" * labels[1], linestyle=:dot, linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,3] .+ quan_cash[:,3], label="合計" * labels[3], linestyle=:dash, linecolor=:blue)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,2] .+ quan_cash[:,2], label="合計" * labels[2], linestyle=:dash, linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,1] .+ quan_cash[:,1], label="合計" * labels[1], linestyle=:dash, linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black)
savefig("SPXL_30y_5parcentile_assets_70_30リバランス時.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string((quan_assets[end, i] * (1 - 0.20315) + quan_cash[end,i]))[1:7] * "万円\n")
end
quan_assets[2513, 2] + quan_cash[2513, 2] # 10年時点の中央値

print("下位5%の最大下落率[%]:" * string(100 - minimum((quan_assets[:,1] .+ quan_cash[:,1]) ./ total_purchases) * 100)[1:4])

# -------------------------------------
# (5)15年かけて購入，15年かけて売却
# -------------------------------------
# 年月日から購入金額を計算
data_len_2 = Int(floor(data_len / 2))
purchases = zeros(data_len)
amount = (5300 / (15 * 12))
purchases[1] = amount # 初回も同額
# 15年目まで購入
for i = 2:data_len_2
    if (Dates.Month(spxl.timestamp[i - 1]) != Dates.Month(spxl.timestamp[i])) # 月が変わったタイミングで投資
        purchases[i] = amount
    end
end

# 購入金額を積算して，総投資額も計算しておく
total_purchases = accumulate(+, purchases)

# 購入株式数を購入金額から計算
buy_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    buy_spxl[:,n] = purchases ./ spxl_30y[:,n]
end

# 売却株式数(きっかり15年で同量の株を売却する)
sale_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    sale_spxl[data_len_2:data_len_2 * 2 - 1,n] = buy_spxl[1:data_len_2,n]
end

# 購入株式数と売却株式数を積算して保有株式数を計算する
total_holdings = copy(buy_spxl)
for c = 1:data_num
    total_buy = accumulate(+, buy_spxl[:,c]) # 購入株式数の積算
    total_sale = accumulate(+, sale_spxl[:,c]) # 売却株式数の積算
    total_holdings[:,c] = total_buy .- total_sale # 保有株式数
end
# 保有株式数にその時の価格を掛け算して，金融資産額を計算
holdings_price = copy(total_holdings)
for c = 1:data_num
    holdings_price[:,c] .*= spxl_30y[:,c]
end

# 売却株式数にその時の価格を掛け算して，売却金額を算出
sale_price_spxl = copy(sale_spxl)
for c = 1:data_num
    sale_price_spxl[:,c] .*= spxl_30y[:,c]
end

# 売却額を積算して現金保有額を計算
cash_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    cash_spxl[:,n] = accumulate(+, sale_price_spxl[:,n])
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_spxl_holdings = zeros(data_len, 3)
for r in 1:data_len
    quan_spxl_holdings[r,:] = quantile(holdings_price[r,:], [0.05, 0.5, 0.95])
end
quan_spxl_cash = zeros(data_len, 3)
for r in 1:data_len
    quan_spxl_cash[r,:] = quantile(cash_spxl[r,:], [0.05, 0.5, 0.95])
end

# 対数グラフだと0が表示できないので，値を微妙に補正
quan_spxl_holdings .+= 10
quan_spxl_cash .+= 10

plot( (1:data_len) .* 30 / data_len, quan_spxl_holdings[:,3], label="金融資産額" * labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_holdings[:,2], label="金融資産額" * labels[2], linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_holdings[:,1], label="金融資産額" * labels[1], linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_cash[:,3], label="売却合計額" * labels[3], linestyle=:dot, linecolor=:blue)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_cash[:,2], label="売却合計額" * labels[2], linestyle=:dot, linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_spxl_cash[:,1], label="売却合計額" * labels[1], linestyle=:dot, linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black)
savefig("SPXL_30y_5parcentile_assets_15年積立15年売却時.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string(quan_spxl_cash[end, i] * (1 - 0.20315) / total_purchases[end])[1:4] * "倍\n")
end
print("下位5%の最大下落率[%]:" * string(100 - minimum((quan_spxl_holdings[:,1] ./ total_purchases)[1:data_len_2]) * 100)[1:4])

# -------------------------------------
# (6)タイミング投資(直近高値より10%下落で購入)
# -------------------------------------
# 最初に200万円投資，残り300万あり．毎月10万＋ボーナス月は＋20万追加．
# 直近高値に対して10%を下回っていた場合，下落幅に対して以下の式で追加投資を毎日する．
# ((下落率-8)^1.1)+12)*0.2[万円]

# 年月日から投資への現金追加額を計算
purchases = zeros(data_len)
purchases[1] = 500 # 初回500万円
mid_flag = false
for i = 2:data_len
    if (Dates.Month(spxl.timestamp[i - 1]) != Dates.Month(spxl.timestamp[i]))
        purchases[i] = 10
        mid_flag = false
    elseif ( (Dates.Month(spxl.timestamp[i]) == Dates.Month(6) || Dates.Month(spxl.timestamp[i]) == Dates.Month(12)) &&
        Dates.Day(spxl.timestamp[i]) > Dates.Day(15) && mid_flag == false)
        purchases[i] = 20
        mid_flag = true
    end
end
# 購入金額を積算して，総投資額も計算しておく
total_purchases = accumulate(+, purchases)

# ナンピン投資を関数化
# price: 株価推移
# purchases: 追加投資現金
function additional(prices, purchases, timestamp, fee)
    # for debug
    # prices = spxl_30y[:,1]
    # timestamp = spxl.timestamp
    # fee = 0.45/100

    # 変数定義と初期化
    recent_high_price = prices[1] # 直近高値
    hold_stock = zeros(length(prices)) # 保有株式数推移
    hold_stock[1] = purchases[1] * 2.0 / 5.0 / prices[1]
    hold_cash = zeros(length(prices)) # 保有現金推移
    hold_cash[1] = purchases[1] * 3.0 / 5.0
    decrease_rate_total = zeros(length(prices)) # 下落幅

    # 2日目以降を計算
    for i in 2:length(prices)
        # 保有現金を更新
        hold_cash[i] = hold_cash[i - 1] + purchases[i]
        # 直近高値を更新=直近1ヶ月の最高値
        recent = 10
        if (i <= recent)
            recent_high_price = maximum(prices[1:i])
        else
            recent_high_price = maximum(prices[i - recent:i])
        end

        if (hold_cash[i] > 0)
            # 高値から10%を下回っていたら
            decrease_rate = 100.0 - prices[i] / recent_high_price * 100.0 # 下落率
            decrease_rate_total[i] = decrease_rate
            if ( decrease_rate > 10.0 )
                # 投資額を計算して，その分の株式を購入
                amount = (((decrease_rate - 8)^1.1) + 12) * 0.2
                if ( amount > hold_cash[i] ) # 保有現金が不足した場合，保有現金分しか買わない
                    amount = hold_cash[i] / (1 + fee)
                end
                hold_stock[i] = amount / prices[i]
                hold_cash[i] -= amount * (1 + fee)
            end
        end
    end

    return (hold_stock, hold_cash)
end

hold_stocks = zeros(size(spxl_30y))
hold_cash = zeros(size(spxl_30y))
for n in 1:data_num
    prices = spxl_30y[:,n]
    timestamp = spxl.timestamp
    fee = 0.45 / 100
    (stock, cash) = additional(prices, purchases, timestamp, fee)
    hold_stocks[:,n] = stock
    hold_cash[:,n] = cash
end

# 購入株式数を積算して保有株式数とする
stocks_spxl = zeros(size(spxl_30y))
for n in 1:data_num
    stocks_spxl[:,n] = accumulate(+, hold_stocks[:,n])
end

# 保有株式数にその時の価格を掛け算して，資産額を算出
assets_spxl = copy(stocks_spxl)
for c = 1:data_num
    assets_spxl[:,c] .*= spxl_30y[:,c]
end

# 5パーセンタイル，50パーセンタイル, 95パーセンタイルを計算してプロット
quan_assets = zeros(data_len, 3)
for r in 1:data_len
    quan_assets[r,:] = quantile(assets_spxl[r,:], [0.05, 0.5, 0.95])
end
quan_cash = zeros(data_len, 3)
for r in 1:data_len
    quan_cash[r,:] = quantile(hold_cash[r,:], [0.05, 0.5, 0.95])
end
# quan_cashが0だと対数グラフで描画できないので1万円足す
quan_cash .+= 1

#=
plot( (1:data_len) .* 30 / data_len, quan_assets[:,3], label="金融資産額" * labels[3], linecolor=:blue, xlabel="年", ylabel="資産額[万円]", yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,2], label="金融資産額" * labels[2], linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,1], label="金融資産額" * labels[1], linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,3], label="現金" * labels[3], linestyle=:dot, linecolor=:blue)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,2], label="現金" * labels[2], linestyle=:dot, linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_cash[:,1], label="現金" * labels[1], linestyle=:dot, linecolor=:red) =#
plot( (1:data_len) .* 30 / data_len, quan_assets[:,3] .+ quan_cash[:,3], label="合計" * labels[3], linecolor=:blue, yscale=:log10)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,2] .+ quan_cash[:,2], label="合計" * labels[2], linecolor=:green)
plot!( (1:data_len) .* 30 / data_len, quan_assets[:,1] .+ quan_cash[:,1], label="合計" * labels[1], linecolor=:red)
plot!( (1:data_len) .* 30 / data_len, total_purchases, label=labels[4], linecolor=:black)
savefig("SPXL_30y_5parcentile_assets_タイミング投資.png")

# 成績の計算
for i = 3:-1:1
    print(labels[i] * ": " * string((quan_assets[end, i] * (1 - 0.20315) + quan_cash[end,i]) / total_purchases[end])[1:4] * "倍\n")
end
(quan_assets[2513, 2] + quan_cash[2513, 2]) / total_purchases[2513]
for i = 3:-1:1
    print(labels[i] * ": " * string((quan_assets[end, i] * (1 - 0.20315) + quan_cash[end,i]))[1:7] * "万円\n")
end
quan_assets[2513, 2] + quan_cash[2513, 2]

print("下位5%の最大下落率[%]:" * string(100 - minimum((quan_assets[:,1] .+ quan_cash[:,1]) ./ total_purchases) * 100)[1:4])


text:
using KernelAbstractions
using AMDGPU
using ROCKernels
using CUDA
using CUDAKernels

# One day for Intel GPUs
# using oneAPI
# using OneKernels

using LinearAlgebra
using Printf
using PowerModels
using LazyArtifacts
using ExaTronKernels
using Test
const KA = KernelAbstractions
include("../example/admm/opfdata.jl")
include("../example/admm/environment.jl")
include("../example/admm/generator_kernel.jl")
include("../example/admm/eval_kernel.jl")
include("../example/admm/polar_kernel.jl")
include("../example/admm/bus_kernel.jl")
include("../example/admm/tron_kernel.jl")
include("../example/admm/acopf_admm_gpu.jl")

CASE = joinpath(artifact"ExaData", "ExaData", "matpower", "case9.m")

function one_level_admm(case::String, device::KA.Device)
    # NB: Need to run almost 2,000 iterations to reach convergence with this
    # set of parameters.
    env = admm_rect_gpu(
        case;
        verbose=1,
        iterlim=2000,
        rho_pq=400.0,
        rho_va=40000.0,
        scale=1e-4,
        use_polar=true,
        device=device
    )
    @test isa(env, AdmmEnv)

    model = env.model
    ngen = model.gen_mod.ngen

    par = env.params
    sol = env.solution

    # Check results
    pg = active_power_generation(env)
    qg = reactive_power_generation(env)

    @test sol.status == HAS_CONVERGED
    # Test with solution returned by PowerModels + Ipopt
    @test sol.objval ≈ 5296.6862 rtol=1e-4
    @test (pg |> Array) ≈ [0.897987, 1.34321, 0.941874] rtol=1e-4
    @test (qg |> Array) ≈ [0.1296564, 0.00031842, -0.226342] rtol=1e-2

    # Test restart API
    admm_restart!(env)
    return nothing
end

ext: jl


text:
module Points

export PointStyle, Is2DPoint, Is3DPoint
export get_x, get_y, get_z, set_x!, set_y!, set_z!
export distance

abstract type PointStyle end
struct Is2DPoint <: PointStyle end
struct Is3DPoint <: PointStyle end
struct IsNotPoint <: PointStyle end

PointStyle(::Type) = IsNotPoint()

# x accessors
get_x(p::T) where {T} = get_x(PointStyle(T), p)
get_x(::Union{Is2DPoint, Is3DPoint}, p) = p.x
get_x(::IsNotPoint, p) = throw(MethodError(get_x, p))

set_x!(p::T, val::Real) where {T} = set_x!(PointStyle(T), p, val)
set_x!(::Union{Is2DPoint, Is3DPoint}, p, val) = begin p.x = val; end
set_x!(::IsNotPoint, p, val) = throw(MethodError(set_x, [p, val]))

# y accessors
get_y(p::T) where {T} = get_y(PointStyle(T), p)
get_y(::Union{Is2DPoint, Is3DPoint}, p) = p.y
get_y(::IsNotPoint, p) = throw(MethodError(get_y, p))

set_y!(p::T, val::Real) where {T} = set_y!(PointStyle(T), p, val)
set_y!(::Union{Is2DPoint, Is3DPoint}, p, val) = begin p.y = val; end
set_y!(::IsNotPoint, p, val) = throw(MethodError(set_y, [p, val]))

# z accessors
get_z(p::T) where {T} = get_z(PointStyle(T), p)
get_z(::Is3DPoint, p) = p.z
get_z(::IsNotPoint, p) = throw(MethodError(get_z, p))

set_z!(p::T, val::Real) where {T} = set_z!(PointStyle(T), p, val)
set_z!(::Is3DPoint, p, val) = begin p.z = val; end
set_z!(::IsNotPoint, p, val) = throw(MethodError(set_z, [p, val]))

# distance function
distance(p1::T, p2::T) where {T} = distance(PointStyle(T), p1, p2)
distance(::Is2DPoint, p1, p2) = begin
    x1 = get_x(p1); y1 = get_y(p1)
    x2 = get_x(p2); y2 = get_y(p2)
    return sqrt((x1 - x2)^2 + (y1 - y2)^2)
end
distance(::Is3DPoint, p1, p2) = begin
    x1 = get_x(p1); y1 = get_y(p1); z1 = get_z(p1)
    x2 = get_x(p2); y2 = get_y(p2); z2 = get_z(p2)
    return sqrt((x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2)
end
distance(::IsNotPoint, p1, p2) = throw(MethodError(distance, [p1, p2]))

end


text:
## Exercise 4-9
## Write an appropriately general set of functions that can draw flowers as in Turtle flowers.
## https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#fig04-2
println("Ans: ")

println("End.")


text:
"""
Used for concise reporting of modeling results.
"""
mutable struct ReactionStatus
    already_present::Bool
    index::Int
    info::String
end


text:
using ModelCheck

# Check qqplot works
ModelCheck.plot_qqplot(model = modPoisson, df = df, y = "y")

typeof(plot_qqplot(model = modPoisson,
                           df = df,
                           y = "y")) == Plot

# Check residual vs fitted plot works
ModelCheck.plot_residual_fitted(model = modPoisson, df = df, y = "y")

# Check scale-location plot works
ModelCheck.plot_scale_location(model = modPoisson, df = df, y = "y")

ModelCheck.plot_model_check(model = modPoisson, df = df, y = "y")


text:
using LinearAlgebra
import SparseArrays: spdiagm

# remove linear time trend
function detrend(y::Vector{Float64})
	n = length(y)
	@assert n >= 4

	x = [ones(n) collect(1:n)]
	c=x\y
	return y-x*c
end

# Hodrick-Prescott Filter
# derived from: http://www.econforge.org/posts/2014/juil./28/cef2014-julia/
# output is (trend component, cyclical component)
function hp_filter(y::Vector{Float64}, lambda::Real;detrendfirst::Bool = false)
	n = length(y)
	@assert n >= 4

	if detrendfirst
		y = detrend(y)
	end
	diag2 = lambda*ones(n-2)
	diag1 = [ -2lambda; -4lambda*ones(n-3); -2lambda ]
	diag0 = [ 1+lambda; 1+5lambda; (1+6lambda)*ones(n-4); 1+5lambda; 1+lambda ]

	D = spdiagm(-2=>diag2, -1=>diag1, 0=>diag0, 1=>diag1, 2=>diag2)
	t=D\y
	t, y-t
end

lastnonzero(X) = (y=filter(!iszero,X); !isempty(y) ? last(y) : 0.0)
dropcols(df::AbstractDataFrame, cols) = df[setdiff(names(df), cols)]
yrqtrfun(dt) = string(year(dt),"Q",quarterofyear(dt))

function rollingfillfirst(fun::Function, data::AbstractArray{T1,1}, windowspan::Int64,filler::T2) where {T1,T2<:Union{Missing,Number}}
	TL = size(data,1)
	FL = windowspan - 1
	outtmp = allowmissing(zeros(TL))
	outtmp[1:FL] .= filler
	outtmp[(1+FL):end]=rolling(fun,data,windowspan)
	return outtmp
end

function rollingfillfirst(fun::Function, data1::AbstractArray{T1,1},data2::AbstractArray{T1,1},windowspan::Int64,filler::T2) where {T1,T2<:Union{Missing,Number}}
	TL = size(data1,1)
	TL2 = size(data2,1)
	@assert TL == TL2
	FL = windowspan - 1
	outtmp = allowmissing(zeros(TL))
	outtmp[1:FL] .= filler
	outtmp[(1+FL):end]=rolling(fun,data1,data2,windowspan)
	return outtmp
end

function tomissing!(d::AbstractDataFrame,x::Union{Integer,Symbol})
	allowmissing!(d,x)
	d[findall(isnan,d[x]),x] = missing
end

function alltomissing!(d::AbstractDataFrame)
	for nm in names(d)
		if eltype(d[nm]) <: Float64
			if any(isnan,d[nm])
				tomissing!(d,nm)
			end
		end
	end
end

function maxdrawdown(R)
	T = length(R)
	P = 100 .* cumprod(1 .+ R)
	DD = zeros(size(R))
	peak = -999.99
	for i = 1:T
		if P[i] > peak
			peak = P[i]
		end
		DD[i] = (P[i]-peak)/peak
	end
	maxDD = minimum(DD)
end


text:
using .JSort

struct Sorter

end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using Arpack
# Sparse positive definite matrix together with a Cholesky factorization object
struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T}
  dim::Int
  mat::S
  chol::CholTypeSparse
  PDSparseMat{T,S}(d::Int,m::AbstractSparseMatrix{T},c::CholTypeSparse) where {T,S} = new{T,S}(d,m,c) #add {T} to CholTypeSparse argument once #14076 is implemented
end

function PDSparseMat(mat::AbstractSparseMatrix,chol::CholTypeSparse)
    d = size(mat, 1)
    size(chol, 1) == d ||
        throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
    PDSparseMat{eltype(mat),typeof(mat)}(d, mat, chol)
end

PDSparseMat(mat::SparseMatrixCSC) = PDSparseMat(mat, cholesky(mat))
PDSparseMat(fac::CholTypeSparse) = PDSparseMat(sparse(fac) |> x -> x*x', fac)

### Conversion
Base.convert(::Type{PDSparseMat{T}}, a::PDSparseMat) where {T<:Real} = PDSparseMat(convert(SparseMatrixCSC{T}, a.mat))

### Basics

dim(a::PDSparseMat) = a.dim
Base.Matrix(a::PDSparseMat) = Matrix(a.mat)
LinearAlgebra.diag(a::PDSparseMat) = diag(a.mat)

### Arithmetics

# add `a * c` to a dense matrix `m` of the same size inplace.
function pdadd!(r::Matrix, a::Matrix, b::PDSparseMat, c)
    @check_argdims size(r) == size(a) == size(b)
    _addscal!(r, a, b.mat, c)
end

*(a::PDSparseMat, c::T) where {T<:Real} = PDSparseMat(a.mat * c)
*(a::PDSparseMat, x::StridedVecOrMat) = a.mat * x
\(a::PDSparseMat{T}, x::StridedVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076

### Algebra

Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat( a\eye(T,a.dim) )
LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol)
LinearAlgebra.eigmax(a::PDSparseMat{T}) where {T<:Real} = convert(T,eigs(convert(SparseMatrixCSC{Float64,Int},a.mat), which=:LM, nev=1, ritzvec=false)[1][1]) #to avoid type instability issues in eigs, see e.g., julia issue #13929
LinearAlgebra.eigmin(a::PDSparseMat{T}) where {T<:Real} = convert(T,eigs(convert(SparseMatrixCSC{Float64,Int},a.mat), which=:SM, nev=1, ritzvec=false)[1][1]) #to avoid type instability issues in eigs, see e.g., julia issue #13929

### whiten and unwhiten

function whiten!(r::StridedVecOrMat, a::PDSparseMat, x::StridedVecOrMat)
    r[:] = sparse(chol_lower(a.chol)) \ x
    return r
end

function unwhiten!(r::StridedVecOrMat, a::PDSparseMat, x::StridedVecOrMat)
    r[:] = sparse(chol_lower(a.chol)) * x
    return r
end

### quadratic forms

quad(a::PDSparseMat, x::StridedVector) = dot(x, a * x)
invquad(a::PDSparseMat, x::StridedVector) = dot(x, a \ x)

function quad!(r::AbstractArray, a::PDSparseMat, x::StridedMatrix)
    for i in 1:size(x, 2)
        r[i] = quad(a, x[:,i])
    end
    return r
end

function invquad!(r::AbstractArray, a::PDSparseMat, x::StridedMatrix)
    for i in 1:size(x, 2)
        r[i] = invquad(a, x[:,i])
    end
    return r
end

### tri products

function X_A_Xt(a::PDSparseMat, x::StridedMatrix)
    z = x * sparse(chol_lower(a.chol))
    z * transpose(z)
end

function Xt_A_X(a::PDSparseMat, x::StridedMatrix)
    z = transpose(x) * sparse(chol_lower(a.chol))
    z * transpose(z)
end

function X_invA_Xt(a::PDSparseMat, x::StridedMatrix)
    z = a.chol \ collect(transpose(x))
    x * z
end

function Xt_invA_X(a::PDSparseMat, x::StridedMatrix)
    z = a.chol \ x
    transpose(x) * z
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ripgrep_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ripgrep")
JLLWrappers.@generate_main_file("ripgrep", UUID("e10fc14b-37cd-5cbc-b289-ad01b12ebaad"))
end  # module ripgrep_jll


text:
# Everything is experimental! If you want to use it, you need to import it yourself.
# If you find problems, please open an issue.

using Optim
using Turing
using LinearAlgebra
using StatsBase

#      Run like
# using Turing
# @model height(heights) = begin
#     μ ~ Normal(178, 20)
#     σ ~ Uniform(0, 50)
#     heights .~ Normal(μ, σ)
# end
# m = height(d2.height)
# res = quap(m)

# Find the MAP via optimization and get the information matrix at that point.
# Look if your solution converged. Sometimes even solutions that didn't converge
# might be pretty good, on the other hand, just because the solver converged that
# doesn't mean you got the point you are looking for; this isn't even global
# optimization. In any case, trying other optimization methods might help.
function quap(model::Turing.Model, args...; kwargs...)
    opt = optimize(model, MAP(), args...; kwargs...)

    coef = opt.values.array
    var_cov_matrix = informationmatrix(opt)
    sym_var_cov_matrix = Symmetric(var_cov_matrix)  # lest MvNormal complains, loudly
    converged = Optim.converged(opt.optim_result)

    distr = if length(coef) == 1
        Normal(coef[1], √sym_var_cov_matrix[1])  # Normal expects stddev
    else
        MvNormal(coef, sym_var_cov_matrix)       # MvNormal expects variance matrix
    end

    params = StatsBase.params(model)
    params_tuple = tuple(Symbol.(params)...)

    (coef = NamedTuple{params_tuple}(coef), vcov = sym_var_cov_matrix, converged = converged,
        distr = distr, params = [params...])
end

function StatsBase.params(model::Turing.Model)
    nt = model |> Turing.VarInfo |> Turing.tonamedtuple
    p = String[]
    for (a, v) in nt
        arr = a[1]
        var = v[1]
        if length(arr) != 1
            append!(p, ["$var[$i]" for i in 1:length(arr)])
        else
            push!(p, var)
        end
    end
    return p
end


text:
function Space()
    Space(ccall(dlsym(libchipmunk, :cpSpaceNew), Ptr{Void}, ()))
end

function free(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceFree), Void, (Ptr{Void},), space.ptr)
end

function get_iterations(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetIterations), Int32, (Ptr{Void},), space.ptr)
end

function set_iterations(space::Space, iterations::Integer)
    ccall(dlsym(libchipmunk, :cpSpaceSetIterations), Void, (Ptr{Void}, Int32,), space.ptr, iterations)
end

function get_gravity(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetGravity), Vect, (Ptr{Void},), space.ptr)
end

function set_gravity(space::Space, gravity::Vect)
    ccall(dlsym(libchipmunk, :cpSpaceSetGravity), Void, (Ptr{Void}, Vect,), space.ptr, gravity)
end

function set_damping(space::Space, damping::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetDamping), Void, (Ptr{Void}, Cdouble,), space.ptr, damping)
end

function get_damping(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetDamping), Cdouble, (Ptr{Void},), space.ptr)
end

function get_idle_speed_threshold(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetIdleSpeedThreshold), Cdouble, (Ptr{Void},), space.ptr)
end

function set_idle_speed_threshold(space::Space, threshold::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetIdleSpeedThreshold), Void, (Ptr{Void}, Cdouble,), space.ptr, threshold)
end

function get_sleep_time_threshold(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetSleepTimeThreshold), Cdouble, (Ptr{Void},), space.ptr)
end

function set_sleep_time_threshold(space::Space, threshold::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetSleepTimeThreshold), Void, (Ptr{Void}, Cdouble,), space.ptr, threshold)
end

function get_collision_slop(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetCollisionSlop), Cdouble, (Ptr{Void},), space.ptr)
end

function set_collision_slop(space::Space, slop::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetCollisionSlop), Void, (Ptr{Void}, Cdouble,), space.ptr, slop)
end

function get_collision_bias(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetCollisionBias), Cdouble, (Ptr{Void},), space.ptr)
end

function set_collision_bias(space::Space, bias::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetCollisionBias), Void, (Ptr{Void}, Cdouble,), space.ptr, bias)
end

function get_collision_persistence(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetCollisionPersistence), Cdouble, (Ptr{Void},), space.ptr)
end

function set_collision_persistence(space::Space, persistence::Real)
    ccall(dlsym(libchipmunk, :cpSpaceSetCollisionPersistence), Void, (Ptr{Void}, Cdouble,), space.ptr, persistence)
end

function get_userdata(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetUserData), Ptr{Void}, (Ptr{Void},), space.ptr)
end

function set_userdata(space::Space, userdata::Ptr{Void})
    ccall(dlsym(libchipmunk, :cpSpaceSetUserData), Void, (Ptr{Void}, Ptr{Void},), space.ptr, userdata)
end

function get_current_timestep(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceGetCurrentTimeStep), Cdouble, (Ptr{Void},), space.ptr)
end

function is_locked(space::Space)
    ccall(dlsym(libchipmunk, :cpSpaceIsLocked), Bool, (Ptr{Void},), space.ptr)
end

function set_default_collision_handler(space)
    CollisionHandler(ccall(dlsym(libchipmunk, :cpSpaceAddDefaultCollisionHandler), Ptr{Void}, (Ptr{Void},), space.ptr))
end

function add_shape(space::Space, shape::Shape)
    shape_type = typeof(shape)
    shape_type(ccall(dlsym(libchipmunk, :cpSpaceAddShape), Ptr{Void}, (Ptr{Void}, Ptr{Void},), space.ptr, shape.ptr))
end

function add_body(space::Space, body::Body)
    Body(ccall(dlsym(libchipmunk, :cpSpaceAddBody), Ptr{Void}, (Ptr{Void}, Ptr{Void},), space.ptr, body.ptr))
end

function add_constraint(space::Space, constraint::Constraint)
    constraint_type = typeof(constraint)
    constraint_type(ccall(dlsym(libchipmunk, :cpSpaceAddConstraint), Ptr{Void}, (Ptr{Void}, Ptr{Void},), space.ptr, constraint.ptr))
end

function remove_shape(space::Space, shape::Shape)
    ccall(dlsym(libchipmunk, :cpSpaceRemoveShape), Void, (Ptr{Void}, Ptr{Void},), space.ptr, shape.ptr)
end

function remove_body(space::Space, body::Body)
    ccall(dlsym(libchipmunk, :cpSpaceRemoveBody), Void, (Ptr{Void}, Ptr{Void},), space.ptr, body.ptr)
end

function remove_constraint(space::Space, constraint::Constraint)
    ccall(dlsym(libchipmunk, :cpSpaceRemoveConstraint), Void, (Ptr{Void}, Ptr{Void},), space.ptr, constraint.ptr)
end

# function point_query(space::Space, point::Vect, max_distance::Cdouble, filter::ShapeFilter, out::PointQueryInfo)

# end

function use_spatial_hash(space::Space, dim::Real, count::Integer)
    ccall(dlsym(libchipmunk, :cpSpaceUseSpatialHash), Void, (Ptr{Void}, Cdouble, Cint,), space.ptr, dim, count)
end

function body_iterator_func(body_ptr::Ptr{Void}, data::Ptr{Void})
    callback = unsafe_pointer_to_objref(data)
    callback(Body(body_ptr))
    return nothing
end

function each_body(space::Space, callback::Function)
    c_callback = cfunction(body_iterator_func, Void, (Ptr{Void}, Ptr{Void},))
    callback_ptr = pointer_from_objref(callback)
    ccall(dlsym(libchipmunk, :cpSpaceEachBody), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), space.ptr, c_callback, callback_ptr)
end

function shape_iterator_func(shape_ptr::Ptr{Void}, data::Ptr{Void})
    callback = unsafe_pointer_to_objref(data)
    callback(Shape(body_ptr))
    return nothing
end

function each_shape(space::Space, callback::Function)
    c_callback = cfunction(shape_iterator_func, Void, (Ptr{Void}, Ptr{Void},))
    callback_ptr = pointer_from_objref(callback)
    ccall(dlsym(libchipmunk, :cpSpaceEachShape), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), space.ptr, c_callback, callback_ptr)
end

function constraint_iterator_func(constraint_ptr::Ptr{Void}, data::Ptr{Void})
    callback = unsafe_pointer_to_objref(data)
    callback(Shape(body_ptr))
    return nothing
end

function each_constraint(space::Space, callback::Function)
    c_callback = cfunction(constraint_iterator_func, Void, (Ptr{Void}, Ptr{Void},))
    callback_ptr = pointer_from_objref(callback)
    ccall(dlsym(libchipmunk, :cpSpaceEachConstraint), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void},), space.ptr, c_callback, callback_ptr)
end

function step(space::Space, dt::Real)
    ccall(dlsym(libchipmunk, :cpSpaceStep), Void, (Ptr{Void}, Cdouble,), space.ptr, dt)
end

export Space, get_iterations, set_iterations, get_gravity, set_gravity, set_damping, get_damping, add_body,
remove_body, step, add_shape, add_constraint, remove_constraint, remove_shape, is_locked, get_current_timestep,
get_userdata, set_userdata, set_collision_persistence, get_collision_persistence, set_collision_bias,
get_collision_bias, get_collision_slop, set_collision_slop, set_idle_speed_threshold, get_idle_speed_threshold,
set_sleep_time_threshold, get_sleep_time_threshold, use_spatial_hash, each_body, each_shape, each_constraint


text:
type Edge <: Container
	# Represents the edge between 2 vertices

	id::UUID.Uuid
	containedobject::Any
	attachedproperties::Dict{String,Any}

	graph::Any
	tail::Vertex
	head::Vertex
	typelabel::String

	partiallyloaded::Bool

	function Edge(tail::Vertex, head::Vertex, properties::Dict{String,Any})
		# Construct an edge from a pair of vertices with a set of property values

		e = new()
		e.id = UUID.v4()
		e.containedobject = UnspecifiedValue
		e.typelabel = "Unspecified"
		e.tail = tail
		e.head = head

		e.partiallyloaded = false

		e.attachedproperties = properties

		return e
	end

	function Edge(tail::Vertex, head::Vertex)
		# Construct an edge from a pair of vertices

		e = Edge(tail, head, Dict{String, Any}())

		return e
	end

	function Edge(typelabel::String, tail::Vertex, head::Vertex, properties::Dict{String,Any}=Dict{String,Any}())
		# Construct an edge from a typelabel, a pair of vertices, and a set of property values

		e = Edge(tail, head, properties)
		e.typelabel = typelabel

		return e
	end
end

function edgeforobject(object::Any, tail::Vertex, head::Vertex, properties::Dict{String,Any}=Dict{String,Any}())
	# Construct an edge for a specified object, and  a pair of vertices with a set of property values

	e = Edge(tail, head, properties)
	e.containedobject = object

	return e
end

function edgeforobject(typelabel::String, object::Any, tail::Vertex, head::Vertex, properties::Dict{String,Any}=Dict{String,Any}())
	# Construct an edge with a typelabel for a specified object, and  a pair of vertices with a set of property values

	e = Edge(typelabel, tail, head, properties)
	e.containedobject = object

	return e
end


text:
# adapted from https://github.com/sadit/SimilaritySearchDemos/blob/main/Glove/create-index-and-umap.jl
# relased under the MIT licensing
#
# Both SimilaritySearch.jl and UMAP.jl take advantage of multithreading parallelism
# run as:
# julia -t64 -L src/umap-embedding-with-common-voc.jl

using SimilaritySearch, UMAP, LinearAlgebra, JLD2, Embeddings, Glob, CSV, DataFrames

function create_index(valid_tokens, modelname)
    E = Embeddings.load_embeddings(FastText_Text, modelname)
    mask = [(token in valid_tokens) for token in E.vocab]
    X = Matrix(E.embeddings[:, mask])
    vocab = E.vocab[mask]
    for c in eachcol(X)
        normalize!(c)
    end
    @show modelname => (size(vocab), size(X))
    dist = NormalizedCosineDistance()
    index = SearchGraph(; dist, db=MatrixDatabase(X))
    index.neighborhood.reduce = SatNeighborhood()
    push!(index.callbacks, OptimizeParameters(; kind=ParetoRecall()))
    index!(index; parallel_block=512)
    optimize!(index, OptimizeParameters(; kind=MinRecall(), minrecall=0.9))
    index, vocab
end

function main(lang)
    edir = "data/$lang/embeddings"
    cclist = [first(split(basename(modelname), '.')) for modelname in glob("$edir/*.vec")]
    common = CSV.read("data/SpanishLang/common-tokens-per-region.tsv.gz", DataFrame, delim='\t')
    subset!(common, :n_regions => n -> n .> 10)
    valid_tokens = Set(common.token)

    for cc in reverse(cclist)
        k = 33
        umapfile = "$edir/umap+index-common-tokens-$cc.k=$k.jld2"
        embfile = "$edir/umap-embeddings-common-tokens-$cc.jld2"
        @info cc
        isfile(embfile) && continue
        U2 = if isfile(umapfile)
            load(umapfile, "U2")
        else
            index, vocab = create_index(valid_tokens, "$edir/$cc.vec")
            U2 = UMAP_(index, 2; n_neighbors=k, init=:random)
            jldsave(umapfile; U2, k, vocab)
            U2
        end

        U2 = UMAP_(index, 2; n_neighbors=k, init=:random)  # spectral layout is too slow for the input-data's size
        U3 = UMAP_(U2, 3; init=:random)  # reuses input data
        jldsave(embfile, e2=U2.embedding, e3=U3.embedding)
    end
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
struct QDevScalingSpec{T} <: SpatQSimSpec
    realization::Int
    qdev_scale::T
    prep_file::String

    function QDevScalingSpec(realization::Int,
                             qdev_scale::T,
                             prep_file::String = "prep.h5") where T
        new{T}(realization, qdev_scale, prep_file)
    end
end

sim_value(spec::QDevScalingSpec) = spec.qdev_scale

comm_qdev_scale(spec::QDevScalingSpec) = sim_value(spec)


text:
@userplot MarginalKDE

@recipe function f(kc::MarginalKDE; levels=10, clip=((-3.0, 3.0), (-3.0, 3.0)))
    x,y = kc.args

    x = vec(x)
    y = vec(y)

    m_x = median(x)
    m_y = median(y)

    dx_l = m_x - quantile(x, 0.16)
    dx_h = quantile(x, 0.84) - m_x

    dy_l = m_y - quantile(y, 0.16)
    dy_h = quantile(y, 0.84) - m_y

    xmin = m_x + clip[1][1]*dx_l
    xmax = m_x + clip[1][2]*dx_h

    ymin = m_y + clip[2][1]*dy_l
    ymax = m_y + clip[2][2]*dy_h

    k = KernelDensity.kde((x, y))
    kx = KernelDensity.kde(x)
    ky = KernelDensity.kde(y)

    ps = [pdf(k, xx, yy) for (xx, yy) in zip(x,y)]

    ls = []
    for p in range(1.0/levels, stop=1-1.0/levels, length=levels-1)
        push!(ls, quantile(ps, p))
    end

    legend --> false
    layout := @layout [topdensity           _
                       contour{0.9w,0.9h} rightdensity]

    @series begin
        seriestype := :contour
        levels := ls
        fill := false
        colorbar := false
        subplot := 2
        xlims := (xmin, xmax)
        ylims := (ymin, ymax)

        (collect(k.x), collect(k.y), k.density')
    end

    ticks := nothing
    xguide := ""
    yguide := ""

    @series begin
        seriestype := :density
        subplot := 1
        xlims := (xmin, xmax)
        ylims := (0, 1.1*maximum(kx.density))

        x
    end

    @series begin
        seriestype := :density
        subplot := 3
        orientation := :h
        xlims := (0, 1.1*maximum(ky.density))
        ylims := (ymin, ymax)

        y
    end
end


text:
#! /usr/bin/env julia

println("Start test 'tests/imageop.jl'")

include("../src/DCellC.jl")
using DCellC

ops = [ Id(),
        Soften(3),
        FlipV(),
        FlipH(),
        Jitter(),
        PixelNoise(0.1),
        StretchV(1.1),
        StretchH(1.5),
        parse(ImageOp, "(Soften(1), 0.5) * (Jitter(2), 1)")
      ]

println("Test imageop.jl: Creation of base image operations works")

ran = (Soften(3), 1) * (Soften(1), 2) * (Jitter(), 3)
pip = Soften(3) * (PixelNoise(0.1), 2) * ran

println("Test imageop.jl: Creation of pipelines and random operations works")

ops = [ops; [ran, pip]]

lmg = synthesize(256, 256, 1000)

for op in ops
  apply(op, lmg)
end

println("Test imageop.jl: Applying operations works")

println("Test imageop.jl completed")

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
module NormalSmoothingSplines

###### Inteface definition
export prepare_approximation, construct_approximation, approximate
export prepare_interpolation, construct_interpolation, interpolate
export evaluate, evaluate_at
export evaluate_gradient, evaluate_derivative
export NormalSpline, RK_H0, RK_H1, RK_H2

export get_epsilon, estimate_epsilon, get_cond, estimate_cond
export estimate_accuracy
######

using LinearAlgebra

abstract type ReproducingKernel end
abstract type ReproducingKernel_0 <: ReproducingKernel end
abstract type ReproducingKernel_1 <: ReproducingKernel_0 end
abstract type ReproducingKernel_2 <: ReproducingKernel_1 end

abstract type AbstractSpline end

"""
`struct NormalSpline{T, RK} <: AbstractSpline where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Define a structure containing full information of a normal spline
# Fields
- `_kernel`: a reproducing kernel spline was built with
- `_compression`: factor of transforming the original node locations into unit hypercube
- `_nodes`: transformed function value interpolation nodes
- `_nodes_b`: transformed function value approximation nodes
- `_values`: function values at interpolation nodes
- `_values_lb`: function lower bound values at approximation nodes
- `_values_ub`: function upper bound values at approximation nodes
- `_d_nodes`: transformed function directional derivative interpolation nodes
- `_d_nodes_b`: transformed function directional derivative approximation nodes
- `_es`: normalized derivative directions at interpolation nodes
- `_es_b`: normalized derivative directions at approximation nodes
- `_d_values`: function directional derivative values at interpolation nodes
- `_d_values_lb`: function lower bound directional derivative values at approximation nodes
- `_d_values_ub`: function upper bound directional derivative values at approximation nodes
- `_min_bound`: minimal bounds of the original node locations area
- `_gram`: Gram matrix of the problem
- `_chol`: Cholesky factorization of the Gram matrix
- `_mu`: spline coefficients
- `_active`: active inequality constraint numbers at solution
- `_cond`: estimation of the Gram matrix condition number
- `_ier`: An integer flag. If it is equal to 0, the optimal solution was found.
          If it is equal to 1, the approximate solution was found. QP algorithm iterations were stopped
          because of small spline norm change.
          If it is equal to 2, the approximate solution was found. QP algorithm iterations were stopped
          because of maximum allowed number of iterations was reached.
          If it is equal to -1, the solution was not calculated.
"""
struct NormalSpline{T, RK} <: AbstractSpline where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    _kernel::RK
    _compression::T
    _nodes::Union{Matrix{T}, Nothing}
    _nodes_b::Union{Matrix{T}, Nothing}
    _values::Union{Vector{T}, Nothing}
    _values_lb::Union{Vector{T}, Nothing}
    _values_ub::Union{Vector{T}, Nothing}
    _d_nodes::Union{Matrix{T}, Nothing}
    _d_nodes_b::Union{Matrix{T}, Nothing}
    _es::Union{Matrix{T}, Nothing}
    _es_b::Union{Matrix{T}, Nothing}
    _d_values::Union{Vector{T}, Nothing}
    _d_values_lb::Union{Vector{T}, Nothing}
    _d_values_ub::Union{Vector{T}, Nothing}
    _min_bound::Union{Vector{T}, Nothing}
    _gram::Union{Matrix{T}, Nothing}
    _chol::Union{Cholesky{T, Matrix{T}}, Nothing}
    _mu::Union{Vector{T}, Nothing}
    _active::Union{Vector{Int}, Nothing}
    _cond::T
    _ier::Int
end

include("./NormalInterpolatingSplines.jl")
include("./ReproducingKernels.jl")
include("./GramMatrix.jl")
include("./Utils.jl")
include("./QP.jl")
include("./Interpolate.jl")
include("./Approximate.jl")

##
#include("./examples/Main.jl")
##

"""
`prepare_approximation(nodes_b::Matrix{T}, kernel::RK = RK_H0())
                       where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the approximating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes_b`: function value approximation nodes.
          This should be an `n×n_1_b` matrix, where `n` is dimension of the sampled space and
          `n_1_b` is the number of function value approximation nodes. It means that each column in the matrix defines one node.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_approximation(nodes_b::Matrix{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(nodes_b, kernel)
     return spline
end

"""
`prepare_approximation(nodes::Matrix{T}, nodes_b::Matrix{T}, kernel::RK = RK_H0())
                       where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the approximating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space and
           `n_1` is the number of function value interpolation nodes. It means that each column in the matrix defines one node.
- `nodes_b`: function value approximation nodes.
          This should be an `n×n_1_b` matrix, where `n` is dimension of the sampled space and
          `n_1_b` is the number of function value approximation nodes. It means that each column in the matrix defines one node.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_approximation(nodes::Matrix{T},
                               nodes_b::Matrix{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(nodes, nodes_b, kernel)
     return spline
end

"""
`construct_approximation(spline::NormalSpline{T, RK},
                         values_lb::Vector{T}, values_ub::Vector{T},
                         maxiter::Int, ftol::T)
                         where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

construct the approximating normal spline by calculating its coefficients and
completely initializing the `NormalSpline` object.
# Arguments
- `spline`: the partly initialized `NormalSpline` object returned by `prepare` function.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object and the number of QP algorithm iterations done.
"""
function construct_approximation(spline::NormalSpline{T, RK},
                                 values_lb::Vector{T},
                                 values_ub::Vector{T},
                                 maxiter::Int,
                                 ftol::T = T(1.e-3)
                                ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    spline, nit_done = _construct_approximation(spline, values_lb, values_ub, maxiter, ftol)
    return spline, nit_done
end

"""
`construct_approximation(spline::NormalSpline{T, RK}, values::Vector{T},
                         values_lb::Vector{T}, values_ub::Vector{T}, maxiter::Int, ftol::T)
                         where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

construct the approximating normal spline by calculating its coefficients and
completely initializing the `NormalSpline` object.
# Arguments
- `spline`: the partly initialized `NormalSpline` object returned by `prepare` function.
- `values`: function values at interpolation nodes.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object and the number of QP algorithm iterations done.
"""
function construct_approximation(spline::NormalSpline{T, RK},
                                 values::Vector{T},
                                 values_lb::Vector{T},
                                 values_ub::Vector{T},
                                 maxiter::Int,
                                 ftol::T = T(1.e-3)
                                ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    spline, nit_done = _construct_approximation(spline, values, values_lb, values_ub, maxiter, ftol)
    return spline, nit_done
end

"""
`approximate(nodes::Matrix{T}, values::Vector{T}, nodes_b::Matrix{T}, values_lb::Vector{T}, values_ub::Vector{T},
             kernel::RK, maxiter::Int, ftol::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the approximating normal spline.
# Arguments
- `nodes_b`: function value approximation nodes.
        This should be an `n×n_1_b` matrix, where `n` is dimension of the sampled space and
        `n_1_b` is the number of function value approximation nodes. It means that each column in the matrix defines one node.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object.
"""
function approximate(nodes_b::Matrix{T},
                     values_lb::Vector{T},
                     values_ub::Vector{T},
                     kernel::RK,
                     maxiter::Int,
                     ftol::T = T(1.e-3)
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(nodes_b, kernel)
     spline, nit_done = _construct_approximation(spline, values_lb, values_ub, maxiter, ftol)
     return spline
end

"""
`approximate(nodes::Matrix{T}, values::Vector{T}, nodes_b::Matrix{T}, values_lb::Vector{T}, values_ub::Vector{T},
             kernel::RK, maxiter::Int, ftol::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the approximating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space
           and `n_1` is the number of function value nodes.
           It means that each column in the matrix defines one node.
- `nodes_b`: function value approximation nodes.
        This should be an `n×n_1_b` matrix, where `n` is dimension of the sampled space and
        `n_1_b` is the number of function value approximation nodes. It means that each column in the matrix defines one node.
- `values`: function values at interpolation nodes.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object.
"""
function approximate(nodes::Matrix{T},
                     values::Vector{T},
                     nodes_b::Matrix{T},
                     values_lb::Vector{T},
                     values_ub::Vector{T},
                     kernel::RK,
                     maxiter::Int,
                     ftol::T = T(1.e-3)
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(nodes, nodes_b, kernel)
     spline, nit_done = _construct_approximation(spline, values, values_lb, values_ub, maxiter, ftol)
     return spline
end

"""
`evaluate(spline::NormalSpline{T, RK}, points::Matrix{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate the spline values at `points` locations.

# Arguments
- `spline: constructed `NormalSpline` object.
- `points`: locations at which spline values are evaluating.
            This should be an `n×m` matrix, where `n` is dimension of the sampled space
            and `m` is the number of locations where spline values are evaluating.
            It means that each column in the matrix defines one location.

Return: `Vector{T}` of the spline values at the locations defined in `points`.
"""
function evaluate(spline::NormalSpline{T, RK}, points::Matrix{T}
                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return _evaluate(spline, points)
end

"""
`evaluate_at(spline::NormalSpline{T, RK}, point::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate the spline value at the `point` location.

# Arguments
- `spline`: constructed `NormalSpline` object.
- `point`: location at which spline value is evaluating.
           This should be a vector of size `n`, where `n` is dimension of the sampled space.

Return: spline value at the location defined in `point`.
"""
function evaluate_at(spline::NormalSpline{T, RK}, point::Vector{T}
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return _evaluate(spline, reshape(point, :, 1))[1]
end

"""
`evaluate_gradient(spline::NormalSpline{T, RK}, point::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate gradient of the spline at the location defined in `point`.

# Arguments
- `spline`: the `NormalSpline` object returned by `interpolate` or `construct` function.
- `point`: location at which gradient value is evaluating.
           This should be a vector of size `n`, where `n` is dimension of the sampled space.

Note: Gradient of spline built with reproducing kernel RK_H0 does not exist at the spline nodes.

Return: `Vector{T}` - gradient of the spline at the `point` location.
"""
function evaluate_gradient(spline::NormalSpline{T, RK},
                           point::Vector{T}
                          ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return _evaluate_gradient(spline, point)
end

########

"""
`get_epsilon(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Get the 'scaling parameter' of Bessel Potential space the spline was built in.
# Arguments
- `spline`: prepared `NormalSpline` object.

Return: `ε` - the 'scaling parameter'.
"""
function get_epsilon(spline::NormalSpline{T, RK}
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return spline._kernel.ε
end

"""
`estimate_epsilon(nodes::Matrix{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Get the estimation of the 'scaling parameter'2 of Bessel Potential space the spline being built in.
It coincides with the result returned by `get_epsilon` function.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space
           and `n_1` is the number of function value nodes.
           It means that each column in the matrix defines one node.
- `kernel`: reproducing kernel of Bessel potential space the normal spline will be constructed in.
           It must be a struct object of the following type:
             `RK_H0` if the spline is constructing as a continuous function,
             `RK_H1` if the spline is constructing as a differentiable function,
             `RK_H2` if the spline is constructing as a twice differentiable function.
Return: estimation of `ε`.
"""
function estimate_epsilon(nodes::Matrix{T},
                          kernel::RK = RK_H0()
                         ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    ε = _estimate_epsilon(nodes, kernel)
    return ε
end

"""
`estimate_epsilon(nodes::Matrix{T}, d_nodes::Matrix{T}, kernel::RK = RK_H1()) where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Get an the estimation of the 'scaling parameter' of Bessel Potential space the spline being built in.
It coincides with the result returned by `get_epsilon` function.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space
           and `n_1` is the number of function value nodes.
           It means that each column in the matrix defines one node.
- `d_nodes`: function directional derivative nodes.
           This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
           `n_2` is the number of function directional derivative nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline will be constructed in.
            It must be a struct object of the following type:
            `RK_H1` if the spline is constructing as a differentiable function,
            `RK_H2` if the spline is constructing as a twice differentiable function.

Return: estimation of `ε`.
"""
function estimate_epsilon(nodes::Matrix{T},
                          d_nodes::Matrix{T},
                          kernel::RK = RK_H1()
                         ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    ε = _estimate_epsilon(nodes, d_nodes, kernel)
    return ε
end

"""
`estimate_cond(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Get an estimation of the Gram matrix condition number. It needs the `spline` object is prepared and requires O(N^2) operations.
(C. Brás, W. Hager, J. Júdice, An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP Vol.20, No.3, 2012.)
# Arguments
- `spline`: prepared `NormalSpline` object.

Return: an estimation of the Gram matrix condition number.
"""
function estimate_cond(spline::NormalSpline{T, RK}
                      ) where {T <: AbstractFloat, RK <: ReproducingKernel}
    return spline._cond
end

"""
`estimate_accuracy(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Assess accuracy of interpolation results by analyzing residuals.
# Arguments
- `spline`: constructed `NormalSpline` object.

Return: estimation of the number of significant digits in the interpolation result.
"""
function estimate_accuracy(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return _estimate_accuracy(spline)
end

############################## One-dimensional case

"""
`prepare_approximation(nodes_b::Vector{T}, kernel::RK = RK_H0())
                       where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the 1D approximating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes_b`: function value approximation nodes.
             This should be an `n_1_b` vector where `n_1_b` is the number of function value approximation nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_approximation(nodes_b::Vector{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}

     spline = _prepare_approximation(Matrix(nodes_b'), kernel)
     return spline
end

"""
`prepare_approximation(nodes::Vector{T}, nodes_b::Vector{T}, kernel::RK = RK_H0())
                       where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the 1D approximating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n_1` vector where `n_1` is the number of function value interpolation nodes.
- `nodes_b`: function value approximation nodes.
           This should be an `n_1_b` vector where `n_1_b` is the number of function value approximation nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_approximation(nodes::Vector{T},
                               nodes_b::Vector{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(Matrix(nodes'), Matrix(nodes_b'), kernel)
     return spline
end

"""
`approximate(values::Vector{T}, nodes_b::Vector{T}, values_lb::Vector{T}, values_ub::Vector{T},
             kernel::RK, maxiter::Int, ftol::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the 1D approximating normal spline.
# Arguments
- `nodes_b`: function value approximation nodes.
          This should be an `n_1_b` vector where `n_1_b` is the number of function value approximation nodes.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object.
"""
function approximate(nodes_b::Vector{T},
                     values_lb::Vector{T},
                     values_ub::Vector{T},
                     kernel::RK,
                     maxiter::Int,
                     ftol::T = T(1.e-3)
               ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(Matrix(nodes_b'), kernel)
     spline, nit_done = _construct_approximation(spline, values_lb, values_ub, maxiter, ftol)
     return spline
end

"""
`approximate(nodes::Vector{T}, values::Vector{T}, nodes_b::Vector{T}, values_lb::Vector{T}, values_ub::Vector{T},
             kernel::RK, maxiter::Int, ftol::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the 1D approximating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
          This should be an `n_1` vector where `n_1` is the number of function value interpolation nodes.
- `values`: function values at interpolation nodes.
- `nodes_b`: function value approximation nodes.
          This should be an `n_1_b` vector where `n_1_b` is the number of function value approximation nodes.
- `values_lb`: function lower bound values at approximation nodes
- `values_ub`: function upper bound values at approximation nodes
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.
- `maxiter`: Maximum allowed number of iterations.
- `ftol`: convergence tolerance. The iteration stops when relative spline norm change is smaller than ftol.

Return: constructed `NormalSpline` object.
"""
function approximate(nodes::Vector{T},
                     values::Vector{T},
                     nodes_b::Vector{T},
                     values_lb::Vector{T},
                     values_ub::Vector{T},
                     kernel::RK,
                     maxiter::Int,
                     ftol::T = T(1.e-3)
               ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_approximation(Matrix(nodes'), Matrix(nodes_b'), kernel)
     spline, nit_done = _construct_approximation(spline, values, values_lb, values_ub, maxiter, ftol)
     return spline
end

"""
`evaluate(spline::NormalSpline{T, RK}, points::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate the 1D spline values at the `points` locations.

# Arguments
- `spline`: constructed `NormalSpline` object.
- `points`: locations at which spline values are evaluating.
            This should be a vector of size `m` where `m` is the number of evaluating points.

Return: spline value at the `point` location.
"""
function evaluate(spline::NormalSpline{T, RK}, points::Vector{T}
                               ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    return _evaluate(spline, Matrix(points'))
end

"""
`evaluate_at(spline::NormalSpline{T, RK}, point::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate the 1D spline value at the `point` location.

# Arguments
- `spline`: constructed `NormalSpline` object.
- `point`: location at which spline value is evaluating.

Return: spline value at the `point` location.
"""
function evaluate_at(spline::NormalSpline{T, RK},
                                       point::T
                                      ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    v_points = Vector{T}(undef, 1)
    v_points[1] = point
    return _evaluate(spline, Matrix(v_points'))[1]
end

"""
`evaluate_derivative(spline::NormalSpline{T, RK}, point::T) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Evaluate the 1D spline derivative at the `point` location.

# Arguments
- `spline`: the `NormalSpline` object returned by `interpolate` or `construct` function.
- `point`: location at which spline derivative is evaluating.

Note: Derivative of spline built with reproducing kernel RK_H0 does not exist at the spline nodes.

Return: spline derivative value at the `point` location.
"""
function evaluate_derivative(spline::NormalSpline{T, RK},
                             point::T
                            ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    v_points = Vector{T}(undef, 1)
    v_points[1] = point
    return _evaluate_gradient(spline, v_points)[1]
end

"""
`estimate_epsilon(nodes::Vector{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Get an the estimation of the 'scaling parameter' of Bessel Potential space the spline being built in.
It coincides with the result returned by `get_epsilon` function.
# Arguments
- `nodes`: function value interpolation nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: estimation of `ε`.
"""
function estimate_epsilon(nodes::Vector{T},
                          kernel::RK = RK_H0()
                         ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    ε = _estimate_epsilon(Matrix(nodes'), kernel)
    return ε
end

"""
`estimate_epsilon(nodes::Vector{T}, d_nodes::Vector{T}, kernel::RK = RK_H1()) where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Get an the estimation of the 'scaling parameter' of Bessel Potential space the spline being built in.
It coincides with the result returned by `get_epsilon` function.
# Arguments
- `nodes`: function value interpolation nodes.
- `d_nodes`: function derivative nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: estimation of `ε`.
"""
function estimate_epsilon(nodes::Vector{T},
                          d_nodes::Vector{T},
                          kernel::RK = RK_H1()
                         ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    ε = _estimate_epsilon(Matrix(nodes'), Matrix(d_nodes'), kernel)
    return ε
end

"""
`get_cond(nodes::Matrix{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Get a value of the Gram matrix spectral condition number. It is obtained by means of the matrix SVD decomposition.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space and
           `n_1` is the number of function value nodes. It means that each column in the matrix defines one node.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: a value of the Gram matrix spectral condition number.
"""
function get_cond(nodes::Matrix{T}, kernel::RK) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
  return _get_cond(nodes, kernel)
end

"""
`get_cond(nodes::Matrix{T}, d_nodes::Matrix{T}, es::Matrix{T}, kernel::RK = RK_H1()) where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Get a value of the Gram matrix spectral condition number. It is obtained by means of the matrix SVD decomposition.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space and
           `n_1` is the number of function value nodes.
            It means that each column in the matrix defines one node.
- `d_nodes`: function directional derivatives nodes.
             This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
             `n_2` is the number of function directional derivative nodes.
- `es`: Directions of the function directional derivatives.
        This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
        `n_2` is the number of function directional derivative nodes.
        It means that each column in the matrix defines one direction of the function directional derivative.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: value of the Gram matrix spectral condition number.
"""
function get_cond(nodes::Matrix{T}, d_nodes::Matrix{T}, es::Matrix{T}, kernel::RK = RK_H1()
                 ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    return _get_cond(nodes, d_nodes, es, kernel)
end

end # module


text:
"""
    PeriodicSelector(period::DatePeriod, stride::DatePeriod=Day(1), offset::DatePeriod=Day(0))

Assign holdout dates by taking a set of size `stride` once per `period`.
The offset is relative to _Monday 1st Jan 1900_, and controls when the selected section starts.

For example, `PeriodicSelector(Week(1), Day(2), Day(1))` will select 2 days per week.
With this selected periods offset by 1 day from 1st Jan 1900.
I.e. if applied to the first two weeks of the year 1900,
it would select 2nd, 3rd, 9th and 8th of Jan 1900.

Note: this cannot be actually used to select days earlier than `offset` after
1st Jan 1900.
"""
struct PeriodicSelector <: DateSelector
    period::DatePeriod
    stride::DatePeriod
    offset::DatePeriod

    function PeriodicSelector(period, stride=Day(1), offset=Day(0))
        period ≥ Day(2) || throw(DomainError(period, "period must be at least 2 Days."))
        stride ≥ Day(1) || throw(DomainError(stride, "stride must be at least 1 Day."))

        if Day(stride) > Day(period)
            throw(ArgumentError(
                "Cannot take a $stride stride within a $period period."
            ))
        end

        return new(period, stride, offset)
    end
end

function Iterators.partition(dates::StepRange{Date, Day}, s::PeriodicSelector)
    initial_time = _initial_date(s, dates)
    sd, ed = extrema(dates)

    #NOTE: you might be thinking that this process that actually checks all dates starting
    # from year 1900 is too slow and we should do something smart with modulo arithmetic
    # but for current decades this takes thousands of a second and even for year 9000
    # is still is well under 1/4 second. so keeping it simple

    holdout_dates = Date[]
    curr_window = initial_time:step(dates):(initial_time + s.stride - step(dates))
    while first(curr_window) <= ed
        # optimization: only creating holdout window if intersect not empty
        if last(curr_window) >= sd
            curr_holdout_window = curr_window ∩ dates
            append!(holdout_dates, curr_holdout_window)
        end
        curr_window = curr_window .+ s.period
    end

    return _getdatesets(dates, holdout_dates)
end


text:
using Base.Test

using Blackbody

custom_handler( r::Test.Success ) = println( "PASS $(r.expr)" )
custom_handler( r::Test.Failure ) = println( "FAIL $(r.expr)" )
custom_handler( r::Test.Error   ) = rethrow( r )

Test.with_handler( custom_handler ) do
    for i = 1:20
        temperature = 1000.0 * i
        @test 1 == 1
        @test_approx_eq_eps quadgk( wl -> planck( wl, temperature, PerAngstrom ), 0.0, 1.0e9 )[ 1 ] / stefan_boltzmann( temperature ) 1.0 1.0e-7
    end
end


text:
using FunctionalCollections: PersistentVector

#######################################
# trace for vector of generator calls #
#######################################

"""

U is the type of the subtrace, R is the return value type for the kernel
"""
struct VectorTrace{T,U}
    subtraces::PersistentVector{U}
    call::CallRecord{PersistentVector{T}}
    is_empty::Bool
end

function VectorTrace{T,U}() where {T,U}
    VectorTrace{T,U}(PersistentVector{U}(), PersistentVector{T}(), true)
end

function get_subtrace(trace::VectorTrace{T,U}, i::Int) where {T,U}
    get(trace.subtraces, i)
end

# TODO need to manage is_empty

# trace API

get_call_record(trace::VectorTrace) = trace.call
has_choices(trace::VectorTrace) = !trace.is_empty
get_assignment(trace::VectorTrace) = VectorTraceAssignment(trace)

struct VectorTraceAssignment <: Assignment
    trace::VectorTrace
end

Base.isempty(assignment::VectorTraceAssignment) = assignment.trace.is_empty
get_address_schema(::Type{VectorTraceAssignment}) = VectorAddressSchema()
has_internal_node(assignment::VectorTraceAssignment, addr) = false

function has_internal_node(assignment::VectorTraceAssignment, addr::Int)
    n = length(assignment.trace.subtraces)
    addr >= 1 && addr <= n
end

function has_internal_node(assignment::VectorTraceAssignment, addr::Pair)
    (first, rest) = addr
    sub_assignment = get_assignment(assignment.trace.subtraces[first])
    has_internal_node(sub_assignment, rest)
end

function get_internal_node(assignment::VectorTraceAssignment, addr::Int)
    get_assignment(assignment.trace.subtraces[addr])
end

function get_internal_node(assignment::VectorTraceAssignment, addr::Pair)
    (first, rest) = addr
    sub_assignment = get_assignment(assignment.trace.subtraces[first])
    get_internal_node(sub_assignment, rest)
end

has_leaf_node(assignment::VectorTraceAssignment, addr) = false

function has_leaf_node(assignment::VectorTraceAssignment, addr::Pair)
    (first, rest) = addr
    sub_assignment = get_assignment(assignment.trace.subtraces[first])
    has_leaf_node(sub_assignment, rest)
end

function get_leaf_node(assignment::VectorTraceAssignment, addr::Pair)
    (first, rest) = addr
    sub_assignment = get_assignment(assignment.trace.subtraces[first])
    get_leaf_node(sub_assignment, rest)
end

function get_internal_nodes(assignment::VectorTraceAssignment)
    [(i, get_assignment(assignment.trace.subtraces[i])) for i=1:length(assignment.trace.subtraces)]
end

get_leaf_nodes(assignment::VectorTraceAssignment) = []

##########################################
# trace for vector of distribution calls #
##########################################

struct VectorDistTrace{T}
    values::PersistentVector{T}
    call::CallRecord{PersistentVector{T}}
end

function VectorDistTrace{T}() where {T}
    VectorDistTrace{T}(PersistentVector{T}())
end

# trace API

get_call_record(trace::VectorDistTrace) = trace.call
has_choices(trace::VectorDistTrace) = length(trace.values) > 0
get_assignment(trace::VectorDistTrace) = VectorDistTraceAssignment(trace)

struct VectorDistTraceAssignment <: Assignment
    trace::VectorDistTrace
end

Base.isempty(assignment::VectorDistTraceAssignment) = !has_choices(assignment.trace)
get_address_schema(::Type{VectorDistTraceAssignment}) = VectorAddressSchema()
has_internal_node(assignment::VectorDistTraceAssignment, addr) = false
has_leaf_node(assignment::VectorDistTraceAssignment, addr) = false

function get_leaf_node(assignment::VectorDistTraceAssignment, addr::Int)
    assignment.trace.values[addr]
end

get_internal_nodes(assignment::VectorDistTraceAssignment) = ()
get_leaf_nodes(assignment::VectorDistTraceAssignment) = assignment.values

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen