text:
module TokyoConfirmed4
using Logging
using CSV
using DataFrames
using Dates
using Statistics: mean, stdm
using LsqFit: curve_fit
#using SingularSpectrumAnalysis: analyze
using RCall
using Printf
using Plots
using Plots.PlotMeasures

export date_to_value, value_to_date

@info "module", @__MODULE__

function has_name(df::DataFrame, name::String)
    return !isnothing(findfirst(isequal(name), names(df)))
end

@info "--------"
@info "作業ディレクトリ"
@info pwd()

@info "--------"
@info "基準日・開始日"
const epoch_date = Date("2021-01-01")
@info "epoch_date", epoch_date
#
date_to_value(d::Date) = Dates.value(d - epoch_date)
value_to_date(v)::Date = epoch_date + Day(v)
#
const start_date = Date("2021-01-01")
@info "start_date", start_date

@info "========"
@info "データ読み込み・変換"

function load_df(csv_filepath)::DataFrame
    @info csv_filepath

    @info "--------"
    @info "CSV 読込み"
    csv_df = CSV.read(csv_filepath, DataFrame; missingstring="NA")
    @info size(csv_df)

    @info "--------"
    @info "データフレーム作成"
    df = DataFrame()
    @info "--------"
    @info "日付 (date_start, date_end, days, date_error, date_value, date_mid)"
    @assert has_name(csv_df, "date")
    dates = csv_df.date
    df.date_start  = dates
    df.date_end    = dates
    df.days        = date_to_value.(df.date_end) .+ 1.0 .- date_to_value.(df.date_start)
    df.date_error  = df.days ./ 2.0
    df.date_value  = date_to_value.(df.date_start) .+ df.date_error
    df.date_mid    = value_to_date.(round.(df.date_value .- 0.5))
    @assert all(df.date_start .≤ df.date_mid .≤ df.date_end)
    @info "df.date_start[end]", df.date_start[end]
    @info "df.date_mid[end]", df.date_mid[end]
    @info "df.date_end[end]", df.date_end[end]
    @info size(df)
    n_rows = nrow(df)
    @info "--------"
    @info "曜日・祝日・休日 (holiday, holiday_next)"
    df.day_of_week = dayofweek.(df.date_mid)
    has_name(csv_df, "h")
    holidays = Vector{Bool}(undef, nrow(df))
    for i in 1:n_rows
        w = df.day_of_week[i]
        h = csv_df.h[i]
        holidays[i] = ((h == "祝") || (h == "休")) && (w != 6) && (w != 7)
    end
    df.holiday = holidays
    df.holiday_next = ([false; df.holiday])[1:end-1]
    @info size(df)
    @info "--------"
    @info "感染確認者数のコピー (count)"
    @assert has_name(csv_df, "count")
    df.count = csv_df.count
    @info size(df)
    @info "--------"
    @info "開始日以降への制限"
    df2 = filter((row) -> (row.date_start ≥ start_date), df)
    @info size(df2)
    return df2
end

const tokyo_confirmed_csv_filename = "東京感染確認者4.csv"
confirmed_df = load_df(tokyo_confirmed_csv_filename)

@info "========"
@info "時系列解析"

const ci95s = 1.95996398454005

const ssa_window = 7
@info "ssa_window", ssa_window
const qc_period = 21
@info "qc_period", qc_period

linreg(x, β) = β[1] .+ β[2] .* x
linreg_inv(y, β) = (y .- β[1]) ./ β[2]
qc_reg(x, α) = α[1] .+ x .* (α[2] .+ x .* α[3])

function ssanalysis!(df::DataFrame)
    n_rows = nrow(df)
    @info "--------"
    @info "予備特異スペクトル解析"
    @info "トレンド (xs0), 周期成分 (ss0), 残差 (rs0) への分解. ps0 は ss0 + rs0"
    cs = df.count
    #=
    (xs0, ss0_array) = analyze(cs, ssa_window; robust=false)
    @info "****", size(ss0_array)
    ps0 = cs .- xs0     # pecular = count - trend
    ss0 = zeros(Float64, n_rows)
    for ss0_col in eachcol(ss0_array)
        ss0 .+= ss0_col
    end
    rs0 = ps0 .- ss0    # residual = pecular - seasonal
    =#
    cs_stl = rcopy(R"""
    cs.ts <- ts($cs, frequency=$ssa_window, start=1)
    cs.stl <- stl(cs.ts, s.window=$ssa_window)
    """)
    cs_rts = cs_stl[:time_series]
    @info size(cs_rts)
    ps0 = cs_rts[:,1] .+ cs_rts[:,3]    # seasonal + residual
    @info "--------"
    @info "超過量補正 (ladjust1)"
    @info "周期的変動の大きさが count に依存することの補正"
    ps0_abs = 2.0 * abs.(ps0)
    β0 = [50.0, 0.1]
    p1fit = curve_fit(linreg, cs, ps0_abs, β0)
    @info p1fit.param
    @assert p1fit.param[1] ≥ 0.0
    lcs1 = log.(linreg(cs, p1fit.param))
    df.ladjust1 = lcs1
    @info size(df)
    #=
    @info "--------"
    @info "2次係数補正 (ladjust2)"
    @info "直近の2次係数が0となるよう補正"
    α0 = [0.0, 0.0, 0.0]
    qc_x = df.date_value[end-qc_period:end]
    qc_y = lcs1[end-qc_period:end]
    p2fit = curve_fit(qc_reg, qc_x, qc_y, α0)
    qc_α2 = p2fit.param[3]
    lcs2 = lcs1 .- qc_α2 .* (df.date_value).^2
    @info "qc_α2", qc_α2
    df.ladjust2 = lcs2
    @info size(df)
    =#
    @info size(df)
    @info "--------"
    @info "特異スペクトル解析 (ltrend1, lseasonal)"
    #=
    (lxs, lss_array) = analyze(lcs2, ssa_window; robust=false)
    @info "****", size(lss_array)
    lps = lcs2 .- lxs
    lss = zeros(Float64, n_rows)
    for lss_col in eachcol(lss_array)
        lss .+= lss_col
    end
    lrs = lps .- lss
    df.ltrend    = lxs
    df.lpeculiar = lps
    df.lseasonal = lss
    df.lresidual = lrs
    =#
    lcs_stl = rcopy(R"""
    lcs.ts <- ts($(df.ladjust1), frequency=$ssa_window, start=1)
    lcs.stl <- stl(lcs.ts, s.window=$ssa_window, l.degree=1, t.degree=1, s.degree=1, robust=T)
    plot(lcs.stl)
    lcs.stl
    """)
    lcs_rts = lcs_stl[:time_series]
    @info size(lcs_rts)
    df.lseasonal = lcs_rts[:,1]
    df.ltrend1   = lcs_rts[:,2]
    df.lresidual = lcs_rts[:,3]
    @info size(df)
    @info "--------"
    @info "残差誤差 (ltrend1_cil, ltrend1_ciu)"
    lrs = filter((row) -> (!row.holiday_next), df).lresidual
    lrm = mean(lrs)
    lre = ci95s * stdm(lrs, lrm)
    @info lrm, lre
    df.ltrend1_cil = df.ltrend1 .- (lre + lrm)
    df.ltrend1_ciu = df.ltrend1 .+ (lre - lrm)
    @info size(df)
    @info "--------"
    @info "補正トレンド (trend, trend_cil, trend_ciu)"
    #=
    df.ltrend1     = df.ltrend2     .+ qc_α2 .* (df.date_value).^2
    df.ltrend1_cil = df.ltrend2_cil .+ qc_α2 .* (df.date_value).^2
    df.ltrend1_ciu = df.ltrend2_ciu .+ qc_α2 .* (df.date_value).^2
    =#
    df.trend     = linreg_inv(exp.(df.ltrend1),     p1fit.param)
    df.trend_cil = linreg_inv(exp.(df.ltrend1_cil), p1fit.param)
    df.trend_ciu = linreg_inv(exp.(df.ltrend1_ciu), p1fit.param)
    for i in 1:n_rows
        @assert df.trend[i] ≥ 0.0
    end
    @info size(df)
    @info "--------"
    @info "拡大係数 (lgrowth, lgrowth_cil, lgrowth_ciu)"
    xs     = df.trend
    xs_cil = df.trend_cil
    xs_ciu = df.trend_ciu
    ξs     = log.(xs)
    ξs_cil = log.(xs_cil)
    ξs_ciu = log.(xs_ciu)
    ξds     = zeros(Float64, n_rows)
    ξds_cil = zeros(Float64, n_rows)
    ξds_ciu = zeros(Float64, n_rows)
    for i in 1:n_rows
        if     i == 1
            ξd = ξs[i+1] - ξs[i]
            ξd_cil = ξd - sqrt((ξs[i+1] - ξs_cil[i+1])^2 + (ξs_ciu[i] - ξs[i])^2)
            ξd_ciu = ξd + sqrt((ξs_ciu[i+1] - ξs[i+1])^2 + (ξs[i] - ξs_cil[i])^2)
        elseif i == n_rows
            ξd = ξs[i] - ξs[i-1]
            ξd_cil = ξd - sqrt((ξs[i] - ξs_cil[i])^2 + (ξs_ciu[i-1] - ξs[i-1])^2)
            ξd_ciu = ξd + sqrt((ξs_ciu[i] - ξs[i])^2 + (ξs[i-1] - ξs_cil[i-1])^2)
        else
            ξd = (ξs[i+1] - ξs[i-1]) / 2.0
            ξd_cil = ξd - sqrt((ξs[i+1] - ξs_cil[i+1])^2 + (ξs_ciu[i-1] - ξs[i-1])^2) / 2.0
            ξd_ciu = ξd + sqrt((ξs_ciu[i+1] - ξs[i+1])^2 + (ξs[i-1] - ξs_cil[i-1])^2) / 2.0
        end
        ξds[i]     = ξd
        ξds_cil[i] = ξd_cil
        ξds_ciu[i] = ξd_ciu
    end
    df.lgrowth     = ξds
    df.lgrowth_cil = ξds_cil
    df.lgrowth_ciu = ξds_ciu
    @info size(df)
end

ssanalysis!(confirmed_df)

@info "========"
@info "曜日変動"
function day_of_week_var(confirmed_df::DataFrame)
    wdf = DataFrame()
    lrs = skipmissing([((df.holiday_next[i]) ? missing : df.lresidual[i]) for i in 1:n_rows])
end

@info "========"
@info "株別推定"
@info "--------"
@info "尤度データ取り込み"
const likelihoods_csv_filepath = "tokyo_variant_analysis4_likelihoods.csv"
@info likelihoods_csv_filepath
likelihoods_df = CSV.read(likelihoods_csv_filepath, DataFrame; missingstring="NA")
@info size(likelihoods_df)
@info "--------"
@info "最大尤度"
likelihood_max = likelihoods_df[1,:]
likelihood_max_α = likelihood_max.α
likelihood_max_β = likelihood_max.β
@info likelihood_max
#
δ_logit(t; α, β) = α + β * t
δ_odds(t; α, β) = exp(δ_logit(t; α=α, β=β))
δ_prop(t; α, β) = 1.0 / (1.0 + 1.0 / δ_odds(t; α=α, β=β))
δ_prop_maxL(t) = δ_prop(t; α=likelihood_max_α, β=likelihood_max_β)
# δ_prop_maxL(t) は、時刻 t における最尤線が示す δ 株の割合
#
const likelihood_lags = [0, 7, 14]
function variants_estimate!(df::DataFrame)
    n_rows = nrow(df)
    @info "--------"
    @info "δ 株と非 δ 株のトレンド (trend_δ_lagXX, trend_nonδ_lagXX)"
    ts = df.date_value
    xs = df.trend
    for lag in likelihood_lags
        ps = δ_prop_maxL.(ts .+ lag)    # 各時刻での δ 株の割合（遅れ lag を見込む）
        xs_δ    = ps .* xs              # δ 株の感染確認者トレンド数
        xs_nonδ = xs .- xs_δ            # δ 株以外の感染確認者トレンド数
        lagname = @sprintf("lag%02d", lag)
        df[!,"trend_δ_"*lagname]    = xs_δ
        df[!,"trend_nonδ_"*lagname] = xs_nonδ
    end
    @info "--------"
    @info "非 δ 株拡大係数 (lgrowth_nonδ_lagXX)"
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_nonδ = df[!,"trend_nonδ_"*lagname]
        ξs_nonδ  = log.(xs_nonδ)
        ξds_nonδ = zeros(Float64, n_rows)
        for i in 1:n_rows
            if     i == 1
                ξds_nonδ[i] = ξs_nonδ[i+1] - ξs_nonδ[i]
            elseif i == n_rows
                ξds_nonδ[i] = ξs_nonδ[i] - ξs_nonδ[i-1]
            else
                ξds_nonδ[i] = (ξs_nonδ[i+1] - ξs_nonδ[i-1]) / 2.0
            end
        end
        df[!,"lgrowth_nonδ_"*lagname] = ξds_nonδ
    end
end

variants_estimate!(confirmed_df)

@info "--------"
@info "解析データ書き出し"
const output_csv_filename = "tokyo_confirmed_analysis4.csv"
@info output_csv_filename
CSV.write(output_csv_filename, confirmed_df; missingstring="NA")
@info size(confirmed_df)

@info "========"
@info "予測"

const predict_date_start = confirmed_df.date_mid[end]
const predict_date_end   = predict_date_start + Day(60)
@info predict_date_start, predict_date_end
function predict_confirmed(confirmed_df::DataFrame; β)
    @info "--------"
    @info "データフレーム作成"
    df = DataFrame()
    n_rows = Dates.value(predict_date_end - predict_date_start) + 1
    @info "n_rows", n_rows
    @info "--------"
    @info "日付 (date, date_value)"
    df.date = predict_date_start:Day(1):predict_date_end
    df.date_value = date_to_value.(df.date)
    @info size(df)
    @info "--------"
    @info "非 δ 株拡大係数設定 (lgrowth_nonδ_lagXX)"
    @info "過去 7 日間の平均. 値は時間に依らない定数"
    c_n_rows = nrow(confirmed_df)
    period = min(c_n_rows, 7)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        ξdname = "lgrowth_nonδ_"*lagname
        ξds_nonδ = confirmed_df[!,ξdname]
        ξd_nonδ_f = mean(ξds_nonδ[(c_n_rows-period+1):c_n_rows])
        @info ξdname, ξd_nonδ_f
        df[!,ξdname] = fill(ξd_nonδ_f, n_rows)
    end
    @info size(df)
    @info "--------"
    @info "拡大係数予測 (lgrowth_lagXX)"
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        ts = df.date_value
        ps = δ_prop_maxL.(ts .+ lag)
        ξds_nonδ = df[!,"lgrowth_nonδ_"*lagname]
        ξds = β .* ps .+ ξds_nonδ
        df[!,"lgrowth_"*lagname] = ξds
    end
    @info size(df)
    @info "--------"
    @info "非 δ 株感染確認数予測 (trend_nonδ_lagXX)"
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        ts = df.date_value
        ξds_nonδ = df[!,"lgrowth_nonδ_"*lagname]
        x_nonδ_1 = confirmed_df[end,"trend_nonδ_"*lagname]
        a_nonδ_1 = x_nonδ_1 / exp(ξds_nonδ[1] * ts[1])
        xs_nonδ = a_nonδ_1 .* exp.(ξds_nonδ .* ts)
        df[!,"trend_nonδ_"*lagname] = xs_nonδ
    end
    @info size(df)
    @info "--------"
    @info "感染確認数予測 (trend_lagXX)"
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        ts = df.date_value
        ps = δ_prop_maxL.(ts .+ lag)
        xs_nonδ = df[!,"trend_nonδ_"*lagname]
        xs = xs_nonδ ./ (1.0 .- ps)
        df[!,"trend_"*lagname] = xs
    end
    @info size(df)
    return df
end

predict_df = predict_confirmed(confirmed_df; β=likelihood_max_β)

@info "========"
@info "プロット"

pyplot(
    titlefont=font("Meiryo",10),
    guidefont=font("Meiryo",10),
    tickfont=font("Meiryo",9)
)
const plot_date_start = confirmed_df.date_mid[end] - Day(50)
const plot_date_end   = plot_date_start + Day(100)
#const plot_date_start = confirmed_df.date_mid[end] - Day(50) - Day(30)
#const plot_date_end   = plot_date_start + Day(100) + Day(30)
const date_ticks_dict = Dict(
    Date("2021-03-15") => "03\u201015",
    Date("2021-03-21") => "緊急事態宣言解除",
    Date("2021-04-01") => "2021\u201004\u201001",
    Date("2021-04-15") => "04\u201015",
    Date("2021-04-25") => "緊急事態宣言発令",
    Date("2021-05-01") => "2021\u201005\u201001",
    Date("2021-05-15") => "05\u201015",
    Date("2021-06-01") => "2021\u201006\u201001",
    Date("2021-06-15") => "06\u201015",
    Date("2021-06-20") => "緊急事態宣言解除",
    Date("2021-07-01") => "2021\u201007\u201001",
    Date("2021-07-04") => "都議会議員選挙",
    Date("2021-07-12") => "緊急事態宣言発令",
    Date("2021-07-15") => "07\u201015",
    Date("2021-07-23") => "オリンピック開会",
    Date("2021-08-01") => "2021\u201008\u201001",
    Date("2021-08-08") => "オリンピック閉会",
    Date("2021-08-15") => "08\u201015",
    Date("2021-08-22") => "宣言解除予定",
    Date("2021-09-01") => "2021\u201009\u201001",
)
const date_ticks_dates = sort(collect(keys(date_ticks_dict)))
const date_ticks_dict_en = Dict(
    Date("2021-03-21") => "End\\ of\\ the\\ SoE",
    Date("2021-04-01") => "2021\u201004\u201001",
    Date("2021-04-15") => "04\u201015",
    Date("2021-04-25") => "Decl.\\ of\\ the\\ SoE",
    Date("2021-05-01") => "2021\u201005\u201001",
    Date("2021-05-15") => "05\u201015",
    Date("2021-06-01") => "2021\u201006\u201001",
    Date("2021-06-15") => "06\u201015",
    Date("2021-06-20") => "End\\ of\\ the\\ SoE",
    Date("2021-07-01") => "2021\u201007\u201001",
    Date("2021-07-04") => "Assembly\\ election",
    Date("2021-07-12") => "Decl.\\ of\\ the\\ SoE",
    Date("2021-07-15") => "07\u201015",
    Date("2021-07-23") => "Olympics\\ starts",
    Date("2021-08-01") => "2021\u201008\u201001",
    Date("2021-08-08") => "Olympics\\ ends",
    Date("2021-08-15") => "08\u201015",
    Date("2021-08-22") => "SoE\\ ends\\ (sched)",
    Date("2021-09-01") => "2021\u201009\u201001",
)
const date_ticks_dates_en = sort(collect(keys(date_ticks_dict_en)))

plots = Dict{String,Plots.Plot}()
@info "--------"
@info "確認感染者数・トレンド"
for lag in likelihood_lags
    lagname = @sprintf("lag%02d", lag)
    title = "東京"
    if lag != 0
        title *= "（時間遅れ $(lag) 日）"
    end
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title=title,
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 1500.0)
    ylabel!(pl, "感染確認者数 [人/日]")
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    xs_nonδ_c = confirmed_df[!,"trend_nonδ_"*lagname]
    plot!(pl, ts_c, cs_c; label="感染確認者数", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="トレンド", lw=2, lc=:red3)
    plot!(pl, ts_c, xs_nonδ_c; label="非デルタ株（推定・予測）",
        frange=zeros(Float64, nrow(confirmed_df)), fa=0.3, c=:cyan)
    plot!(pl, ts_c, xs_c; label="デルタ株（推定・予測）",
        frange=xs_nonδ_c, fa=0.3, c=:red3)
    plot!(pl, ts_c, xs_nonδ_c; label=:none, lc=:cyan)
    plot!(pl, ts_c, xs_c; label=:none, lw=2, lc=:red3)
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    # 感染確認者予測
    ts_p = predict_df.date_value
    xs_p = predict_df[!,"trend_"*lagname]
    xs_nonδ_p = predict_df[!,"trend_nonδ_"*lagname]
    plot!(pl, ts_p, xs_nonδ_p; label=:none,
        frange=zeros(Float64, nrow(predict_df)), fa=0.1, c=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none,
        frange=xs_nonδ_p, fa=0.1, c=:red3, lw=2, ls=:dot)
    plot!(pl, ts_p, xs_nonδ_p; label=:none, lc=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none, lw=2, lc=:red3, ls=:dot)
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, 100, text("""
        予測は、直近 1 週間の拡大係数をもとに同等の対応が継続する（非デルタ株に対し一定）と想定
        デルタ株確認者数の予測は最尤線をたどるとした場合""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "confirmed_"*lagname
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "in English"
for lag in likelihood_lags
    lagname = @sprintf("lag%02d", lag)
    title = "Tokyo"
    if lag != 0
        title *= "(time lag: $(lag) days)"
    end
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title=title,
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates_en)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict_en[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 1500.0)
    ylabel!(pl, "Confirmed cases [cases/day]")
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    xs_nonδ_c = confirmed_df[!,"trend_nonδ_"*lagname]
    plot!(pl, ts_c, cs_c; label="Confirmed cases", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="Trend", lw=2, lc=:red3)
    plot!(pl, ts_c, xs_nonδ_c; label="non-δ cases (est./pred.)",
        frange=zeros(Float64, nrow(confirmed_df)), fa=0.3, c=:cyan)
    plot!(pl, ts_c, xs_c; label="δ cases (est./pred.)",
        frange=xs_nonδ_c, fa=0.3, c=:red3)
    plot!(pl, ts_c, xs_nonδ_c; label=:none, lc=:cyan)
    plot!(pl, ts_c, xs_c; label=:none, lw=2, lc=:red3)
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    # 感染確認者予測
    ts_p = predict_df.date_value
    xs_p = predict_df[!,"trend_"*lagname]
    xs_nonδ_p = predict_df[!,"trend_nonδ_"*lagname]
    plot!(pl, ts_p, xs_nonδ_p; label=:none,
        frange=zeros(Float64, nrow(predict_df)), fa=0.1, c=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none,
        frange=xs_nonδ_p, fa=0.1, c=:red3, lw=2, ls=:dot)
    plot!(pl, ts_p, xs_nonδ_p; label=:none, lc=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none, lw=2, lc=:red3, ls=:dot)
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, 100, text("""
        The predictions are based on the expansion coefficient for the recent week
        and assume that the equivalent response would continue.
        The prediction of the δ cases assumes that the maximum likelihood line is followed.""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "confirmed_"*lagname*"_en"
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "確認感染者数・トレンド（ラグ統合）"
for once in 1:1
    pl = plot(
        size=(640,700),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title="東京",
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 3500.0)
    ylabel!(pl, "感染確認者数 [人/日]")
    yticks!(pl,
        [0.0, 500.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 3500.0],
        ["0", "500", "1000", "1500", "2000", "2500", "3000", "3500"]
    )
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    plot!(pl, ts_c, cs_c; label="感染確認者数", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="トレンド", lw=2, lc=:red3)
    xs_nonδ_c = confirmed_df.trend_nonδ_lag00
    plot!(pl, ts_c, xs_nonδ_c; label="非デルタ株（推定・予測）",
        frange=zeros(Float64, nrow(confirmed_df)), fa=0.3, c=:cyan)
    plot!(pl, ts_c, xs_c; label="デルタ株（推定・予測）",
        frange=xs_nonδ_c, fa=0.3, c=:red3)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_nonδ_c_l = confirmed_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan)
        else
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan, la=0.7)
        end
    end
    plot!(pl, ts_c, xs_c; label=:none, lw=2, lc=:red3)
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    # 感染確認者予測
    ts_p = predict_df.date_value
    xs_p = predict_df.trend_lag00
    xs_nonδ_p = predict_df.trend_nonδ_lag00
    plot!(pl, ts_p, xs_nonδ_p; label=:none,
        frange=zeros(Float64, nrow(predict_df)), fa=0.1, c=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none,
        frange=xs_nonδ_p, fa=0.1, c=:red3, lw=2, ls=:dot)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_p_l = predict_df[!,"trend_"*lagname]
        xs_nonδ_p_l = predict_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot)
            plot!(pl, ts_p, xs_p_l; label=:none, lw=2, lc=:red3, ls=:dot)
        else
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot, la=0.7)
            plot!(pl, ts_p, xs_p_l; label=:none, lc=:red3, ls=:dot, la=0.7)
        end
    end
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, 100, text("""
        予測は、直近 1 週間の拡大係数をもとに同等の対応が継続する（非デルタ株に対し一定）と想定
        デルタ株確認者数の予測は最尤線をたどるとした場合""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "confirmed"
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "確認感染者数・トレンド（片対数グラフ）"
for once in 1:1
    pl = plot(
        size=(640,700),
#        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title="東京",
        legend=:topleft,
        fontfamily="Meiryo",
        yscale=:log
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 150.0, 8000.0)
#    ylims!(pl, 150.0, 4000.0)
    ylabel!(pl, "感染確認者数 [人/日]（対数目盛）")
    yticks!(pl,
        [150.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1500.0, 2000.0, 4000.0, 6000.0, 8000.0],
        ["150", "200", "400", "600", "800", "1000", "1500", "2000", "4000", "6000", "8000"]
#        [150.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1500.0, 2000.0, 4000.0],
#        ["150", "200", "400", "600", "800", "1000", "1500", "2000", "4000"]
    )
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    plot!(pl, ts_c, cs_c; label="感染確認者数", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="トレンド", lw=2, lc=:red3)
    xs_nonδ_c = confirmed_df.trend_nonδ_lag00
    plot!(pl, ts_c, xs_nonδ_c; label="非デルタ株（推定・予測）",
        frange=zeros(Float64, nrow(confirmed_df)), fa=0.3, c=:cyan)
    plot!(pl, ts_c, xs_c; label="デルタ株（推定・予測）",
        frange=xs_nonδ_c, fa=0.3, c=:red3)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_nonδ_c_l = confirmed_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan)
        else
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan, la=0.7)
        end
    end
    plot!(pl, ts_c, xs_c; label=:none, lw=2, lc=:red3)
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    # 感染確認者予測
    ts_p = predict_df.date_value
    xs_p = predict_df.trend_lag00
    xs_nonδ_p = predict_df.trend_nonδ_lag00
    plot!(pl, ts_p, xs_nonδ_p; label=:none,
        frange=zeros(Float64, nrow(predict_df)), fa=0.1, c=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none,
        frange=xs_nonδ_p, fa=0.1, c=:red3, lw=2, ls=:dot)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_p_l = predict_df[!,"trend_"*lagname]
        xs_nonδ_p_l = predict_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot)
            plot!(pl, ts_p, xs_p_l; label=:none, lw=2, lc=:red3, ls=:dot)
            annotate!(pl, date_to_value(plot_date_end) - 2.0, xs_nonδ_p_l[end-14],
                text("変異株報告に\n追加の遅延がないとき", font("Meiryo",7), RGB(0.3,0.6,0.6), :right))
        else
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot, la=0.7)
            plot!(pl, ts_p, xs_p_l; label=:none, lc=:red3, ls=:dot, la=0.7)
            annotate!(pl, date_to_value(plot_date_end) - 2.0, xs_nonδ_p_l[end-14],
                text("追加の遅延 $lag 日", font("Meiryo",7), RGB(0.3,0.6,0.6), :right))
        end
    end
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, 170, text("""
        予測は、直近 1 週間の拡大係数をもとに同等の対応が継続する（非デルタ株に対し一定）と想定
        デルタ株確認者数の予測は最尤線をたどるとした場合""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "confirmed_log"
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "in English"
for once in 1:1
    pl = plot(
        size=(640,750),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title="Tokyo",
        legend=:topleft,
        fontfamily="Meiryo",
        yscale=:log
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict_en[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 150.0, 6000.0)
    ylabel!(pl, "Confirmed cases [cases/day] (logscale)")
    yticks!(pl,
        [150.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1500.0, 2000.0, 4000.0, 6000.0],
        ["150", "200", "400", "600", "800", "1000", "1500", "2000", "4000", "6000"]
    )
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    plot!(pl, ts_c, cs_c; label="Confirmed cases", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="Trend", lw=2, lc=:red3)
    xs_nonδ_c = confirmed_df.trend_nonδ_lag00
    plot!(pl, ts_c, xs_nonδ_c; label="non-δ cases (est./pred.)",
        frange=zeros(Float64, nrow(confirmed_df)), fa=0.3, c=:cyan)
    plot!(pl, ts_c, xs_c; label="δ cases (est./pred.)",
        frange=xs_nonδ_c, fa=0.3, c=:red3)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_nonδ_c_l = confirmed_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan)
        else
            plot!(pl, ts_c, xs_nonδ_c_l; label=:none, lc=:cyan, la=0.7)
        end
    end
    plot!(pl, ts_c, xs_c; label=:none, lw=2, lc=:red3)
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    # 感染確認者予測
    ts_p = predict_df.date_value
    xs_p = predict_df.trend_lag00
    xs_nonδ_p = predict_df.trend_nonδ_lag00
    plot!(pl, ts_p, xs_nonδ_p; label=:none,
        frange=zeros(Float64, nrow(predict_df)), fa=0.1, c=:cyan, ls=:dot)
    plot!(pl, ts_p, xs_p; label=:none,
        frange=xs_nonδ_p, fa=0.1, c=:red3, lw=2, ls=:dot)
    for lag in likelihood_lags
        lagname = @sprintf("lag%02d", lag)
        xs_p_l = predict_df[!,"trend_"*lagname]
        xs_nonδ_p_l = predict_df[!,"trend_nonδ_"*lagname]
        if lag == 0
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot)
            plot!(pl, ts_p, xs_p_l; label=:none, lw=2, lc=:red3, ls=:dot)
            annotate!(pl, date_to_value(plot_date_end) - 2.0, xs_nonδ_p_l[end-14],
                text("no additional delay\nin reporting mutations", font("Meiryo",7), RGB(0.3,0.6,0.6), :right))
        else
            plot!(pl, ts_p, xs_nonδ_p_l; label=:none, lc=:cyan, ls=:dot, la=0.7)
            plot!(pl, ts_p, xs_p_l; label=:none, lc=:red3, ls=:dot, la=0.7)
            annotate!(pl, date_to_value(plot_date_end) - 2.0, xs_nonδ_p_l[end-14],
                text("$lag-day delay", font("Meiryo",7), RGB(0.3,0.6,0.6), :right))
        end
    end
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, 180, text("""
        The predictions are based on the expansion coefficient for the recent week
        and assume that the equivalent response would continue.
        The prediction of the δ cases assumes that the maximum likelihood line is followed.""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "confirmed_log_en"
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "拡大速度係数"
for lag in likelihood_lags
    lagname = @sprintf("lag%02d", lag)
    title = "東京"
    if lag != 0
        title *= "（変異株報告の追加の遅延 $(lag) 日）"
    end
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title=title,
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, -0.075, 0.150)
    ylabel!(pl, "拡大係数（時定数の逆数）[/日]")
    # グリッド
    plot!(pl, [sdate_value, edate_value], [0.0, 0.0]; label=:none, lc=RGB(0.3,0.3,0.3))
    # 拡大係数
    ts_c = confirmed_df.date_value
    ξds_c = confirmed_df.lgrowth
    ξds_nonδ_c = confirmed_df[!,"lgrowth_nonδ_"*lagname]
    plot!(pl, ts_c, ξds_c; label="拡大係数", lw=2, lc=:red3)
    plot!(pl, ts_c, ξds_nonδ_c; label="拡大係数（非デルタ株）推定", lc=:cyan)
    plot!(pl, ts_c, ξds_c; label=:none, lw=2, lc=:red3)
    # 拡大係数予測
    ts_p = predict_df.date_value
    ξds_p = predict_df[!,"lgrowth_"*lagname]
    ξds_nonδ_p = predict_df[!,"lgrowth_nonδ_"*lagname]
    plot!(pl, ts_p, ξds_nonδ_p; label="拡大係数（非デルタ株）想定値", lw=2, lc=:cyan, ls=:dot)
    plot!(pl, ts_p, ξds_p; label="デルタ株を加えた拡大係数予測", lw=2, lc=:red3, ls=:dot)
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, -0.06, text("""
        ここで拡大係数とは確認者数に指数関数をあてはめたときの時定数の逆数. 正で拡大、負で減少
        拡大係数が正のとき、その逆数の約 0.7 倍が拡大の倍加時間にあたる（0.05 なら倍加時間約 14 日）
        予測は、直近 1 週間の拡大係数をもとに同等の対応が継続する（非デルタ株に対し一定）と想定
        デルタ株拡大の予測は最尤線をたどるとした場合""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "growth_"*lagname
    plots[plname] = pl
    @info plname
end
@info "--------"
@info "拡大速度係数（週あたり拡大率）"
for lag in likelihood_lags
    lagname = @sprintf("lag%02d", lag)
    title = "東京"
    if lag != 0
        title *= "（変異株報告の追加の遅延 $(lag) 日）"
    end
    pl = plot(
        size=(640,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title=title,
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(plot_date_start)
    edate_value = date_to_value(plot_date_end)
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates[i]
        if plot_date_start ≤ date ≤ plot_date_end
            push!(vdates, date_to_value(date))
            push!(labels, date_ticks_dict[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, -0.075, 0.150)
    ylabel!(pl, "週あたり拡大率")
    yt_wrs = [0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8]
    yticks!(pl,
        [log(wr)/7.0 for wr in yt_wrs],
        ["$(wr)" for wr in yt_wrs]
    )
    # グリッド
    plot!(pl, [sdate_value, edate_value], [0.0, 0.0]; label=:none, lc=RGB(0.3,0.3,0.3))
    # 拡大係数
    ts_c = confirmed_df.date_value
    ξds_c = confirmed_df.lgrowth
    ξds_nonδ_c = confirmed_df[!,"lgrowth_nonδ_"*lagname]
    plot!(pl, ts_c, ξds_c; label="拡大係数", lw=2, lc=:red3)
    plot!(pl, ts_c, ξds_nonδ_c; label="拡大係数（非デルタ株）推定", lc=:cyan)
    plot!(pl, ts_c, ξds_c; label=:none, lw=2, lc=:red3)
    # 拡大係数予測
    ts_p = predict_df.date_value
    ξds_p = predict_df[!,"lgrowth_"*lagname]
    ξds_nonδ_p = predict_df[!,"lgrowth_nonδ_"*lagname]
    plot!(pl, ts_p, ξds_nonδ_p; label="拡大係数（非デルタ株）想定値", lw=2, lc=:cyan, ls=:dot)
    plot!(pl, ts_p, ξds_p; label="デルタ株を加えた拡大係数予測", lw=2, lc=:red3, ls=:dot)
    # アノテーション
    annotate!(pl, date_to_value(plot_date_start) + 2, -0.06, text("""
        ここで拡大係数とは確認者数に指数関数をあてはめたときの時定数の逆数. 正で拡大、負で減少
        拡大係数が正のとき、その逆数の約 0.7 倍が拡大の倍加時間にあたる（0.05 なら倍加時間約 14 日）
        予測は、直近 1 週間の拡大係数をもとに同等の対応が継続する（非デルタ株に対し一定）と想定
        デルタ株拡大の予測は最尤線をたどるとした場合""",
        font("Meiryo",7), RGB(0.3,0.3,0.3), :left))
    plname = "growth_m_"*lagname
    plots[plname] = pl
    @info plname
end

#=
@info "--------"
@info "ゴールデン・ウィーク"
const date_ticks_dict_gw = Dict(
    Date("2021-04-03") => "2021\u201004\u201003（土）",
    Date("2021-04-04") => "2021\u201004\u201004（日）",
    Date("2021-04-10") => "2021\u201004\u201010（土）",
    Date("2021-04-11") => "2021\u201004\u201011（日）",
    Date("2021-04-17") => "2021\u201004\u201017（土）",
    Date("2021-04-18") => "2021\u201004\u201018（日）",
    Date("2021-04-24") => "2021\u201004\u201024（土）",
    Date("2021-04-25") => "2021\u201004\u201025（日）",
    Date("2021-04-29") => "昭和の日（祝）",
    Date("2021-05-01") => "2021\u201005\u201001（土）",
    Date("2021-05-02") => "2021\u201005\u201002（日）",
    Date("2021-05-03") => "憲法記念日（祝）",
    Date("2021-05-04") => "みどりの日（祝）",
    Date("2021-05-05") => "こどもの日（祝）",
    Date("2021-05-08") => "2021\u201005\u201008（土）",
    Date("2021-05-09") => "2021\u201005\u201009（日）",
    Date("2021-05-15") => "2021\u201005\u201015（土）",
    Date("2021-05-16") => "2021\u201005\u201016（日）",
    Date("2021-05-22") => "2021\u201005\u201022（土）",
    Date("2021-05-23") => "2021\u201005\u201023（日）",
    Date("2021-05-29") => "2021\u201005\u201029（土）",
    Date("2021-05-30") => "2021\u201005\u201030（日）",
)
const date_ticks_dates_gw = sort(collect(keys(date_ticks_dict_gw)))
for once in 1:1
    pl = plot(
        size=(720,600),
        framestyle=:box,
        bottom_margin=8px, left_margin=16px, top_margin=16px, right_margin=32px,
        title="東京（ゴールデン・ウィーク連休時の休日効果）",
        legend=:topleft,
        fontfamily="Meiryo"
    )
    # x 軸設定
    sdate_value = date_to_value(Date("2021-04-14"))
    edate_value = date_to_value(Date("2021-05-27"))
    xlims!(pl, sdate_value, edate_value)
    nticks = length(date_ticks_dates_gw)
    vdates = Float64[]
    labels = String[]
    for i in 1:nticks
        date = date_ticks_dates_gw[i]
        if sdate_value ≤ date_to_value(date) ≤ edate_value
            push!(vdates, date_to_value(date) + 0.5)
            push!(labels, date_ticks_dict_gw[date])
        end
    end
    xticks!(pl, vdates, labels; xtickfontrotation=90.0)
    # y 軸設定
    ylims!(pl, 0.0, 1400.0)
    ylabel!(pl, "感染確認者数 [人/日]")
    yticks!(pl,
        [0.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1200.0, 1400.0],
        ["0", "200", "400", "600", "800", "1000", "1200", "1400"]
    )
    # 感染確認者数
    ts_c = confirmed_df.date_value
    cs_c = confirmed_df.count
    xs_c = confirmed_df.trend
    plot!(pl, ts_c, cs_c; label="感染確認者数", m=:circle, lc=:royalblue, mc=:royalblue1)
    plot!(pl, ts_c, xs_c; label="トレンド", lw=2, lc=:red3)
    xs_nonδ_c = confirmed_df.trend_nonδ_lag00
    plot!(pl, ts_c, cs_c; label=:none, m=:circle, lc=:royalblue, mc=:royalblue1)
    plname = "gw_holidays"
    plots[plname] = pl
    @info plname
end
=#

end #module

max_stars_repo_name: muimuiz/TokyoVariant

max_issues_repo_name: muimuiz/TokyoVariant

max_forks_repo_name: muimuiz/TokyoVariant


text:
using Qwt
using QuickStructs

function cpos(t)
  hyp = t
  θ = t * π / 20.
  x = cos(θ) * hyp
  y = sin(θ) * hyp
  P2(x,y)
end

function demo_anim()

  scene = currentScene()
  empty!(scene)
  background!(:gray)

  c = circle!(10)
  brush!(c, :black)

  T = 0.0 : 1.0 : 400.0
  for t in T
    position!(c, cpos(t))
    sleep(0.002)
  end

  scene
end

function demo_anim_lsm()

  scene = currentScene()
  empty!(scene)
  background!(:white)

  # create circles
  n = 100
  # colors = [makecolor(0,0,0,0) for i in 1:n]
  circles = SceneItem[pen!(brush!(circle!(10), :white), 0) for i in 1:n]

  T = 0.0 : 1.0 : 400.0
  for (i,t) in enumerate(T)

    # adjust position/alpha for circles
    for (j,c) in enumerate(circles)
      j > i && break
      tj = T[i-j+1]
      alpha = max(0., min(exp((tj-t) / 10), 1.0))
      # brush!(c, 0, 0, 0, alpha)
      # alpha!(c, alpha)
      # colors[j][:setAlphaF](alpha)
      # brush!(c, colors[j])
      brush!(c, 1-alpha, 1-alpha, 1-alpha)
      position!(c, cpos(tj))
    end

    sleep(0.000001)
  end

  scene
end

ext: jl


text:
@testsuite "interface" (AT, eltypes)->begin
    AT <: AbstractGPUArray || return

    N = 10
    x = AT(Vector{Int}(undef, N))
    x .= 0
    gpu_call(x) do ctx, x
        x[linear_index(ctx)] = 2
        return
    end
    @test all(x-> x == 2, Array(x))

    gpu_call(x; elements=N) do ctx, x
        x[linear_index(ctx)] = 2
        return
    end
    @test all(x-> x == 2, Array(x))
    gpu_call(x; threads=2, blocks=(N ÷ 2)) do ctx, x
        x[linear_index(ctx)] = threadidx(ctx)
        return
    end
    @test Array(x) == [1,2,1,2,1,2,1,2,1,2]

    gpu_call(x; threads=2, blocks=(N ÷ 2)) do ctx, x
        x[linear_index(ctx)] = blockidx(ctx)
        return
    end
    @test Array(x) == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5]
    x2 = AT([0])
    gpu_call(x2; threads=2, blocks=(N ÷ 2), target=x) do ctx, x
        x[1] = blockdim(ctx)
        return
    end
    @test Array(x2) == [2]

    gpu_call(x2; threads=2, blocks=(N ÷ 2), target=x) do ctx, x
        x[1] = griddim(ctx)
        return
    end
    @test Array(x2) == [5]

    gpu_call(x2; threads=2, blocks=(N ÷ 2), target=x) do ctx, x
        x[1] = global_size(ctx)
        return
    end
    @test Array(x2) == [10]
end


text:
# Solves the scalar, purely nonlinear, eigenvalue problem
# A(λ) = 0.2*sqrt(λ) - 0.6*sin(2*λ) with both a polynomial
# and fully rational approach

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

function nleigs_scalar()
    as_matrix(x::Number) = (M = Matrix{eltype(x)}(undef,1,1); M[1] = x; M)
    n = 1
    C = [as_matrix(0.2), as_matrix(-0.6)]
    f = [λ -> sqrt(λ), λ -> sin(2*λ)]
    nep = SPMF_NEP([C[1], C[2]], [f[1], f[2]])

    Σ = complex([0.01, 4])

    @bench @testset "Polynomial" begin
        lambda, X = nleigs(nep, Σ, logger=displaylevel, maxit=100, v=ones(n).+0im, leja=2, isfunm=false)

        # single eigenvalue converges
        verify_lambdas(1, nep, lambda, X)
    end

    @bench @testset "Fully rational" begin
        # set of poles candidates
        Ξ = -10 .^ range(-6, stop = 5, length = 10000)

        lambda, X = nleigs(nep, Σ, Ξ=Ξ, logger=displaylevel, maxit=100, v=ones(n).+0im, leja=2, isfunm=false)

        # three eigenvalues converge
        verify_lambdas(3, nep, lambda, X)
    end
end

@testset "NLEIGS: Scalar" begin
    nleigs_scalar()
end


text:
@nospecialize
Print = (
    indent = function Print_indent(p)
        function (io::IO, prefix::AbstractString)
            prefix = "  " * prefix
            p(io, prefix)
        end
    end,
    line = function Print_line(io::IO, prefix::AbstractString)
        println(io)
        print(io, prefix)
    end,
    w = function Print_word(s::AbstractString)
        function write(io::IO, ::AbstractString)
            print(io, s)
        end
    end,
    seq = function Print_seq(ps...)
        function (io::IO, prefix::AbstractString)
            prefix = prefix
            for p in ps
                p(io, prefix)
            end
        end
    end,
    run = function Print_run(io::IO, builder)
        builder(io, "")
    end,
)
@specialize


text:
#ConventionalApp: Develop conventional applications in Julia
#-------------------------------------------------------------------------------
module ConventionalApp

using Pkg
#using LibGit2

include("base.jl")
include("filegen.jl")
#include("install.jl")

export Project, setup_env, activeproject, @include_startup

end #ConventionalApp

max_stars_repo_name: ma-laforge/JuliaApplications

max_issues_repo_name: ma-laforge/JuliaApplications

max_forks_repo_name: ma-laforge/JuliaApplications


text:
module TestBind
using MXNet
using Base.Test

using ..Main: rand_dims, reldiff

################################################################################
# Test Implementations
################################################################################
function test_arithmetic{T <: mx.DType}(::Type{T}, uf, gf)
  shape = rand_dims()
  info("Bind::arithmetic::$T::$uf::dims = $shape")

  lhs = mx.Variable(:lhs)
  rhs = mx.Variable(:rhs)
  ret = uf(lhs, rhs)
  @test mx.list_arguments(ret) == [:lhs, :rhs]

  lhs_arr  = mx.NDArray(rand(T, shape))
  rhs_arr  = mx.NDArray(rand(T, shape))
  lhs_grad = mx.empty(T, shape)
  rhs_grad = mx.empty(T, shape)

  exec2 = mx.bind(ret, mx.Context(mx.CPU), [lhs_arr, rhs_arr], args_grad=[lhs_grad, rhs_grad])
  exec3 = mx.bind(ret, mx.Context(mx.CPU), [lhs_arr, rhs_arr])
  exec4 = mx.bind(ret, mx.Context(mx.CPU), Dict(:lhs=>lhs_arr, :rhs=>rhs_arr),
                  args_grad=Dict(:rhs=>rhs_grad, :lhs=>lhs_grad))

  mx.forward(exec2)
  mx.forward(exec3)
  mx.forward(exec4)

  out1 = uf(copy(lhs_arr), copy(rhs_arr))
  out2 = copy(exec2.outputs[1])
  out3 = copy(exec3.outputs[1])
  out4 = copy(exec4.outputs[1])
  @test isapprox(out1, out2)
  @test isapprox(out1, out3)
  @test isapprox(out1, out4)

  # test gradients
  out_grad = mx.NDArray(ones(T, shape))
  lhs_grad2, rhs_grad2 = gf(copy(out_grad), copy(lhs_arr), copy(rhs_arr))
  mx.backward(exec2, out_grad)
  @test isapprox(copy(lhs_grad), lhs_grad2)
  @test isapprox(copy(rhs_grad), rhs_grad2)

  # reset grads
  lhs_grad[:] = 0
  rhs_grad[:] = 0
  # compute using another binding
  mx.backward(exec4, out_grad)
  @test isapprox(copy(lhs_grad), lhs_grad2)
  @test isapprox(copy(rhs_grad), rhs_grad2)
end

function test_arithmetic()
  for T in [mx.fromTypeFlag(TF) for TF in instances(mx.TypeFlag)]
    test_arithmetic(T, .+, (g,x,y) -> (g,g))
    test_arithmetic(T, .-, (g,x,y) -> (g,-g))
    test_arithmetic(T, .*, (g,x,y) -> (y.*g, x.*g))
    if T <: Integer || T == Float16
      warn("Not running division test for $T")
    else
      test_arithmetic(T, ./, (g,x,y) -> (g ./ y, -x .* g ./ (y.^2)))
    end
  end
end

################################################################################
# Run tests
################################################################################
test_arithmetic()

end


text:
using SparseArrays
using ViewerGL
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#Octahedron DEBUG
octahedron = Lar.octahedron(1)
V,EV,FV,CV = Lar.struct2lar(octahedron)
GL.VIEW([ GL.GLGrid(V,FV), GL.GLFrame]);
#EV=collect(Set(EV))
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);

W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE]);
FVs = convert(Array{Lar.Cells}, triangulated_faces);
V = convert(Lar.Points, V');
GL.VIEW( GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1) );

EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
GL.VIEW( GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1) );

#Octahedron whit arbitrry height
octahedron = Lar.octahedron(1,1)
V,EV,FV,CV = Lar.struct2lar(octahedron)
GL.VIEW([ GL.GLGrid(V,FV, GL.Point4d(1,1,1,0.1)), GL.GLFrame]);

#Regular base octahedron
octahedron = Lar.rectangularbase_octahedron(1,2,1)
V,EV,FV,CV = Lar.struct2lar(octahedron)
GL.VIEW([ GL.GLGrid(V,FV, GL.Point4d(1,1,1,0.1)), GL.GLFrame]);


text:
# Perform actomyosin network simulations
# Alex Tam, 12/10/2020

# Load packages
using Revise # Modify code without restarting Julia
using Parameters # Tools for storing parameters in data structures
using LinearAlgebra # Matrix and vector operations
using StatsBase # Sampling without replacement
using Plots # Plotting library
using LaTeXStrings # Display LaTeX output
using ForwardDiff # Automatic differentiation tools
using Optim # Optimisation routines
using LineSearches # Line searches for optimisation
using Printf # C-style printing macros
using Distributions # Probability distributions
using DelimitedFiles # Write data to files

# Include code from external files
include("model_parameters.jl")
include("actomyosin_network.jl")
include("State.jl")
include("initial_condition.jl")
include("Actin_Filament.jl")
include("periodic.jl")
include("Cross_Link.jl")
include("Myosin_Motor.jl")
include("intersection_search.jl")
include("network_statistics.jl")
include("network_force.jl")
include("draw_network.jl")
include("optimise_network.jl")
include("energy.jl")

# Specify parameters
parN = Numerical_Parameters(); # Initialise struct of numerical parameters
parA = Actin_Properties(nSeg = 1, LSeg = 1.0); # Initialise struct of actin filament properties
parM = Myosin_Properties(); # Initialise struct of myosin motor properties

# Run simulations
@time state, af, mm, Force, Curvature, Index = actomyosin_network(parN, parA, parM);


text:
using Test

# QuadTS is for numerical integration in [-1, 1]
@testset "QuadTS" begin
    include("quadts.jl")
end

# QuadES is for numerical integration in [0, ∞]
@testset "QuadES" begin
    include("quades.jl")
end

# QuadSS is for numerical integration in [-∞, ∞]
@testset "QuadSS" begin
    include("quadss.jl")
end

# QuadDE is the general purpose numerical integrator
# NOTE: If one wants to know how to use, the only files they should check are
#       the below items.
@testset "QuadDE" begin
    @testset "quadde DEMO" begin
        include("quadde/Mori.jl")
    end

    @testset "Varius intervals" begin
        include("quadde/intervals.jl")
    end

    @testset "Basic principles" begin
        include("quadde/principles.jl")
    end

    @testset "Kahaner's problems" begin
        include("quadde/Kahaner.jl")
    end
end

# quaddeo is for the numerical integration of oscillatory functions
@testset "quaddeo" begin
    include("quaddeo.jl")
end


text:
# helper.jl
# A lot of helper functions

using LinearAlgebra
import OffsetArrays
using ImageFiltering: BorderArray, Fill, Pad
using FFTW

Power2 = x -> 2^(ceil(Int, log2(x)))
_padup(nx, nx_psf) = ceil(Int, (Power2(nx + nx_psf - 1) - nx) / 2)
_paddown(nx, nx_psf) = floor(Int, (Power2(nx + nx_psf - 1) - nx) / 2)
_padleft(nz, nz_psf) = ceil(Int, (Power2(nz + nz_psf - 1) - nz) / 2)
_padright(nz, nz_psf) = floor(Int, (Power2(nz + nz_psf - 1) - nz) / 2)

"""
    padzero!(output, img, pad_x, pad_y)
Mutating version of padding a 2D image by filling zeros.
Output has size `(size(img, 1) + padsize[1] + padsize[2], size(img, 2) + padsize[3] + padsize[4])`.
"""
Base.@propagate_inbounds function padzero!(
    output::AbstractMatrix{T},
    img::AbstractMatrix,
    padsize::NTuple{4, <:Int}, # up, down, left, right
) where {T}

    @boundscheck size(output) ==
        size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size")
    M, N = size(img)
    output .= zero(T)
    for j = padsize[3] + 1:padsize[3] + N, i = padsize[1] + 1:padsize[1] + M
        @inbounds output[i, j] = img[i - padsize[1], j - padsize[3]]
    end
    return output
end

"""
    padrepl!(output, img, padsize)
Pad with replication from `img` into `output`
"""
Base.@propagate_inbounds function padrepl!(
    output::AbstractMatrix,
    img::AbstractMatrix,
    padsize::NTuple{4, <:Int}, # up, down, left, right
)

    @boundscheck size(output) ==
        size(img) .+ (padsize[1] + padsize[2], padsize[3] + padsize[4]) || throw("size")
    M, N = size(img)
    for j = 1:size(output, 2), i = 1:size(output, 1)
        @inbounds output[i, j] = img[clamp(i - padsize[1], 1, M), clamp(j - padsize[3], 1, N)]
    end
    return output
end

"""
    pad2sizezero!(output, img, padsize)

Non-allocating version of padding:
`output[pad_dims[1]+1 : pad_dims[1]+dims[1],
        pad_dims[2]+1 : pad_dims[2]+dims[2]]) .= img
"""
Base.@propagate_inbounds function pad2sizezero!(
    output::AbstractMatrix{T},
    img::AbstractMatrix,
    padsize::Tuple{<:Int, <:Int},
) where {T}
    @boundscheck size(output) == padsize || throw("size")

    dims = size(img)
    pad_dims = ceil.(Int, (padsize .- dims) ./ 2)
    output .= zero(T)
    for j = pad_dims[2]+1:pad_dims[2]+dims[2], i = pad_dims[1]+1:pad_dims[1]+dims[1]
        @inbounds output[i, j] = img[i - pad_dims[1], j - pad_dims[2]]
    end
    return output
end

"""
    pad_it!(X, padsize)
Zero-pad `X` to `padsize`
"""
function pad_it!(X::AbstractArray{T,D}, padsize::NTuple{D,<:Int}) where {D, T <: Number}
    dims = size(X)
    return OffsetArrays.no_offset_view(
        BorderArray(X,
            Fill(zero(T),
               (ceil.(Int, (padsize .- dims) ./ 2)),
               (floor.(Int, (padsize .- dims) ./ 2)),
            )
        )
    )
end

fftshift!(dst::AbstractArray, src::AbstractArray) = circshift!(dst, src, size(src) .÷ 2)

ifftshift!(dst::AbstractArray, src::AbstractArray) = circshift!(dst, src, size(src) .÷ -2)

"""
    fftshift2!(dst, src)
Same as `fftshift` in 2d, but non-allocating
"""
Base.@propagate_inbounds function fftshift2!(
    dst::AbstractMatrix,
    src::AbstractMatrix,
)
    @boundscheck (iseven(size(src, 1)) && iseven(size(src, 2))) || throw("odd")
    @boundscheck size(src) == size(dst) || throw("size")
    m, n = div.(size(src), 2)
    for j = 1:n, i = 1:m
        @inbounds dst[i, j] = src[m+i, n+j]
    end
    for j = n+1:2n, i = 1:m
        @inbounds dst[i, j] = src[m+i, j-n]
    end
    for j = 1:n, i = m+1:2m
        @inbounds dst[i, j] = src[i-m, j+n]
    end
    for j = n+1:2n, i = m+1:2m
        @inbounds dst[i, j] = src[i-m, j-n]
    end
    return dst
end

"""
    plus1di!(mat2d, mat1d)
Non-allocating `mat2d[i, :] += mat1d`
"""
Base.@propagate_inbounds function plus1di!(
    mat2d::AbstractMatrix,
    mat1d::AbstractVector,
    i::Int,
)
    @boundscheck (size(mat1d, 1) == size(mat2d, 2) || throw("size2"))
    @boundscheck (1 ≤ i ≤ size(mat2d, 1) || throw("bad i"))
    for m in 1:size(mat2d, 2)
        @inbounds mat2d[i, m] += mat1d[m]
    end
    return mat2d
end

"""
    plus1dj!(mat2d, mat1d)
Non-allocating `mat2d[:, j] += mat1d`
"""
Base.@propagate_inbounds function plus1dj!(
    mat2d::AbstractMatrix,
    mat1d::AbstractVector,
    j::Int,
)
    @boundscheck (size(mat1d, 1) == size(mat2d, 1) || throw("size1"))
    @boundscheck (1 ≤ j ≤ size(mat2d, 2) || throw("bad j"))
    for n in 1:size(mat2d, 1)
        @inbounds mat2d[n, j] += mat1d[n]
    end
    return mat2d
end

"""
    plus2di!(mat1d, mat2d, i)
Non-allocating `mat1d += mat2d[i,:]`
"""
Base.@propagate_inbounds function plus2di!(
    mat1d::AbstractVector,
    mat2d::AbstractMatrix,
    i::Int,
)
    @boundscheck (size(mat1d, 1) == size(mat2d, 2) || throw("size2"))
    @boundscheck (1 ≤ i ≤ size(mat2d, 1) || throw("bad i"))
    for m in 1:size(mat2d, 2)
        @inbounds mat1d[m] += mat2d[i, m]
    end
    return mat1d
end

"""
    plus2dj!(mat1d, mat2d, j)
Non-allocating `mat1d += mat2d[:,j]`
"""
Base.@propagate_inbounds function plus2dj!(
    mat1d::AbstractVector,
    mat2d::AbstractMatrix,
    j::Int,
)
    @boundscheck (size(mat1d, 1) == size(mat2d, 1) || throw("size1"))
    @boundscheck (1 ≤ j ≤ size(mat2d, 2) || throw("bad j"))
    for n in 1:size(mat2d, 1)
        @inbounds mat1d[n] += mat2d[n, j]
    end
    return mat1d
end

"""
    plus3di!(mat2d, mat3d, i)
Non-allocating `mat2d += mat3d[i,:,:]`
"""
Base.@propagate_inbounds function plus3di!(
    mat2d::AbstractMatrix,
    mat3d::AbstractArray{<:Any, 3},
    i::Int,
)
    @boundscheck (size(mat2d, 1) == size(mat3d, 2) || throw("size2"))
    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
    @boundscheck (1 ≤ i ≤ size(mat3d, 1) || throw("bad i"))
    for n in 1:size(mat2d, 2), m in 1:size(mat2d, 1)
        @inbounds mat2d[m, n] += mat3d[i, m, n]
    end
    return mat2d
end

"""
    plus3dj!(mat2d, mat3d, j)
Non-allocating `mat2d += mat3d[:,j,:]`
"""
Base.@propagate_inbounds function plus3dj!(
    mat2d::AbstractMatrix,
    mat3d::AbstractArray{<:Any, 3},
    j::Int,
)
    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
    for n in 1:size(mat2d, 2), m in 1:size(mat2d, 1)
        @inbounds mat2d[m, n] += mat3d[m, j, n]
    end
    return mat2d
end

"""
    plus3dk!(mat2d, mat3d, k)
Non-allocating `mat2d += mat3d[:,:,k]`
"""
Base.@propagate_inbounds function plus3dk!(
    mat2d::AbstractMatrix,
    mat3d::AbstractArray{<:Any, 3},
    k::Int,
)
    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
    @boundscheck (size(mat2d, 2) == size(mat3d, 2) || throw("size2"))
    @boundscheck (1 ≤ k ≤ size(mat3d, 3) || throw("bad k"))
    for n in 1:size(mat2d, 2), m in 1:size(mat2d, 1)
        @inbounds mat2d[m, n] += mat3d[m, n, k]
    end
    return mat2d
end

"""
    scale3dj!(mat2d, mat3d, j, s)
Non-allocating `mat2d = s * mat3d[:,j,:]`
"""
Base.@propagate_inbounds function scale3dj!(
    mat2d::AbstractMatrix,
    mat3d::AbstractArray{<:Any, 3},
    j::Int,
    s::RealU,
)
    @boundscheck (size(mat2d, 1) == size(mat3d, 1) || throw("size1"))
    @boundscheck (size(mat2d, 2) == size(mat3d, 3) || throw("size3"))
    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
    for n in 1:size(mat2d, 2), m in 1:size(mat2d, 1)
        @inbounds mat2d[m, n] = s * mat3d[m, j, n]
    end
    return mat2d
end

"""
    mul3dj!(mat3d, mat2d, j)
Non-allocating `mat3d[:,j,:] *= mat2d`
"""
Base.@propagate_inbounds function mul3dj!(
    mat3d::AbstractArray{<:Any, 3},
    mat2d::AbstractMatrix,
    j::Int,
)
    @boundscheck (size(mat3d, 1) == size(mat2d, 1) || throw("size1"))
    @boundscheck (size(mat3d, 3) == size(mat2d, 2) || throw("size3"))
    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
    for n in 1:size(mat3d, 3), m in 1:size(mat3d, 1)
        @inbounds mat3d[m, j, n] *= mat2d[m, n]
    end
    return mat3d
end

"""
    copy3dj!(mat2d, mat3d, j)
Non-allocating `mat2d .= mat3d[:,j,:]`
"""
Base.@propagate_inbounds function copy3dj!(
    mat2d::AbstractMatrix,
    mat3d::AbstractArray{<:Any, 3},
    j::Int,
)
    @boundscheck (size(mat3d, 1) == size(mat2d, 1) || throw("size1"))
    @boundscheck (size(mat3d, 3) == size(mat2d, 2) || throw("size3"))
    @boundscheck (1 ≤ j ≤ size(mat3d, 2) || throw("bad j"))
    for n in 1:size(mat3d, 3), m in 1:size(mat3d, 1)
        @inbounds mat2d[m,n] = mat3d[m, j, n]
    end
    return mat2d
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "quadruped", "visuals.jl"))
vis = Visualizer()
open(vis)
# const ContactImplicitMPC = Main

s = get_simulation("quadruped", "flat_2D_lc", "flat")
model = s.model
env = s.env

ref_traj_ = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt))
ref_traj = deepcopy(ref_traj_)
update_friction_coefficient!(ref_traj, model, env)

# time
H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 1000 #4000 #3000

# barrier parameter
κ_mpc = 2.0e-4

obj = TrackingObjective(model, env, H_mpc,
    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.dim.q-3)]) for t = 1:H_mpc],
    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
    n_opts = NewtonOptions(
		solver = :lu_solver,
		r_tol = 3e-4,
		max_iter = 5,
		max_time = ref_traj.h, # HARD REAL TIME
		),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        # altitude_update = true,
        # altitude_impact_threshold = 0.05,
        # altitude_verbose = true,
        ),
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					γ_reg = 0.1,
					κ_tol = κ_mpc,
					r_tol = 1.0e-8,
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1000.0,
					),
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

sim = simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	ip_opts = InteriorPointOptions(
		undercut = Inf,
		γ_reg = 0.0,
		# verbose = true,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = SimulatorOptions(warmstart = true),
    )

@time status = ContactImplicitMPC.simulate!(sim, verbose = true)

################################################################################
# Timing result
################################################################################
process!(sim)
# Time budget
ref_traj.h
# Time used on average
sim.stats.μ_dt
# Speed ratio
H_sim * h_sim / sum(sim.stats.dt)

# @elapsed status = ContactImplicitMPC.simulate!(sim)
# @profiler status = ContactImplicitMPC.simulate!(sim)

plot_lines!(vis, model, sim.traj.q[1:1:end])
plot_surface!(vis, env, ylims=[0.3, -0.05])
anim = visualize_meshrobot!(vis, model, sim.traj, sample=5)
# anim = visualize_robot!(vis, model, sim.traj, anim=anim)
anim = visualize_force!(vis, model, env, sim.traj, anim=anim, h=h_sim, sample = 5)

# Display ghosts
t_ghosts = [1, 1333, 2666]
mvis_ghosts = []
for (i,t) in enumerate(t_ghosts)
    α = i/(length(t_ghosts)+1)
    name = Symbol("ghost$i")
    mvis = build_meshrobot!(vis, model, name=name, α=α)
    push!(mvis_ghosts, mvis)
end

for (i,t) in enumerate(t_ghosts)
    name = Symbol("ghost$i")
    set_meshrobot!(vis, mvis_ghosts[i], model, sim.traj.q[t], name=name)
end

# R = Dict{Symbol,Vector}(
# 	:dyn => [],
# 	:imp => [],
# 	:mdp => [],
# 	:fri => [],
# 	:bimp => [],
# 	:bmdp => [],
# 	:bfri => [],
# 	)
# nz = num_var(model, env)
# dyn = index_dyn(model, env, quat = true)
# imp = index_imp(model, env, quat = true)
# mdp = index_mdp(model, env, quat = true)
# fri = index_fri(model, env, quat = true)
# bimp = index_bimp(model, env, quat = true)
# bmdp = index_bmdp(model, env, quat = true)
# bfri = index_bfri(model, env, quat = true)
# for t = 1:ref_traj.H
#
# 	z = deepcopy(ref_traj.z[t])
# 	θ = deepcopy(ref_traj.θ[t])
# 	r = zeros(nz)
# 	r = residual(model, env, z, θ, 0.0)
# 	s.res.r!(r, z, θ, 1e-6)
# 	push!(R[:dyn], norm(r[dyn]))
# 	push!(R[:imp], norm(r[imp]))
# 	push!(R[:mdp], norm(r[mdp]))
# 	push!(R[:fri], norm(r[fri]))
# 	push!(R[:bimp], norm(r[bimp]))
# 	push!(R[:bmdp], norm(r[bmdp]))
# 	push!(R[:bfri], norm(r[bfri]))
# end
# plt = plot()
# plot!(plt, R[:dyn], label = "dyn")
# plot!(plt, R[:imp], label = "imp")
# plot!(plt, R[:mdp], label = "mdp")
# plot!(plt, R[:fri], label = "fri")
# plot!(plt, R[:bimp], label = "bimp")
# plot!(plt, R[:bmdp], label = "bmdp")
# plot!(plt, R[:bfri], label = "bfri")


text:
abstract type AbstractThunk <: AbstractTangent end

struct MutateThunkException <: Exception end

function Base.showerror(io::IO, e::MutateThunkException)
    print(io, "Tried to mutate a thunk, this is not supported. `unthunk` it first.")
    return nothing
end

Base.Broadcast.broadcastable(x::AbstractThunk) = broadcastable(unthunk(x))

@inline function Base.iterate(x::AbstractThunk)
    val = unthunk(x)
    element, state = iterate(val)
    return element, (val, state)
end

@inline function Base.iterate(::AbstractThunk, (underlying_object, state))
    next = iterate(underlying_object, state)
    next === nothing && return nothing
    element, new_state = next
    return element, (underlying_object, new_state)
end

Base.:(==)(a::AbstractThunk, b::AbstractThunk) = unthunk(a) == unthunk(b)
Base.:(==)(a::AbstractThunk, b) = unthunk(a) == b
Base.:(==)(a, b::AbstractThunk) = a == unthunk(b)

Base.:(-)(a::AbstractThunk) = -unthunk(a)
Base.:(-)(a::AbstractThunk, b) = unthunk(a) - b
Base.:(-)(a, b::AbstractThunk) = a - unthunk(b)
Base.:(/)(a::AbstractThunk, b) = unthunk(a) / b
Base.:(/)(a, b::AbstractThunk) = a / unthunk(b)

Base.real(a::AbstractThunk) = real(unthunk(a))
Base.imag(a::AbstractThunk) = imag(unthunk(a))
Base.Complex(a::AbstractThunk) = Complex(unthunk(a))
Base.Complex(a::AbstractThunk, b::AbstractThunk) = Complex(unthunk(a), unthunk(b))

Base.mapreduce(f, op, a::AbstractThunk; kws...) = mapreduce(f, op, unthunk(a); kws...)
function Base.mapreduce(f, op, itr, a::AbstractThunk; kws...)
    return mapreduce(f, op, itr, unthunk(a); kws...)
end
Base.sum(a::AbstractThunk; kws...) = sum(unthunk(a); kws...)
Base.sum!(r, A::AbstractThunk; kws...) = sum!(r, unthunk(A); kws...)

Base.fill(a::AbstractThunk, b::Integer) = fill(unthunk(a), b)
Base.vec(a::AbstractThunk) = vec(unthunk(a))
Base.reshape(a::AbstractThunk, args...) = reshape(unthunk(a), args...)
Base.getindex(a::AbstractThunk, args...) = getindex(unthunk(a), args...)
Base.setindex!(a::AbstractThunk, value, key...) = throw(MutateThunkException())
Base.selectdim(a::AbstractThunk, args...) = selectdim(unthunk(a), args...)

LinearAlgebra.Array(a::AbstractThunk) = Array(unthunk(a))
LinearAlgebra.Matrix(a::AbstractThunk) = Matrix(unthunk(a))
LinearAlgebra.Diagonal(a::AbstractThunk) = Diagonal(unthunk(a))
LinearAlgebra.LowerTriangular(a::AbstractThunk) = LowerTriangular(unthunk(a))
LinearAlgebra.UpperTriangular(a::AbstractThunk) = UpperTriangular(unthunk(a))
LinearAlgebra.Symmetric(a::AbstractThunk, uplo=:U) = Symmetric(unthunk(a), uplo)
LinearAlgebra.Hermitian(a::AbstractThunk, uplo=:U) = Hermitian(unthunk(a), uplo)

function LinearAlgebra.diagm(
    kv::Pair{<:Integer,<:AbstractThunk}, kvs::Pair{<:Integer,<:AbstractThunk}...
)
    return diagm((k => unthunk(v) for (k, v) in (kv, kvs...))...)
end
function LinearAlgebra.diagm(
    m, n, kv::Pair{<:Integer,<:AbstractThunk}, kvs::Pair{<:Integer,<:AbstractThunk}...
)
    return diagm(m, n, (k => unthunk(v) for (k, v) in (kv, kvs...))...)
end

LinearAlgebra.tril(a::AbstractThunk) = tril(unthunk(a))
LinearAlgebra.tril(a::AbstractThunk, k) = tril(unthunk(a), k)
LinearAlgebra.triu(a::AbstractThunk) = triu(unthunk(a))
LinearAlgebra.triu(a::AbstractThunk, k) = triu(unthunk(a), k)
LinearAlgebra.tr(a::AbstractThunk) = tr(unthunk(a))
LinearAlgebra.cross(a::AbstractThunk, b) = cross(unthunk(a), b)
LinearAlgebra.cross(a, b::AbstractThunk) = cross(a, unthunk(b))
LinearAlgebra.cross(a::AbstractThunk, b::AbstractThunk) = cross(unthunk(a), unthunk(b))
LinearAlgebra.dot(a::AbstractThunk, b) = dot(unthunk(a), b)
LinearAlgebra.dot(a, b::AbstractThunk) = dot(a, unthunk(b))
LinearAlgebra.dot(a::AbstractThunk, b::AbstractThunk) = dot(unthunk(a), unthunk(b))

LinearAlgebra.ldiv!(a, b::AbstractThunk) = throw(MutateThunkException())
LinearAlgebra.rdiv!(a::AbstractThunk, b) = throw(MutateThunkException())

LinearAlgebra.mul!(A, B::AbstractThunk, C) = mul!(A, unthunk(B), C)
LinearAlgebra.mul!(C::AbstractThunk, A, B, α, β) = throw(MutateThunkException())
function LinearAlgebra.mul!(C::AbstractThunk, A::AbstractThunk, B, α, β)
    return throw(MutateThunkException())
end
function LinearAlgebra.mul!(C::AbstractThunk, A, B::AbstractThunk, α, β)
    return throw(MutateThunkException())
end
function LinearAlgebra.mul!(C::AbstractThunk, A::AbstractThunk, B::AbstractThunk, α, β)
    return throw(MutateThunkException())
end
LinearAlgebra.mul!(C, A::AbstractThunk, B, α, β) = mul!(C, unthunk(A), B, α, β)
LinearAlgebra.mul!(C, A, B::AbstractThunk, α, β) = mul!(C, A, unthunk(B), α, β)
function LinearAlgebra.mul!(C, A::AbstractThunk, B::AbstractThunk, α, β)
    return mul!(C, unthunk(A), unthunk(B), α, β)
end

function LinearAlgebra.BLAS.ger!(alpha, x::AbstractThunk, y, A)
    return LinearAlgebra.BLAS.ger!(alpha, unthunk(x), y, A)
end
function LinearAlgebra.BLAS.ger!(alpha, x, y::AbstractThunk, A)
    return LinearAlgebra.BLAS.ger!(alpha, x, unthunk(y), A)
end
function LinearAlgebra.BLAS.gemv!(tA, alpha, A, x::AbstractThunk, beta, y)
    return LinearAlgebra.BLAS.gemv!(tA, alpha, A, unthunk(x), beta, y)
end
function LinearAlgebra.BLAS.gemv(tA, alpha, A, x::AbstractThunk)
    return LinearAlgebra.BLAS.gemv(tA, alpha, A, unthunk(x))
end
function LinearAlgebra.BLAS.scal!(n, a::AbstractThunk, X, incx)
    return LinearAlgebra.BLAS.scal!(n, unthunk(a), X, incx)
end

function LinearAlgebra.LAPACK.trsyl!(transa, transb, A, B, C::AbstractThunk, isgn=1)
    return throw(MutateThunkException())
end

"""
    @thunk expr

Define a [`Thunk`](@ref) wrapping the `expr`, to lazily defer its evaluation.
"""
macro thunk(body)
    # Basically `:(Thunk(() -> $(esc(body))))` but use the location where it is defined.
    # so we get useful stack traces if it errors.
    func = Expr(:->, Expr(:tuple), Expr(:block, __source__, body))
    return :(Thunk($(esc(func))))
end

"""
    unthunk(x)

On `AbstractThunk`s this removes 1 layer of thunking.
On any other type, it is the identity operation.
"""
@inline unthunk(x) = x

Base.conj(x::AbstractThunk) = @thunk(conj(unthunk(x)))
Base.adjoint(x::AbstractThunk) = @thunk(adjoint(unthunk(x)))
Base.transpose(x::AbstractThunk) = @thunk(transpose(unthunk(x)))

#####
##### `Thunk`
#####

"""
    Thunk(()->v)
A thunk is a deferred computation.
It wraps a zero argument closure that when invoked returns a differential.
`@thunk(v)` is a macro that expands into `Thunk(()->v)`.

To evaluate the wrapped closure, call [`unthunk`](@ref) which is a no-op when the
argument is not a `Thunk`.

```jldoctest
julia> t = @thunk(3)
Thunk(var"#4#5"())

julia> unthunk(t)
3
```

### When to `@thunk`?
When writing `rrule`s (and to a lesser exent `frule`s), it is important to `@thunk`
appropriately.
Propagation rules that return multiple derivatives may not have all deriviatives used.
 By `@thunk`ing the work required for each derivative, they then compute only what is needed.

#### How do thunks prevent work?
If we have `res = pullback(...) = @thunk(f(x)), @thunk(g(x))`
then if we did `dx + res[1]` then only `f(x)` would be evaluated, not `g(x)`.
Also if we did `ZeroTangent() * res[1]` then the result would be `ZeroTangent()` and `f(x)` would not be evaluated.

#### So why not thunk everything?
`@thunk` creates a closure over the expression, which (effectively) creates a `struct`
with a field for each variable used in the expression, and call overloaded.

Do not use `@thunk` if this would be equal or more work than actually evaluating the expression itself.
This is commonly the case for scalar operators.

For more details see the manual section [on using thunks effectively](http://www.juliadiff.org/ChainRulesCore.jl/dev/writing_good_rules.html#Use-Thunks-appropriately-1)
"""
struct Thunk{F} <: AbstractThunk
    f::F
end

@inline unthunk(x::Thunk) = x.f()

Base.show(io::IO, x::Thunk) = print(io, "Thunk($(repr(x.f)))")

Base.convert(::Type{<:Thunk}, a::AbstractZero) = @thunk(a)

"""
    InplaceableThunk(add!::Function, val::Thunk)

A wrapper for a `Thunk`, that allows it to define an inplace `add!` function.

`add!` should be defined such that: `ithunk.add!(Δ) = Δ .+= ithunk.val`
but it should do this more efficently than simply doing this directly.
(Otherwise one can just use a normal `Thunk`).

Most operations on an `InplaceableThunk` treat it just like a normal `Thunk`;
and destroy its inplacability.
"""
struct InplaceableThunk{T<:Thunk,F} <: AbstractThunk
    add!::F
    val::T
end

unthunk(x::InplaceableThunk) = unthunk(x.val)

function Base.show(io::IO, x::InplaceableThunk)
    return print(io, "InplaceableThunk($(repr(x.val)), $(repr(x.add!)))")
end


text:
# ----------------------------------------------------
# --- Digital Communication module
# ----------------------------------------------------
__precompile__()

module DigitalComm
# ----------------------------------------------------
# --- Modules
# ----------------------------------------------------
# Only submodules are considered so don't need to overload too much packages here
# --> Go to next section with submodules loading
using FFTW
using Statistics
using DSP

# ----------------------------------------------------
# --- Submodules inclusion
# ----------------------------------------------------
# --- Binary managment
include("genBitSequence.jl");
# Export
export genBitSequence!  , genBitSequence;
export genByteSequence! , genByteSequence;

# --- QAM bit mapping
include("bitMapping.jl");
export bitMappingQAM! , bitMappingQAM;

# --- QAM Hard demapper
include("bitDeMapping.jl");
export bitDemappingQAM! ,  bitDemappingQAM;
include("hardConstellation.jl");
export hardConstellation! , hardConstellation;

# --- QAM soft demapper
include("symbolDemapper.jl");
export symbolDemappingQAM! , symbolDemappingQAM;
export llrToHardBits! , llrToHardBits;
export llrToUInt! , llrToUInt;

# --- AWGN channel model
include("Channel/addNoise.jl");
export addNoise, addNoise!
include("Channel/rayleighChan.jl");
export rayleighChan;
include("Channel/getChannel.jl");
export initChannel
export getChannel
export getFIRResponse
export applyChannel
include("Channel/addCFO.jl")
export addCFO
export addCFO!

# --- Raised cosine filters
include("raisedCosine.jl")
export raisedCosine
export sqrtRaisedCosine

# ----------------------------------------------------
# --- Function definition
# ----------------------------------------------------
# --- Q function definition
import SpecialFunctions.erfc
"""
---
Returns the Q function used for Bit error rate computation in digital system
 Q(x)= 1/2 erfc(x/sqrt(2))
erfc is the Complexmplementary error function, i.e. the accurate version of 1-erf(x) for large x
erfc is inherited from DSP
# --- Syntax
y = qfunc(x)
# --- Input parameters
- x: Input [Float64]
# --- Output parameters
- y: Q(x)[Float64]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function qFunc(x)
	return 1/2 * erfc.( x / sqrt(2));
end
export qFunc;

"""
---
Returns the Signal to interference ratio expressed in dB (or in linear) between a obersvation signal d(n) and a reference signal u(n)
The ratio is expressed as 10*log10( E[ || d(n) - u(n) || / E[||u(n)||^2]  )
 with E the expectation wrt to time
 The 2 vectors d and u should have the same length L
# --- Syntax
  sir = getSIR( d, u , type="dB")
# ---  Input parameter
- d	: Observation signal [Array{Any}]
- u	: Reference signal [Array{Any}]
- type: Output unit [String]: "dB" or "Linear" (default, "dB")
# --- Output parameters
- sir	: Signal to interference ratio in unit `type`
# ---
# v 1.0 - Robin Gerzaguet.
"""
function getSIR(d,u,type="dB")
	# --- Setting type
	if lowercase(type) == "linear"
		# --- flag for casting
		castype = 1;
	else
		# --- By default SIR is in dB
		castype = 0;
	end
	# --- Moving vector to column system
	ur	= u[:];
	dr  = d[:];
	# --- Checking length
	if length(ur) != length(dr)
		# --- Raise an error
		error("getSIR.jl: Observation and reference signals should have the same length")
	end
	# --- Calculating SIR
	sirLinear	= mean( (abs2.( ur - dr))) / mean( (abs2.(dr)));
	# --- Setting output
	if castype == 1
		# --- Return linear SIR
		return sirLinear
	else
		# --- Return SIR in dB scale
		return 10*log10(sirLinear)
	end
end
export getSIR;

"""
---
Calculate the average power of the input signal
σ	= 1 / N Σ | x[n] | ^2
# --- Syntax
σ	= avgPower(x);
# --- Input parameters
- x	  : Input signal [Array{Any}]
# --- Output parameters
- σ	  : Estimated power [Float64]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function avgPower(x)
	return 1/length(x) * sum( abs2.(x) ) ;
end
export avgPower

## --- Complex Convolution definition
## In DSP conv is only defined as Re * Re -> Extend to C
#import DSP.conv
##function conv(x::Array{Complex{T}},h::Array{Complex{T}}) where T
	 ##y = conv(real(x),real(h)).-conv(imag(x),imag(h)).+ 1im*(conv(imag(x),real(h)).+conv(real(x),imag(h)));
##end
#function conv(x::Array{Complex{T}},h) where T
	#y = conv(real(x),h).+1im*conv(imag(x),h);
#end
#function conv(x,h::Array{Complex{T}}) where T
	#y = conv(x,real(h)).+1im*conv(x,imag(h));
#end

# ----------------------------------------------------
# --- Waveform type definition
# ----------------------------------------------------
"""
---
Abstract type gathering all waveform configuration
# v 1.0 - Robin Gerzaguet.
"""
abstract type Waveform end
export Waveform

# --- Waveforms functions
# ZC sequence
include("./Waveforms/genZCSequence.jl");
export genZCSequence;
# LTE Allocation
include("./Waveforms/getLTEAlloc.jl");
export getLTEAlloc;

# OFDM generation
include("./Waveforms/OFDM/ofdmSigGen.jl");
include("./Waveforms//OFDM/ofdmSigDecode.jl");
export ofdmSigGen, initOFDM;
export ofdmSigGen!;
export ofdmSigDecode;
export ofdmSigDecode!;

# SC-FDMA generation
include("./Waveforms//SCFDMA/scfdmaSigGen.jl");
include("./Waveforms//SCFDMA/scfdmaSigDecode.jl");
export scfdmaSigGen, initSCFDMA;
export scfdmaSigDecode, scfdmaPostProcessing;

# UF-OFDM generation
include("./Waveforms//UFOFDM/filterUFOFDM.jl")
include("./Waveforms//UFOFDM/ufofdmSigGen.jl");
include("./Waveforms//UFOFDM/ufofdmSigDecode.jl");
export ufofdmSigGen, initUFOFDM;
export ufofdmSigDecode;

# BF-OFDM generation
include("./Waveforms/BFOFDM/BFOFDM_filter.jl")
include("./Waveforms/BFOFDM/carrierManipulation.jl");
include("./Waveforms/BFOFDM/bfofdmSigGen.jl");
include("./Waveforms/BFOFDM/bfofdmSigDecode.jl");
export getBFOFDMFilter;
export getBFOFDM_subCarrierFromFBMCCarriers, getBFOFDM_carrierFromSubcarriers,getBFOFDM_oversampledGridSubcarriers;
export bfofdmSigGen, initBFOFDM;
export bfofdmSigDecode;

# WOLA generation
include("./Waveforms/WOLA/getWolaWindow.jl");
export getWolaWindow;
include("./Waveforms/WOLA/wolaSigGen.jl");
export initWOLA,wolaSigGen
include("./Waveforms/WOLA/wolaSigDecode.jl");
export wolaSigDecode

# FBMC Generation
include("./Waveforms/FBMC/fbmcSigGen.jl");
export fbmcSigGen
export getFBMCFilter
export oqamMapping
export initFBMC
include("./Waveforms/FBMC/fbmcSigDecode.jl");
export oqamDemapping
export fbmcSigDecode

# --- Global waveform alias
include("./Waveforms/genSig.jl");
export genSig
export getWaveformName;
export initWaveforms
export decodeSig

end


text:
struct Trans2BoolVariate <: AbstractBoolVariate
    name::String
    baseboolvar1::AbstractBoolVariate
    baseboolvar2::AbstractBoolVariate
    transform::Function
end

Base.length(var::Trans2BoolVariate) = length(var.baseboolvar1)

function slice(boolvar::Trans2BoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    baseboolvar1 = boolvar.baseboolvar1
    baseboolvar2 = boolvar.baseboolvar2
    f = boolvar.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip2(slice(baseboolvar1, fromobs, toobs, slicelength), slice(baseboolvar2, fromobs, toobs, slicelength))
    mapslice2(f, slices, slicelength, Bool)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
# Pkg.clone("https://github.com/bcbi/BCBI_base.jl.git")
using BCBI_base

Pkg.add("Cairo")
#it will fail due to a fontconfig issue in homebre see https://github.com/JuliaPackaging/Homebrew.jl/issues/186
#but just moving on seem ok
Pkg.build("Cairo")

const julia_pkgs = [ "VegaLite", "PlotlyJS", "Gadfly"]
BCBI_base.add(julia_pkgs)

const python_pkgs = [ "PyCall", "Pandas", "PyPlot", "Seaborn"]
BCBI_base.add(python_pkgs)

Pkg.add("IJulia")

const datasets_pkgs = ["RDatasets", "VegaDatasets"]
BCBI_base.add(datasets_pkgs)


text:
module metrics
export l1_loss

using Statistics: mean

function l1_loss(x::AbstractArray, y::AbstractArray)
    return mean(x-y)
end


text:
module ShaleDrillingLikelihood_DynamicDrillModelTest

DOBTIME = false
PRINTSTUFF = false
DOPROFILE = false

using ShaleDrillingLikelihood
using Test
using SparseArrays
using BenchmarkTools
using Calculus
using Random
using Profile
# using ProfileView
using InteractiveUtils

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: DCDPEmax,
    DCDPTmpVars,
    ubVfull    ,
    dubVfull   ,
    dubVfullperm,
    q          ,
    lse        ,
    tmp        ,
    tmp_cart   ,
    Πψtmp      ,
    IminusTEVp,
    ΠsumdubVj,
    dev0tmpj,
    ubV,
    dubV,
    dubVperm,
    flow!,
    end_ex0,
    _βΠψ!,
    _βΠψdθρ!,
    check_dΠψ,
    _dρdθρ,
    _dρdσ,
    theta_ρ,
    reward,
    peturb,
    vw_revenue,
    idx_ρ,
    revenue,
    _dψ1dρ,
    _dψ1dθρ,
    _ψ1,
    _ρ,
    dp1space,
    sprimes,
    _horizon,
    vfit!,
    tmpEVj,
    dEV,
    EV,
    anticipate_t1ev,
    _dmax,
    solve_inf_vfit!,
    solve_inf_pfit!,
    solve_inf_vfit_pfit!,
    gradinf!,
    solve_vf_terminal!,
    solve_vf_infill!,
    learningUpdate!,
    solve_vf_explore!,
    solve_vf_all!,
    solve_vf_all_timing!,
    solve_vf_infill_timing!,
    _nSexp,
    update_static_payoffs!,
    NoValueFunction,
    ValueFunction,
    ValueFunctionArrayOnly,
    value_function

println("print to keep from blowing up")

@testset "Dynamic Drilling Model" begin

    Random.seed!(1234)
    f = DrillReward(DrillingRevenue(Constrained(),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant(), ScrapValue_Constant())
    theta = randn(_nparm(f))
    # @show theta

    nψ, dmx, nz =  13, 3, 11

    # psi
    ψs = range(-4.5; stop=4.5, length=nψ)

    # z transition
    zs = ( range(-1.5; stop=1.5, length=nz), )
    nz = prod(length.(zs))
    ztrans = spdiagm(-1 => fill(1.0/3, nz-1), 0 => fill(1.0/3, nz), 1 => fill(1.0/3, nz-1) )
    ztrans[1,1] = ztrans[1,2] = ztrans[end,end-1] = ztrans[end,end] = 0.5

    # wp = LeasedProblemContsDrill(dmx,4,5,3,2)
    wp = LeasedProblem(dmx,4,5,3,2)

    # ichars
    ichar = (2.0, 0.25,)

    # ddm object
    ddm_no_t1ev   = DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, false)
    ddm_with_t1ev = DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, true)

    if false
        @code_warntype ValueFunction(f, 0.9, wp, zs, ztrans, ψs)
        @code_warntype DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, true, ValueFunction)
        @code_warntype DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, true, ValueFunctionArrayOnly)
        @code_warntype DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, true, NoValueFunction)
    end

    @testset "VF Structs" begin
        vfao = ValueFunctionArrayOnly(ddm_no_t1ev)
        vf = ValueFunction(ddm_no_t1ev)
        @test EV(vfao) == EV(vf)
        @test dEV(vfao) == dEV(vf)
    end

    tmpv = DCDPTmpVars(ddm_no_t1ev)
    @test tmpv isa ShaleDrillingLikelihood.DCDPTmpVarsArray
    @test view(tmpv, 1:2) isa ShaleDrillingLikelihood.DCDPTmpVarsView

    fd = similar(dubVfull(tmpv))
    ubvminus = zero(ubVfull(tmpv))
    ubvplus  = zero(ubVfull(tmpv))
    thetaminus = similar(theta)
    thetaplus = similar(theta)

    @testset "Fill per-period flow payoffs" begin

        ddm = ddm_with_t1ev
        tmpvminus = DCDPTmpVars(ubvminus, dubVfull(tmpv), dubVfullperm(tmpv), q(tmpv), lse(tmpv), tmp(tmpv), tmp_cart(tmpv), Πψtmp(tmpv), IminusTEVp(tmpv), tmpEVj(tmpv), ΠsumdubVj(tmpv), dev0tmpj(tmpv) )
        tmpvplus  = DCDPTmpVars(ubvplus,  dubVfull(tmpv), dubVfullperm(tmpv), q(tmpv), lse(tmpv), tmp(tmpv), tmp_cart(tmpv), Πψtmp(tmpv), IminusTEVp(tmpv), tmpEVj(tmpv), ΠsumdubVj(tmpv), dev0tmpj(tmpv) )

        for sidx in 1:length(wp)
            T = eltype(theta)
            idxd = dp1space(wp,sidx)

            for i in OneTo(length(theta))
                thetaminus .= theta
                thetaplus .= theta
                h = peturb(theta[i])
                thetaplus[i] += h
                thetaminus[i] -= h
                twoh = thetaplus[i] - thetaminus[i]
                @test maximum(abs.(thetaplus .- thetaminus)) ≈ twoh

                fill!(dubVfull(tmpv), 0)
                fill!(ubVfull(tmpvminus), 0)
                fill!(ubVfull(tmpvplus), 0)

                update_static_payoffs!(view(tmpvminus, idxd), ddm, thetaminus, sidx, ichar, false)
                update_static_payoffs!(view(tmpvplus , idxd), ddm, thetaplus, sidx, ichar, false)
                @test all(dubVfull(tmpv) .== 0)

                fd[i,:,:,:] .= (ubVfull(tmpvplus) .- ubVfull(tmpvminus))./ twoh
            end

            update_static_payoffs!(view(tmpv, idxd), ddm, theta, sidx, ichar, true)

            fdnosig = fd[1:end-1,:,:,:] .- dubVfull(tmpv)[1:end-1,:,:,:]
            fdsig = fd[end,:,:,:] .- dubVfull(tmpv)[end,:,:,:]

            if PRINTSTUFF
                @views maxv, idx = findmax(abs.(fdnosig))
                println("worst value is $maxv at $(CartesianIndices(fdnosig)[idx])")

                @views maxv, idx = findmax(abs.(fdsig))
                println("worst value is $maxv at $(CartesianIndices(fdsig)[idx])")
            end

            @test fd ≈ dubVfull(tmpv)
        end

        if DOBTIME
            @btime update_static_payoffs!($tmpv, $ddm, $theta, 2, (2.0, 0.25,), true)
        end
    end

    @testset "Learning transition" begin
        ddm = ddm_with_t1ev
        check_dΠψ(ddm, theta)
        _βΠψ!(tmpv, ddm, theta)
        _βΠψdθρ!(tmpv, ddm, theta)

        @testset "learning transition small functions" begin
            for x0 ∈ (0.25, 0.5, 0.75,)
                vspace = range(-2.0, stop=2.0, length=5)
                @test Calculus.derivative(_ρ, x0) ≈ _dρdθρ(x0)

                for u in vspace
                    for v in vspace
                        @test Calculus.derivative((z) -> _ψ1(u,v,_ρ(z)), x0) ≈ _dψ1dθρ(u,v,_ρ(x0), x0)
                        @test Calculus.derivative((z) -> _ψ1(u,v,z), x0) ≈ _dψ1dρ(u,v,x0)
                    end
                end
            end
        end

        @testset "Frechet derivative of ψ transition matrix" begin

            ddm = ddm_with_t1ev
            rwrd = reward(ddm)
            σv = theta_ρ(rwrd, theta)
            h = peturb(σv)
            σ1 = σv - h
            σ2 = σv + h
            hh = σ2 - σ1

            theta1 = copy(theta)
            theta2 = copy(theta)
            setindex!(vw_revenue(rwrd, theta1), σ1, idx_ρ(revenue(rwrd)))
            setindex!(vw_revenue(rwrd, theta2), σ2, idx_ρ(revenue(rwrd)))
            hh = theta_ρ(reward(ddm), theta2) - theta_ρ(reward(ddm), theta1)
            @test hh > 0

            pipsi = Πψtmp(tmpv)
            fdpi = similar(pipsi)
            @test fdpi isa Matrix
            @test !(fdpi === pipsi)
            pi1 = similar(fdpi)
            pi2 = similar(fdpi)

            _βΠψ!(tmpv, ddm, theta1)
            pi1 .= pipsi
            _βΠψ!(tmpv, ddm, theta2)
            pi2 .= pipsi

            fdpi .= (pi2 .- pi1) ./ hh
            @test any(fdpi .!= 0)
            _βΠψdθρ!(tmpv, ddm, theta)

            @views maxv, idx = findmax(abs.(pipsi .- fdpi))
            sub = CartesianIndices(fdpi)[idx]
            PRINTSTUFF && println("worst value is $maxv at $sub for β dΠ/dσ")

            @test 0.0 < maxv < 1e-5
            @test fdpi ≈ pipsi

        end
    end

    @testset "VF interpolation" begin

    end

    @testset "VF Iteration" begin

        thetaminus = similar(theta)
        thetaplus = similar(theta)

        for ddm in (ddm_with_t1ev, ddm_no_t1ev)

            evs = ValueFunction(ddm)
            fdEV = zero(dEV(evs))

            # println("\n-----------------------\nanticipate_e = $(anticipate_t1ev(ddm))\n---------------------------------------")

            @testset "Finite horizon VF gradient with anticipate_e = $(anticipate_t1ev(ddm))" begin
                for i in 1:length(statespace(ddm))
                    PRINTSTUFF && println("state $i of $(length(statespace(ddm)))")
                    idxd  = dp1space(wp,i)
                    idxs  = collect(sprimes(statespace(ddm),i))
                    horzn = _horizon(wp,i)

                    tmp_vw = view(tmpv, idxd)

                    EV0   = view(EV(evs) ,    :, :, i)
                    EV1   = view(EV(evs) ,    :, :, idxs)
                    dEV0  = view(dEV(evs), :, :, :, i)
                    dEV1  = view(dEV(evs), :, :, :, idxs)
                    fdEV0 = view(fdEV    , :, :, :, i)

                    fill!(fdEV0, 0)

                    for k in OneTo(length(theta))
                        thetaminus .= theta
                        thetaplus .= theta
                        h = peturb(theta[k])
                        thetaminus[k] -= h
                        thetaplus[k] += h
                        hh = thetaplus[k] - thetaminus[k]

                        fill!(evs, 0)
                        update_static_payoffs!(tmp_vw, ddm, thetaminus, i, ichar, false)
                        vfit!(EV0, tmp_vw, ddm)
                        ubV(tmp_vw) .+= discount(ddm) .* EV1
                        fdEV0[:,:,k] .-= EV0

                        fill!(evs, 0)
                        update_static_payoffs!(tmp_vw, ddm, thetaplus, i, ichar, false)
                        vfit!(EV0, tmp_vw, ddm)
                        ubV(tmp_vw) .+= discount(ddm) .* EV1
                        fdEV0[:,:,k] .+= EV0

                        fdEV0[:,:,k] ./= hh
                    end

                    fill!(evs, 0)
                    update_static_payoffs!(tmp_vw, ddm, theta, i, ichar, true)
                    ubV(tmp_vw)  .+= discount(ddm) .* EV1
                    dubVperm(tmp_vw) .+= discount(ddm) .* dEV1
                    vfit!(EV0, dEV0, tmp_vw, ddm)

                    PRINTSTUFF && println("\textrema(dEV0) = $(extrema(dEV0))")
                    PRINTSTUFF && println("\textrema(fdEV0) = $(extrema(fdEV0))")

                    @views maxv, idx = findmax(abs.(fdEV0 .- dEV0))
                    @views sub = CartesianIndices(fdEV0)[idx]
                    PRINTSTUFF && println("\tworst value is $maxv at $sub for dθ")

                    @test 0.0 <= maxv < 1.0
                    @test all(isfinite.(dEV0))
                    @test all(isfinite.(fdEV0))
                    @test fdEV0 ≈ dEV0
                end
            end # test finite horizon

            @testset "Check vfit/pfit for infinite horizon, anticipate_e = $(anticipate_t1ev(ddm))" begin
                i = length(statespace(ddm))-1

                fdEV = zero(dEV(evs))

                idxd  = dp1space(wp,i)
                idxs  = collect(sprimes(statespace(ddm),i))
                horzn = _horizon(wp,i)
                @test horzn == :Infinite

                EV0   = view(EV(evs) ,    :, :, i)
                EV1   = view(EV(evs) ,    :, :, idxs)
                dEV0  = view(dEV(evs), :, :, :, i)
                dEV1  = view(dEV(evs), :, :, :, idxs)
                fdEV0 = view(fdEV    , :, :, :, i)

                tmp_vw = view(tmpv, idxd)

                vfEV0 = zeros(size(EV0))

                fill!(evs, 0)
                update_static_payoffs!(tmp_vw, ddm, theta, i, ichar, true)
                ubV(tmp_vw)  .+= discount(ddm) .* EV1
                dubVperm(tmp_vw) .+= discount(ddm) .* dEV1
                converged, iterv, bnds = solve_inf_vfit!(EV0, tmp_vw, ddm; maxit=1000, vftol=1e-10)
                PRINTSTUFF && println("vfit done. converged = $converged after $iterv iterations. error bounds are $bnds")
                vfEV0 .= EV0

                fill!(evs, 0)
                update_static_payoffs!(tmp_vw, ddm, theta, i, ichar, true)
                ubV(tmp_vw)  .+= discount(ddm) .* EV1
                dubVperm(tmp_vw) .+= discount(ddm) .* dEV1
                converged, iterp, bnds = solve_inf_vfit_pfit!(EV0, tmp_vw, ddm; vftol=1e-10, maxit0=20, maxit1=40)
                PRINTSTUFF && println("pfit done. converged = $converged after $iterp iterations. error bounds are $bnds")

                @test iterv > iterp

                @views maxv, idx = findmax(abs.(vfEV0 .- EV0))
                @views sub = CartesianIndices(EV0)[idx]
                PRINTSTUFF &&  println("worst value is $maxv at $sub for vfit vs pfit")
                @test EV0 ≈ vfEV0

            end # test infinite horizon

            @testset "Check gradient for infinite horizon, anticipate_e = $(anticipate_t1ev(ddm))" begin

                i = length(statespace(ddm))-1

                fdEV = zero(dEV(evs))

                idxd  = dp1space(wp,i)
                idxs  = collect(sprimes(statespace(ddm),i))
                horzn = _horizon(wp,i)
                @test horzn == :Infinite

                EV0   = view(EV(evs) ,    :, :, i)
                EV1   = view(EV(evs) ,    :, :, idxs)
                dEV0  = view(dEV(evs), :, :, :, i)
                dEV1  = view(dEV(evs), :, :, :, idxs)
                fdEV0 = view(fdEV    , :, :, :, i)

                tmp_vw = view(tmpv, idxd)

                vfEV0 = zeros(size(EV0))

                nk = length(theta)
                for k in OneTo(nk)
                    thetaminus .= theta
                    thetaplus .= theta
                    h = peturb(theta[k])
                    thetaminus[k] -= h
                    thetaplus[k] += h
                    hh = thetaplus[k] - thetaminus[k]

                    fill!(evs, 0)
                    update_static_payoffs!(tmp_vw, ddm, thetaminus, i, ichar, false)
                    ubV(tmp_vw) .+= discount(ddm) .* EV1
                    solve_inf_vfit_pfit!(EV0, tmp_vw, ddm; vftol=1e-11, maxit0=10, maxit1=40)
                    fdEV0[:,:,k] .-= EV0

                    fill!(evs, 0)
                    update_static_payoffs!(tmp_vw, ddm, thetaplus, i, ichar, false)
                    ubV(tmp_vw) .+= discount(ddm) .* EV1
                    solve_inf_vfit_pfit!(EV0, tmp_vw, ddm; vftol=1e-11, maxit0=10, maxit1=40)
                    fdEV0[:,:,k] .+= EV0

                    fdEV0[:,:,k] ./= hh
                end
                @test any(fdEV0 .!= 0)

                fill!(evs, 0.0)
                update_static_payoffs!(tmp_vw, ddm, theta, i, ichar, true)
                ubV(tmp_vw)      .+= discount(ddm) .* EV1
                dubVperm(tmp_vw) .+= discount(ddm) .* dEV1

                solve_inf_vfit_pfit!(EV0, tmp_vw, ddm; vftol=1e-11, maxit0=10, maxit1=40)
                ubV(tmp_vw)[:,:,1] .= discount(ddm) .* EV0
                gradinf!(dEV0, tmp_vw, ddm)   # note: destroys ubV & dubV
                @test any(dEV0 .!= 0)

                PRINTSTUFF && println("extrema(dEV0) = $(extrema(dEV0))")
                PRINTSTUFF && println("extrema(fdEV0) = $(extrema(fdEV0))")
                @views maxv, idx = findmax(abs.(fdEV0 .- dEV0))
                @views sub = CartesianIndices(fdEV0)[idx]
                PRINTSTUFF && println("worst value is $maxv at $sub for dθ")

                @test 0.0 < maxv < 1.0
                @test all(isfinite.(dEV0))
                @test all(isfinite.(fdEV0))
                @test fdEV0 ≈ dEV0

            end

            # ---------------- Regime 1 VFI and PFI - test gradient ------------------

            @testset "Gradient of EV for infill, anticipate_e = $(anticipate_t1ev(ddm))" begin
                θ1 = similar(theta)
                θ2 = similar(theta)
                fdEV = similar(dEV(evs))
                itype = (4.7, 0.25,)

                nk = length(theta)

                for k = 1:nk
                    h = peturb(theta[k])
                    θ1 .= theta
                    θ2 .= theta
                    θ1[k] -= h
                    θ2[k] += h
                    hh = θ2[k] - θ1[k]

                    fill!(evs, 0)
                    solve_vf_terminal!(ddm, false)
                    @test all(EV(evs) .== 0)
                    @test all(dEV(evs) .== 0)
                    solve_vf_infill!(tmpv, ddm, θ1, itype, false)
                    fdEV[:,:,k,:] .= -EV(evs)

                    @test any( EV(evs) .!= 0)
                    @test all(dEV(evs) .== 0)

                    solve_vf_terminal!(ddm, false)
                    solve_vf_infill!(tmpv, ddm, θ2, itype, false)
                    fdEV[:,:,k,:] .+= EV(evs)
                    fdEV[:,:,k,:] ./= hh
                end

                solve_vf_terminal!(ddm, true)
                solve_vf_infill!(tmpv, ddm, theta, itype, true)
                @test any(dEV(evs) .!= 0)

                @views maxv, idx = findmax(abs.(fdEV[:,:,2:end,:].-dEV(evs)[:,:,2:end,:]))
                @views sub = CartesianIndices(fdEV[:,:,2:end,:])[idx]
                PRINTSTUFF && println("worst value is $maxv at $sub for dθ[2:end]")

                maxv, idx = findmax(abs.(fdEV.-dEV(evs)))
                sub = CartesianIndices(fdEV)[idx]
                PRINTSTUFF && println("worst value is $maxv at $sub")
                @test 0.0 < maxv < 1.0
                @test all(isfinite.(dEV(evs)))
                @test all(isfinite.(fdEV))
                @test fdEV ≈ dEV(evs)
            end

            @testset "Gradient of EV for Exploratory + learning, anticipate_e = $(anticipate_t1ev(ddm))" begin

                θ1 = similar(theta)
                θ2 = similar(theta)
                fdEV = similar(dEV(evs))
                itype = (4.7, 0.25,)

                nsexp = _nSexp(statespace(ddm))
                nk = length(theta)

                for k = 1:nk
                    h = peturb(theta[k])
                    θ1 .= theta
                    θ2 .= theta
                    θ1[k] -= h
                    θ2[k] += h
                    hh = θ2[k] - θ1[k]

                    fill!(evs, 0)

                    solve_vf_all!(tmpv, ddm, θ1, itype, false)
                    fdEV[:,:,k,:] .= -EV(evs)

                    solve_vf_all!(tmpv, ddm, θ2, itype, false)
                    fdEV[:,:,k,:] .+= EV(evs)

                    fdEV[:,:,k,:] ./= hh
                end

                # ----------------- analytic -----------------

                solve_vf_all!(tmpv, ddm, theta, itype, true)

                # check infill + learning portion gradient
                @views a =     fdEV[:,:, :, nsexp+1:end]
                @views b = dEV(evs)[:,:, :, nsexp+1:end]
                maxv, idx = findmax(abs.(a .- b))
                PRINTSTUFF && println("worst value is $maxv at $(CartesianIndices(a)[idx]) for infill")

                # check exploration portion of gradient
                @views a =     fdEV[:,:, 1:end-1, 1:nsexp]
                @views b = dEV(evs)[:,:, 1:end-1, 1:nsexp]
                maxv, idx = findmax(abs.(a .- b))
                PRINTSTUFF && println("worst value is $maxv at $(CartesianIndices(a)[idx]) for infill")

                # check exploration portion of gradient
                @views a =     fdEV[:,:, end, 1:nsexp]
                @views b = dEV(evs)[:,:, end, 1:nsexp]
                maxv, idx = findmax(abs.(a .- b))
                PRINTSTUFF && println("worst value is $maxv at $(CartesianIndices(a)[idx]) for infill")

                @test 0.0 < maxv < 0.1
                @test all(isfinite.(dEV(evs)))
                @test all(isfinite.(fdEV))

                @test fdEV ≈ dEV(evs)
                # println("dEV/dθ looks ok! :)")
            end

            if DOBTIME
                itype = (4.7, 0.25,)
                println("Timing solve_vf_all! with anticipate_e = $(anticipate_t1ev(ddm))")
                fill!(value_function(ddm), 0)
                @btime solve_vf_terminal!(            $ddm,                 true)
                @btime solve_vf_infill_timing!($tmpv, $ddm, $theta, $itype, true)
                @btime learningUpdate!(        $tmpv, $ddm, $theta, $itype, true)
                @btime solve_vf_explore!(      $tmpv, $ddm, $theta, $itype, true)

                @btime solve_vf_all_timing!($tmpv, $ddm, $theta, $itype, true)
            end

            if DOPROFILE
                println("Profiling solve_vf_all_timing! with anticipate_e = $(anticipate_t1ev(ddm))")
                itype = (4.7, 0.25,)

                @code_warntype solve_vf_terminal!(      ddm,               true)
                @code_warntype solve_vf_infill!(    tmpv, ddm, theta, itype, true)
                @code_warntype learningUpdate!(     tmpv, ddm, theta, itype, true)
                @code_warntype solve_vf_explore!(   tmpv, ddm, theta, itype, true)
                fill!(evs, 0)
                @code_warntype solve_vf_all!(       tmpv, ddm, theta, itype, true)
                @code_warntype solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                fill!(evs, 0)
                # solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # Profile.clear()
                # @profile solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # @profile solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # @profile solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # @profile solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # @profile solve_vf_all_timing!(tmpv, ddm, theta, itype, true)
                # Juno.profiletree()
                # Juno.profiler()
                # Profile.print(format=:flat)
                # ProfileView.view()
                # pprof()
            end

        end # ddm
    end # VF Iteration

end # testset
end # module