text:
# S&P500指数の値動きからSPXLを計算するためのデータを作成する
using CSV, DataFrames, Dates, Plots
pyplot(fmt=:svg, fontfamily="meiryo") # plotの設定

# 正解値SPXL読み込み
dfmt = dateformat"yyyy/mm/dd HH:MM:SS"
spxl = DataFrame(CSV.File("./data/SPXL_20210430.csv", dateformat=dfmt, header=true))
dt = spxl.timestamp .- Dates.Hour(14) # Yahooで取得をしたデータは日本時間なので，USに変換する
spxl.timestamp = dt

# S&P500のDividendYield(配当利回り)の読み込み--------------------------------
dfmt = dateformat"yyyy/mm/dd"
sp500_yield = DataFrame(CSV.File("./data/S&P500_DividendYield.csv", dateformat=dfmt, header=true))
sort!(sp500_yield) # 昇順にソートsort

# SPYのdividendsは四半期単位なので，サンプリングをSPXLに合わせて前値ホールドする
spy_dint = copy(spxl)
select!(spy_dint, Not(:open))
rename!(spy_dint, "low" => "yield") # その年の配当金利回り
rename!(spy_dint, "high" => "dtb3") # 米国3ヶ月債金利
rename!(spy_dint, "close" => "spxl") # SPXL価格
rename!(spy_dint, "volume" => "spxl_pred") # 予測SPXL

d_temp = 1
for y in 1:length(sp500_yield.yield)
    println("y = " * string(y))
    for d in d_temp:length(spy_dint.timestamp)
        # println("d = " * string(d))
        yield = sp500_yield[y,"yield"] # 配当利回りを取得
        spy_dint.yield[d] = yield # 次のyearに該当するまで，前値を当てはめる
        if (Dates.Date(spy_dint.timestamp[d]) > sp500_yield.date[y])
            d_temp = d
            break
        end
    end
end

# DTB3(金利)の読み込み--------------------------------
dfmt = dateformat"yyyy/mm/dd"
dtb3 = DataFrame(CSV.File("./data/DTB3.csv", dateformat=dfmt, header=true))
dfmt = dateformat"yyyy-mm-dd"
dt = [DateTime(x, dfmt) for x in dtb3[:,"DATE"]]
dtb3.timestamp = dt
select!(dtb3, Not(:DATE))

# 値を入れていく
for y in 1:length(spy_dint.timestamp)
    for d in 1:length(dtb3.timestamp)
        if (Date(spy_dint.timestamp[y]) == Date(dtb3.timestamp[d]))
            spy_dint.dtb3[y] = dtb3.DTB3[d]
            # println("d = " * string(d))
            break
        end
        spy_dint.dtb3[y] = NaN
    end
end
spy_dint = spy_dint[ .!(isnan.(spy_dint.dtb3)), :] # NaNのある行の削除

# S&P500の読み込み--------------------------------
dfmt = dateformat"yyyy/mm/dd HH:MM:SS"
sp500 = DataFrame(CSV.File("./data/S&P500.csv", dateformat=dfmt, header=true))
dt = sp500.timestamp .- Dates.Hour(14) # Yahooで取得をしたデータは日本時間なので，USに変換する
sp500.timestamp = dt
select!(sp500, Not(:open))
select!(sp500, Not(:low))
select!(sp500, Not(:high))
select!(sp500, Not(:volume))
rename!(sp500, "close" => "sp500") # SPY価格

spy_dint = innerjoin(spy_dint, sp500, on=:timestamp)

# 疑似SPXLを計算する元データを保存
CSV.write("./data/spxl_prediction_data.csv", spy_dint)

# 疑似SPXLを計算-----------------------------------------------
# 1. 単純計算が実際のSPXLより大きいことを確認
dif = (spy_dint.sp500[2:end] - spy_dint.sp500[1:end - 1]) # 1日の値動き幅
dif_rate = dif ./ spy_dint.sp500[1:end - 1] # 値動き割合
dif_rate_y = dif_rate + spy_dint.yield[1:end - 1] ./ 365 # 配当込み1日値動き割合
dif_rate3 = dif_rate_y * 3
spxl_pred1 = ones(length(spy_dint.sp500), 1) * 4.25
for i = 2:length(spxl_pred1)
    spxl_pred1[i] = spxl_pred1[i - 1] * dif_rate3[i - 1] + spxl_pred1[i - 1]
end

plot(spy_dint.timestamp, spy_dint.spxl,
    xlabel="年月日", ylabel="価格[\$]", label="SPXL")
plot!(spy_dint.timestamp, spxl_pred1, label="SPXL(S&P500値動き3倍)")
savefig("SPXL_dif_rate3.png")

# 2. 0.95% 考慮でも実際のSPXLより大きいことを確認
dif_rate3_c095 = dif_rate3 .- (0.95 / 100 / 365)
spxl_pred2 = ones(length(spy_dint.sp500), 1) * 4.25
for i = 2:length(spxl_pred2)
    spxl_pred2[i] = spxl_pred2[i - 1] * dif_rate3_c095[i - 1] + spxl_pred2[i - 1]
end

plot(spy_dint.timestamp, spy_dint.spxl,
    xlabel="年月日", ylabel="価格[\$]", label="SPXL")
plot!(spy_dint.timestamp, spxl_pred2, label="SPXL(S&P500値動き3倍+コスト0.95%)")
savefig("SPXL_dif_rate3_c095.png")

spxl_err2 = (spxl_pred2 .- spy_dint.spxl) ./ spy_dint.spxl
plot(spy_dint.timestamp, spxl_err2,
    xlabel="年月日", ylabel="SPXL予測の乖離率",
    ylims=(-0.2, 0.4))
savefig("SPXL_dif_rate3_c095_err.png")

# 3. 計算式に割り当てる
dif_rate3_dtb3 = dif_rate3 - ( spy_dint.dtb3[1:end - 1] / 100 / 365 * 2)
dif_rate3_dtb3_c215 = dif_rate3_dtb3 .- (2.15 / 100 / 365)
spxl_pred3 = ones(length(spy_dint.sp500), 1) * 4.25
for i = 2:length(spxl_pred3)
    spxl_pred3[i] = spxl_pred3[i - 1] * dif_rate3_dtb3_c215[i - 1] + spxl_pred3[i - 1]
end

plot(spy_dint.timestamp, spy_dint.spxl,
    xlabel="年月日", ylabel="価格[\$]", label="SPXL")
plot!(spy_dint.timestamp, spxl_pred3, label="SPXL(S&P500値動き3倍+先物金利+コスト2.15%)")
savefig("SPXL_dif_rate3_dtb3_c215.png")

spxl_err3 = (spxl_pred3 .- spy_dint.spxl) ./ spy_dint.spxl
plot(spy_dint.timestamp, spxl_err3,
    xlabel="年月日", ylabel="SPXL予測の乖離率",
    ylims=(-0.2, 0.4))
savefig("SPXL_dif_rate3_dtb3_c215_err.png")

# 3.2 固定コストを上げる
dif_rate3_dtb3_c315 = dif_rate3_dtb3 .- (3.15 / 100 / 365)
spxl_pred4 = ones(length(spy_dint.sp500), 1) * 4.25
for i = 2:length(spxl_pred4)
    spxl_pred4[i] = spxl_pred4[i - 1] * dif_rate3_dtb3_c315[i - 1] + spxl_pred4[i - 1]
end

plot(spy_dint.timestamp, spy_dint.spxl,
    xlabel="年月日", ylabel="価格[\$]", label="SPXL")
plot!(spy_dint.timestamp, spxl_pred4, label="SPXL(S&P500値動き3倍+先物金利+コスト3.15%)")
savefig("SPXL_dif_rate3_dtb3_c315.png")

spxl_err4 = (spxl_pred4 .- spy_dint.spxl) ./ spy_dint.spxl
plot(spy_dint.timestamp, spxl_err4,
    xlabel="年月日", ylabel="SPXL予測の乖離率",
    ylims=(-0.2, 0.4))
savefig("SPXL_dif_rate3_dtb3_c315_err.png")

# 1954以降の疑似SPXLとS&P 500を計算-----------------------------------------------
# SP500, SP500利回り，米国3ヶ月債金利を用意
spxl_pred = sp500[sp500.timestamp.>Dates.Date("1954-1-4"), :]
spxl_pred.timestamp = Dates.Date.(spxl_pred.timestamp)
dtb3.timestamp = Dates.Date.(dtb3.timestamp)
spxl_pred = innerjoin(spxl_pred, dtb3, on=:timestamp)

spxl_pred.sp500_yield = ones(length(spxl_pred.timestamp))
d_temp = 1
for y in 1:length(sp500_yield.yield)
    println("y = " * string(y))
    for d in d_temp:length(spxl_pred.timestamp)
        # println("d = " * string(d))
        yield = sp500_yield[y,"yield"] # 配当利回りを取得
        spxl_pred.sp500_yield[d] = yield # 次のyearに該当するまで，前値を当てはめる
        if (Dates.Date(spxl_pred.timestamp[d]) > sp500_yield.date[y])
            d_temp = d
            break
        end
    end
end

# 1954年のSP500と同一金額からスタートした場合のSPXLを計算
spxl_pred.spxl = ones(length(spxl_pred.timestamp))*spxl_pred.sp500[1]

dif = (spxl_pred.sp500[2:end] - spxl_pred.sp500[1:end - 1]) # 1日の値動き幅
dif_rate = dif ./ spxl_pred.sp500[1:end - 1] # 値動き割合
dif_rate_y = dif_rate + spxl_pred.sp500_yield[1:end - 1] ./ 365 # 配当込み1日値動き割合
dif_rate3 = dif_rate_y * 3 # 配当込値動き割合の3倍
dif_rate3_dtb3 = dif_rate3 - ( spxl_pred.DTB3[1:end - 1] / 100 / 365 * 2) # 3ヶ月債の2倍の金利を引く
dif_rate3_dtb3_c315 = dif_rate3_dtb3 .- (3.15 / 100 / 365) # 固定コスト3.15%を上乗せ

for i = 2:length(spxl_pred.timestamp)
    spxl_pred.spxl[i] = spxl_pred.spxl[i - 1] * dif_rate3_dtb3_c315[i - 1] + spxl_pred.spxl[i - 1]
end

# 比較対象としてSP500の配当込み指数を計算．
spxl_pred.sp500_yidx = ones(length(spxl_pred.timestamp))*spxl_pred.sp500[1]
dif_rate_y_c = dif_rate_y .- (0.141 / 100 / 365) # 固定コスト0.141%を上乗せ

for i = 2:length(spxl_pred.timestamp)
    spxl_pred.sp500_yidx[i] = spxl_pred.sp500_yidx[i - 1] * dif_rate_y_c[i - 1] + spxl_pred.sp500_yidx[i - 1]
end

CSV.write("spxl_pred_1954.csv", spxl_pred)

# トータルリターンの計算
sp500_ret = spxl_pred.sp500_yidx[spxl_pred.timestamp.==Dates.Date("2020-7-1")]/spxl_pred.sp500_yidx[spxl_pred.timestamp.==Dates.Date("1954-7-1")]
spxl_ret = spxl_pred.spxl[spxl_pred.timestamp.==Dates.Date("2020-7-1")]/spxl_pred.spxl[spxl_pred.timestamp.==Dates.Date("1954-7-1")]

plot(spxl_pred.timestamp, spxl_pred.sp500_yidx,
    xlabel="年月日", ylabel="価格[\$]", label="仮想インデックス(S&P500配当込み指数+コスト0.141%)", yscale=:log10)
plot!(spxl_pred.timestamp, spxl_pred.spxl, label="仮想SPXL(S&P500値動き3倍+先物金利+コスト3.15%)", yscale=:log10)
savefig("SP500_SPXL_1954.png")


text:
using SpinWaveTheory
using Documenter

DocMeta.setdocmeta!(SpinWaveTheory, :DocTestSetup, :(using SpinWaveTheory); recursive=true)

makedocs(;
    modules=[SpinWaveTheory],
    authors="waltergu <waltergu1989@gmail.com> and contributors",
    repo="https://github.com/Quantum-Many-Body/SpinWaveTheory.jl/blob/{commit}{path}#{line}",
    sitename="SpinWaveTheory.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Quantum-Many-Body.github.io/SpinWaveTheory.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "examples/Introduction.md",
            "examples/SquareLatticeFerromagnet.md",
            "examples/SquareLatticeAntiFerromagnet.md",
            "examples/CrBr3.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/Quantum-Many-Body/SpinWaveTheory.jl",
)

ext: jl


text:
# model.jl

"""
    abstract type AbstractModel end

Implements a statistical model. Has so far only SDEModel as subtype.
"""
abstract type AbstractModel end

"""
    type SDEModel{T<:Union{Float64, <:Function}} <: AbstractModel
    σ::T
    beginvalue::Float64
    end

Implements an SDE model with begin value beginvalue and σ.

# Warning
It is assumed that
for every b under consideration, the laws of dX_t=b(X_t)dt+σ(X_t)dW_t are
equivalent.

# Example

```julia
SDEModel(1.0, 0.0)
# Implements the model dX_t=b(X_t)dt+σ(X_t)dW_t, with X_0=0.0.
#
SDEModel(identity, 0.0)
# Implements the model dX_t=b(X_t)+X_tdW_t, with X_0=0.0.
```
"""
type SDEModel{T<:Union{Float64, <:Function}} <: AbstractModel
    σ::T
    beginvalue::Float64
end


text:
EXAMPLES_DIR = "../examples/"

function run_deepening_mixed_layer_example(arch)
    example_filepath = joinpath(EXAMPLES_DIR, "deepening_mixed_layer.jl")
    txt = read(example_filepath, String)

    arch == GPU() && (txt = replace(txt, "arch = CPU()" => "arch = GPU()"))

    txt = replace(txt, "N = 32" => "N = 16")
    txt = replace(txt, "model.clock.time < tf" => "model.clock.time < 0.5")
    txt = replace(txt, "time_step!(model, 10, wizard.Δt)" => "time_step!(model, 1, wizard.Δt)")

    test_script_filepath = "deepening_mixed_layer_example_cpu_test.jl"
    open(test_script_filepath, "w") do f
        write(f, txt)
    end

    try
        include(test_script_filepath)
    catch err
        @error sprint(showerror, err)
        rm(test_script_filepath)
        return false
    end

    rm(test_script_filepath)
    return true
end

@testset "Examples" begin
    println("Testing examples...")

    for arch in archs
        @testset "Deepening mixed layer example [$(typeof(arch))]" begin
            println("  Testing deepening mixed layer example [$(typeof(arch))]")
            @test run_deepening_mixed_layer_example(arch)
        end
    end
end


text:
export TabularRandomPolicy

"""
    TabularRandomPolicy(;table=Dict{Int, Float32}(), rng=Random.GLOBAL_RNG)

Use a `Dict` to store action distribution.
"""
Base.@kwdef struct TabularRandomPolicy{S,T,R} <: AbstractPolicy
    table::Dict{S,T} = Dict{Any,Vector{Float32}}()
    rng::R = Random.GLOBAL_RNG
end

TabularRandomPolicy{S}(; rng = Random.GLOBAL_RNG) where {S} =
    TabularRandomPolicy{S,Vector{Float32}}(; rng = rng)
TabularRandomPolicy{S,T}(; rng = Random.GLOBAL_RNG) where {S,T} =
    TabularRandomPolicy(Dict{S,T}(), rng)

RLBase.prob(p::TabularRandomPolicy, env::AbstractEnv) = prob(p, ChanceStyle(env), env)

function RLBase.prob(p::TabularRandomPolicy, ::ExplicitStochastic, env::AbstractEnv)
    if current_player(env) == chance_player(env)
        prob(env)
    else
        prob(p, DETERMINISTIC, env)  # treat it just like a normal one
    end
end

function RLBase.prob(t::TabularRandomPolicy, ::RLBase.AbstractChanceStyle, env::AbstractEnv)
    prob(t, ActionStyle(env), env)
end

function RLBase.prob(t::TabularRandomPolicy, ::FullActionSet, env::AbstractEnv)
    get!(t.table, state(env)) do
        m = legal_action_space_mask(env)
        m ./ sum(m)
    end
end

function RLBase.prob(t::TabularRandomPolicy, ::MinimalActionSet, env::AbstractEnv)
    get!(t.table, state(env)) do
        n = length(action_space(env))
        fill(1 / n, n)
    end
end

function RLBase.prob(t::TabularRandomPolicy, env::AbstractEnv, action)
    prob(t, env, action_space(env), action)
end

function RLBase.prob(t::TabularRandomPolicy, env::AbstractEnv, action_space, action)
    prob(t, env)[findfirst(==(action), action_space)]
end

function RLBase.prob(
    t::TabularRandomPolicy,
    env::AbstractEnv,
    action_space::Base.OneTo,
    action,
)
    prob(t, env)[action]
end

# function RLBase.prob(t::TabularRandomPolicy, env::AbstractEnv, action_space::ZeroTo, action)
#     prob(t, env)[action+1]
# end

function RLBase.prob(t::TabularRandomPolicy, state, action)
    # assume table is already initialized
    t.table[state][action]
end

(p::TabularRandomPolicy)(env::AbstractEnv) =
    sample(p.rng, action_space(env), Weights(prob(p, env), 1.0))

# !!! Assumeing table is already initialized
(p::TabularRandomPolicy{S})(state::S) where {S} =
    sample(p.rng, Weights(p.table[state], 1.0))

"""
    update!(p::TabularRandomPolicy, state => value)

You should manually check `value` sum to `1.0`.
"""
function RLBase.update!(p::TabularRandomPolicy, experience::Pair)
    s, dist = experience
    if haskey(p.table, s)
        p.table[s] .= dist
    else
        p.table[s] = dist
    end
end


text:
using Pkg
Pkg.activate(".")

module EnzymeMFABC
	include("example.jl")
end
using Plots, StatsBase, Statistics, LaTeXStrings
using FileIO, JLD2

ABC_stats = load("ABC_stats.jld2")
MFABC_output = load("MFABC_output.jld2")
BSL_stats = load("BSL_stats.jld2")
MFBSL_output = load("MFBSL_output.jld2")

c_abc = ABC_stats["c"]
v_abc = ABC_stats["v"]

c_mfabc = MFABC_output["c"]
v_mfabc = MFABC_output["v"]
w_mfabc = MFABC_output["w_mf"]
μ_mfabc = MFABC_output["μ"]
μ_t_mfabc = MFABC_output["μ_t"]
ν_t_mfabc = MFABC_output["ν_t"]

c_bsl = BSL_stats["c"]
v_bsl = BSL_stats["v"]

c_mfbsl = MFBSL_output["c"]
v_mfbsl = MFBSL_output["v"]
w_mfbsl = MFBSL_output["w_mf"]
μ_mfbsl = MFBSL_output["μ"]
μ_t_mfbsl = MFBSL_output["μ_t"]
ν_t_mfbsl = MFBSL_output["ν_t"]

function big_fig(c, c_mf, v, v_mf, w_mf_sparse, μ, μ_t, ν_t; alg::String, iterticks=:auto, kwargs...)

    # Comparison
    compare = plot(;
        xscale=:log10,
        xticks= 10.0 .^ (1:0.2:3),
        yscale=:log10,
        yticks= 10.0 .^ (-6:0.2:-3.0),
        xlabel="Simulation cost (s)",
        ylabel="Variance",
        title="(a) Performance scaling: "*alg*" comparison",
        fontfamily="Computer Modern",
        guidefontsize=10,
        tickfontsize=10,
        legendfontsize=10,
        titlefontsize=12,
        dpi=400,
        size=(450,300),
    )
    plot!(compare, c, v, label=alg, markershape=:square)
    plot!(compare, c_mf, v_mf, label="MF-"*alg, markershape=:star)

    wmffig = scatter(selectdim(w_mf_sparse, 2, 1), selectdim(w_mf_sparse, 2, 2),
            title="(b) Multifidelity weights",
            ylabel=L"w_{mf}",
            xlabel="Iteration",
            markershape=:circle,
            markercolor=:black,
            markersize=3,
            markeralpha=0.6,
            legend=:none,
            xticks=iterticks,
            fontfamily="Computer Modern",
            guidefontsize=10,
            tickfontsize=10,
            titlefontsize=12,
            dpi=400,
            size=(450,300),
    )

    D_idx = 100 .* (1:size(μ_t,1))
    mufig = plot(D_idx, μ_t,
        title="(c) Mean function values",
        ylabel=L"\log(\nu_k)",
        xlabel="Iteration",
        legend=:none,
        xticks=iterticks,
        fontfamily="Computer Modern",
        guidefontsize=10,
        tickfontsize=10,
        titlefontsize=12,
        dpi=400,
        size=(450,300),
    )

    compare_nu(j) = hcat(μ_t[:, j], log.(ν_t[:, j]))
    nunu = compare_nu(1)
    adaptfig = plot(D_idx, nunu;
        label=["Adaptive mean" "Optimal mean estimate"],
        title="(d) Adaptive mean function",
        xlabel="Iteration",
        ylabel=L"\log(\nu_1)",
        xticks=iterticks,
        fontfamily="Computer Modern",
        guidefontsize=10,
        tickfontsize=10,
        legendfontsize=10,
        titlefontsize=12,
        dpi=400,
        size=(450,300),
    )

    f = open("mu_"*alg*".tex", "w")
	EnzymeMFABC.print_tree(f, μ)
	close(f)

    figout = plot(
        compare,
        wmffig,
        mufig,
        adaptfig;
        layout=(2,2),
        size=(900,600),
        bottom_margin=6*Plots.mm,
        kwargs...
    )
    return figout
end

figABC = big_fig(c_abc, c_mfabc, v_abc, v_mfabc, w_mfabc, μ_mfabc, μ_t_mfabc, ν_t_mfabc; alg="ABC", iterticks=[200_000, 400_000, 600_000])
figBSL = big_fig(c_bsl, c_mfbsl, v_bsl, v_mfbsl, w_mfbsl, μ_mfbsl, μ_t_mfbsl, ν_t_mfbsl; alg="BSL", iterticks=[20_000, 40_000, 60_000])

traj_data = load("eg_traj.jld2")
trajfig = scatter(
    EnzymeMFABC.DATA,
    10:10:100;
    yticks=0:20:100,
    yminorticks=2,
    yminorgrid=true,
    xticks=:auto,
    c=:black,
    label="Data",
    legend=:topleft,
    title="(a) Simulated product accumulation",
)

plot!(trajfig,
    traj_data["t_lo"],
    selectdim(traj_data["x_lo"], 1, 2),
    c=1,
    label="Lo",
    seriestype=:steppost,
    fontfamily="Computer Modern",
    xlabel="Time",
    ylabel="Product molecule count",
    guidefontsize=10,
    tickfontsize=10,
    legendfontsize=10,
    titlefontsize=12,
    dpi=400,
    size=(450,300),
)

R = traj_data["R"]
uncoupled = traj_data["uncoupled"]
coupled = traj_data["coupled"]
for r in 1:R
    t_hi_uc, x_hi_uc = uncoupled[r]
    t_hi_cp, x_hi_cp = coupled[r]

    label_uc = r==1 ? "Hi (uncoupled)" : ""
    label_cp = r==1 ? "Hi (coupled)" : ""
    plot!(trajfig, t_hi_uc, selectdim(x_hi_uc, 1, 4), alpha=0.5, c=2, label=label_uc)
    plot!(trajfig, t_hi_cp, selectdim(x_hi_cp, 1, 4), alpha=0.5, c=3, label=label_cp)
end
trajfig

function get_summary(t, x)
    idx = [findfirst(==(p), x[end,:]) for p in 10:10:100]
    return t[idx]
end

offset=0.15
sumfig = scatter(1:10, EnzymeMFABC.DATA, c=:black, legend=:topleft, label="Data")
scatter!(sumfig, -offset .+ (1:10), get_summary(traj_data["t_lo"], traj_data["x_lo"]), c=1, label="Lo")
for r in 1:R
    _cp = coupled[r]
    _uc = uncoupled[r]
    label_cp = r==1 ? "Hi (coupled)" : ""
    label_uc = r==1 ? "Hi (uncoupled)" : ""
    scatter!(sumfig, offset .+ (1:10), get_summary(_uc...), c=2, label=label_uc)
    scatter!(sumfig, -2*offset .+ (1:10), get_summary(_cp...), c=3, label=label_cp)
end
plot!(sumfig;
 xticks=(1:10, [latexstring("y_{$(i)}") for i in 1:10]),
 yticks=EnzymeMFABC.DATA,
 title = "(b) Summarised trajectories",
 fontfamily="Computer Modern",
 xlabel="Data dimension",
 ylabel="Value",
 guidefontsize=10,
 tickfontsize=10,
 legendfontsize=10,
 titlefontsize=12,
 dpi=400,
 size=(450,300),
)

datafig = plot(trajfig, sumfig, layout=(1,2), size=(900,300), bottom_margin=5*Plots.mm, left_margin=5*Plots.mm)

savefig(figABC, "figABC.pdf")
savefig(figBSL, "figBSL.pdf")
savefig(datafig, "trajfig.pdf")

max_stars_repo_name: tpprescott/mf-lf

max_issues_repo_name: tpprescott/mf-lf

max_forks_repo_name: tpprescott/mf-lf


text:
# A simple example reading a StochOptFormat file and solving it via Benders
# decomposition.
#
# The code is demonstration only, and does not contain robust checks or many
# nice error messages, preferring to throw assertion errors.

using JuMP
import Clp
import JSON
import JSONSchema

"""
    mathoptformat_to_jump(node)

Convert a MathOptFormat model in `node["subproblem"]` into a JuMP equivalent
using the Clp optimizer.
"""
function mathoptformat_to_jump(node)
    # Read the problem into a MOI model first.
    model = MOI.FileFormats.Model(format = MOI.FileFormats.FORMAT_MOF)
    io = IOBuffer()
    write(io, JSON.json(node["subproblem"]))
    seekstart(io)
    MOI.read!(io, model)
    # Then copy it to a JuMP model.
    subproblem = Model(Clp.Optimizer)
    MOI.copy_to(subproblem, model)
    set_silent(subproblem)
    node["subproblem"] = subproblem
    return node
end

"""
    solve_first_stage(node)

Solve the first stage problem.
"""
function solve_first_stage(node)
    sp = node["subproblem"]
    optimize!(sp)
    if termination_status(sp) != MOI.OPTIMAL
        error("Unable to solve first stage to optimality!")
    end
    return (
        obj = objective_value(sp),
        x = Dict(
            name => value(variable_by_name(sp, s["out"]))
            for (name, s) in node["state_variables"]
        ),
        sol = Dict(name(x) => value(x) for x in all_variables(sp))
    )
end

"""
    solve_second_stage(node, state, noise)
"""
function solve_second_stage(node, state, noise)
    sp = node["subproblem"]
    for (name, s) in node["state_variables"]
        fix(variable_by_name(sp, s["in"]), state[name]; force = true)
    end
    for (name, w) in noise["support"]
        fix(variable_by_name(sp, name), w; force = true)
    end
    optimize!(sp)
    return Dict(
        "probability" => noise["probability"],
        "objective" => objective_value(sp),
        "pi" => Dict(
            name => shadow_price(FixRef(variable_by_name(sp, s["in"])))
            for (name, s) in node["state_variables"]
        )
    )
end

function add_cut(first_stage, x, ret_second)
    sp = first_stage["subproblem"]
    cut_term = @expression(
        sp,
        sum(
            ret["probability"] * (
                ret["objective"] +
                sum(
                    ret["pi"][name] * (variable_by_name(sp, s["out"]) - x[name])
                    for (name, s) in first_stage["state_variables"]
                )
            ) for ret in ret_second
        )
    )
    if objective_sense(sp) == MOI.MAX_SENSE
        @constraint(sp, first_stage["theta"] <= cut_term)
    else
        @constraint(sp, first_stage["theta"] >= cut_term)
    end
    return
end

function load_two_stage_problem(filename)
    data = JSON.parsefile(filename)
    @assert(data["version"]["major"] == 0)
    @assert(data["version"]["minor"] == 1)
    @assert(length(data["nodes"]) == 2)
    @assert(length(data["edges"]) == 2)
    nodes = Dict(
        name => mathoptformat_to_jump(node) for (name, node) in data["nodes"]
    )
    first_stage, second_stage = nothing, nothing
    for edge in data["edges"]
        if edge["from"] == data["root"]["name"]
            first_stage = nodes[edge["to"]]
        else
            second_stage = nodes[edge["to"]]
        end
    end
    sp = first_stage["subproblem"]
    for (name, init) in data["root"]["state_variables"]
        x = variable_by_name(sp, first_stage["state_variables"][name]["in"])
        fix(x, init["initial_value"]; force = true)
    end
    first_stage["theta"] = @variable(sp, -1e6 <= theta <= 1e6)
    set_objective_function(sp, objective_function(sp) + first_stage["theta"])
    return first_stage, second_stage
end

function benders(first_stage, second_stage, iteration_limit = 20)
    bounds = Tuple{Float64, Float64}[]
    # Collect the names of the outgoing state variables in the first stage.
    x_out = Set(s["out"] for s in values(first_stage["state_variables"]))
    for iter = 1:iteration_limit
        # Forward pass. Solve the first stage problem.
        ret_first = solve_first_stage(first_stage)
        ret_second = [
            solve_second_stage(second_stage, ret_first.x, noise)
            for noise in second_stage["realizations"]
        ]
        stat_bound = ret_first.obj - value(first_stage["theta"]) +
            sum(ret["probability"] * ret["objective"] for ret in ret_second)
        add_cut(first_stage, ret_first.x, ret_second)
        push!(bounds, (ret_first.obj, stat_bound))
        if abs(ret_first.obj - stat_bound) < 1e-6
            break
        end
    end
    return bounds, solve_first_stage(first_stage)
end

function validate(filename; schema_filename = "../sof.schema.json")
    schema = JSONSchema.Schema(JSON.parsefile(schema_filename))
    return JSONSchema.validate(JSON.parsefile(filename), schema)
end

validate("news_vendor.sof.json")
first_stage, second_stage = load_two_stage_problem("news_vendor.sof.json")
ret = benders(first_stage, second_stage)

# Check solution!
ret = solve_first_stage(first_stage)
@show ret
@assert(ret.x["x"] ≈ 10)

max_stars_repo_name: joaquimg/StochOptFormat

max_issues_repo_name: joaquimg/StochOptFormat

max_forks_repo_name: joaquimg/StochOptFormat


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

JLLWrappers.@generate_main_file_header("SBML")
JLLWrappers.@generate_main_file("SBML", UUID("bb12108a-f4ef-5f88-8ef3-0b33ff7017f1"))
end  # module SBML_jll


text:
using SurveyDataWeighting
using Test
using DelimitedFiles

#
# This is a simple test using real data, trying
# to match up a sample of FRS data from 2018 to some Scottish Targets
#
# This is some target data from Scotland 2020
# There is discussion of these targets in thhe ScottishTaxBenefit Model docs.
#
TARGETS = [
    1_340_609.0, # 1 - M- Total in employment- aged 16+ - empl, unempl by sex
    60_635, # 2 - M- Total unemployed- aged 16+
    1_301_248, # 3 - F- Total in employment- aged 16+
    59_302, # 4 - F- Total unemployed- aged 16+
    370_502, # 5 - private rented+rent free -- hhlds by tenure type, ommitting O-Os
    282_482, # 6 - housing association
    314_433, # 7 - las etc rented
    139_982, # 8 - M- 0 - 4 -- population in 5-year age / sex bands
    153_297, # 9 - M- 5 - 9
    150_487, # 10 - M- 10 – 14
    144_172, # 11 - M- 15 - 19
    176_066, # 12 - M- 20 - 24
    191_145, # 13 - M- 25 - 29
    182_635, # 14 -  M- 30 - 34
    172_624, # 15 -  M- 35 - 39
    156_790, # 16 - M- 40 - 44
    174_812, # 17 -  M- 45 - 49
    193_940, # 18 -  M- 50 - 54
    190_775, # 19 -  M- 55 - 59
    166_852, # 20 -  M- 60 - 64
    144_460, # 21 -  M- 65 - 69
    132_339, # 22 -  M- 70 - 74
    87_886, # 23 -  M- 75 - 79
    104_741, # 24 - M- 80’+
    131_733, # 25 - F- 0 - 4
    146_019, # 26 - F- 5 - 9
    144_187, # 27 - F- 10 - 14
    137_786, # 28 - F- 15 - 19
    171_390, # 29 - F- 20 - 24
    191_110, # 30 - F- 25 - 29
    186_828, # 31 -  F- 30 - 34
    179_898, # 32 -  F- 35 - 39
    162_642, # 33 - F- 40 - 44
    186_646, # 34 -  F- 45 - 49
    207_150, # 35 -  F- 50 - 54
    202_348, # 36 -  F- 55 - 59
    177_841, # 37 -  F- 60 - 64
    154_984, # 38 -  F- 65 - 69
    146_517, # 39 -  F- 70 - 74
    108_065, # 40 -  F- 75 - 79
    165_153, # 41 -  F- 80+
    439_000, # 42 - 1 adult: male -- household compositions
    467_000, # 43 - 1 adult: female
    797_000, # 44 - 2 adults
    70_000, # 45 - 1 adult 1 child
    66_000, # 46 - 1 adult 2+ children
    448_000, # 47 - 2+ adults 1+ children
    190_000, # 48 - 3+ adults
    77_842, # 49 - CARER’S ALLOWANCE - in receipt
    127_307, # 50 - AA - in receipt
    431_461 ] # 51 PIP or DLA

#
# subsample of scottish FRS 2018
#
data = readdlm( "../data/scotmat.csv")

@testset "Weighting Tests using 2018 Scottish FRS Subset" begin

    nr,nc = size(data)

    @test nc == size(TARGETS)[1]

    hhlds_in_popn = sum( TARGETS[42:48]) # total num hhlds

    initial_weights = ones(nr)*hhlds_in_popn/nr
    @test sum(initial_weights) ≈ hhlds_in_popn
    initial_weighted_popn = (initial_weights' * data)'
    println( "initial-weighted_popn vs targets" )
    for c in 1:nc
        diffpc = 100*(initial_weighted_popn[c]-TARGETS[c])/TARGETS[c]
        println( "$c $(TARGETS[c]) $(initial_weighted_popn[c]) $diffpc%")
    end

    wchi = do_chi_square_reweighting( data, initial_weights, TARGETS )
    println( "direct chi-square results $(wchi)")

    weighted_popn_chi = (wchi' * data)'
    # println( "wchisq; got $weighted_popn_chi")
    @test weighted_popn_chi ≈ TARGETS

    lower_multiple = 0.25 # any smaller min and d_and_s_constrained fails on this dataset
    upper_multiple = 4.8
    for m in [constrained_chi_square] #instances( DistanceFunctionType ) all other methods fail!
      println( "on method $m")
      rw = do_reweighting(
            data               = data,
            initial_weights    = initial_weights,
            target_populations = TARGETS,
            functiontype       = m,
            lower_multiple     = lower_multiple,
            upper_multiple     = upper_multiple,
            tolx               = 0.000001,
            tolf               = 0.000001 )
      println( "results for method $m = $(rw.rc)" )
      weights = rw.weights
      weighted_popn = (weights' * data)'
      println( "weighted_popn = $weighted_popn" )
      @test weighted_popn ≈ TARGETS
      if m != chi_square
         for w in weights # check: what's the 1-liner for this?
            @test w > 0.0
         end
      else
         @test weights ≈ wchi # chisq the direct way should match chisq the iterative way
      end
      if m in [constrained_chi_square, d_and_s_constrained ]
         # check the constrainted methods keep things inside ll and ul
         for r in 1:nrows
            @test weights[r] <= initial_weights[r]*upper_multiple
            @test weights[r] >= initial_weights[r]*lower_multiple
         end
      end
    end # meth loop

end # testset

max_stars_repo_name: grahamstark/SurveyDataWeighting

max_issues_repo_name: grahamstark/SurveyDataWeighting

max_forks_repo_name: grahamstark/SurveyDataWeighting


text:
#!/usr/bin/env julia
#
# Copyright 2022 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
module assignment9

using Plots

function z(Ppr, Tpr)

    t = 1 / Tpr

    a₁  = 0.317842
    a₂  = 0.382216
    a₃  = -7.768354
    a₄  = 14.290531
    a₅  = 0.000002
    a₆  = -0.004693
    a₇  = 0.096254
    a₈  = 0.166720
    a₉  = 0.966910
    a₁₀ = 0.063069
    a₁₁ = -1.966847
    a₁₂ = 21.0581
    a₁₃ = -27.0246
    a₁₄ = 16.23
    a₁₅ = 207.783
    a₁₆ = -488.161
    a₁₇ = 176.29
    a₁₈ = 1.88453
    a₁₉ = 3.05921

    A = a₁ * t * exp(a₂ * (1 - t) ^ 2) * Ppr
    B = a₃ * t + a₄ * t ^ 2 + a₅ * t ^ 6 * Ppr ^ 6
    C = a₉ + a₈ * t * Ppr + a₇ * t ^ 2 * Ppr ^ 2 + a₆ * t ^ 3 * Ppr ^ 3
    D = a₁₀ * t * exp(a₁₁ * (1 - t) ^ 2)
    E = a₁₂ * t + a₁₃ * t ^ 2 + a₁₄ * t ^ 3
    F = a₁₅ * t + a₁₆ * t ^ 2 + a₁₇ * t ^ 3
    G = a₁₈ + a₁₉ * t

    y = D * Ppr / ((1 + A ^ 2) / C - A ^ 2 * B / C ^ 3)

    D * Ppr * (1 + y + y ^ 2 - y ^ 3) / (D * Ppr + E * y ^ 2 - F * y ^ G) / (1 - y) ^ 3
end

function z_plot()

    p = #Add plotting code

    return p #Make sure to return the plot
end

function save_png()
    ENV["GKSwstype"]="nul"
    p = z_plot();
    savefig(p, "z.png");
end

function save_tikz()
    p = z_plot();
    savefig(p, "z.tikz");
end

export save_tikz, save_png, z_plot

end


text:
# row_major_reshape(X::AbstractArray, size...) = permutedims(reshape(X, reverse([size...])...), length(size):-1:1)

"""
    row_major_reshape(Q::AbstractArray, shapes)

Reshaping arrays in a different order.
"""
function row_major_reshape(Q::AbstractArray, shapes)
    shapes = reverse(shapes)
    perm = collect(length(shapes):-1:1)
    return permutedims(reshape(Q, shapes), perm)
end

function chop_op(O::AbstractArray, tol = 1e-8)
    tmp_r = (abs.(real.(O)) .> tol) .* real.(O)
    tmp_i = (abs.(imag.(O)) .> tol) .* imag.(O)
    return tmp_r + 1im .* tmp_i
end

"""
    gaussian(x, mu, sig)

Gaussian function.
"""
function gaussian(x, mu, sig)
    return exp.(- 0.5 * (x .- mu).^2 / sig^2)
end

function gaussian_derivative(x, mu, sig)
    return - (x .- mu) ./ sig^2 .* exp.(- 0.5 * (x .- mu).^2 / sig^2)
end

function trunc_op(op::AbstractArray, states)
    N_trunc = size(states)[1]
    # qstates = [qtp.Qobj(states[i], dims = [[N_s, N_s], [1, 1]]) for i in range(N_trunc)]
    res = spzeros(N_trunc, N_trunc)
    for i in range(1, N_trunc, step = 1)
        for j in range(1, N_trunc, step = 1)
            res += (adjoint(states[i, :]) * op * states[j, :]) * projection(N_trunc, i, j, true)
        end
    end
    return chop_op(res)
end

function eigensystem(A::AbstractArray; k = 6, v0 = nothing, sigma = nothing, krylovdim = 30)
    is_A_hermitian = ishermitian(A)

    if v0 === nothing
        v0 = rand(eltype(A), size(A, 1))
    end
    v0 /= norm(v0)

    if sigma === nothing
        vals, vecs, info = eigsolve(A, v0, k, ishermitian = is_A_hermitian)
    else
        # fac = factorize(A - sigma * I)
        # vals, vecs, info = eigsolve(x -> fac \ x, v0, k, ishermitian = is_A_hermitian)
        # vals = (1 .+ sigma * vals) ./ vals

        A_s = A - sigma * I

        P_cpu = ilu(A_s, τ = 0.01)
        vals, vecs, info = eigsolve(x -> cg(A_s, x; Pl = P_cpu, maxiter = 500), v0, k, ishermitian = true, krylovdim = krylovdim)
        vals = (1 .+ sigma * vals) ./ vals
    end
    if is_A_hermitian
        vals = real.(vals)
    end
    vals = vals[1:k]
    vecs = hcat(vecs...)[:, 1:k]
    return vals, vecs
end

function ptrace(Q::AbstractArray, sel, rd)
    nd = length(rd)
    dkeep = rd[sel]
    qtrace = filter(e->e∉sel,1:nd)
    dtrace = rd[qtrace]

    vmat = row_major_reshape(Q, (prod(rd), 1))
    vmat = row_major_reshape(vmat, rd)
    topermute = []
    append!(topermute, sel)
    append!(topermute, qtrace)
    vmat = permutedims(vmat, topermute)
    vmat = permutedims(vmat, nd:-1:1)
    vmat = row_major_reshape(vmat, (prod(dkeep), prod(dtrace)))
    return vmat * adjoint(vmat)
end

function entropy_vn(rho::AbstractArray, base = 0, tol = 1e-15)
    vals = eigvals(rho)
    indexes = abs.(vals) .> tol
    if 1 ∈ indexes
        nzvals = vals[indexes]
    #     nzvals = vals[vals .!= 0]
        if base != 0
            logvals = log.(base, Complex.(nzvals))
        else
            logvals = log.(Complex.(nzvals))
        end
        return -real(sum(nzvals .* logvals))
    else
        return 0
    end
end

function entanglement(psi, subspaces, size)
    norm = sqrt(sum(abs.(psi).^2))
    ψ = deepcopy(psi)
    if norm > 1e-5
        ψ = ψ ./ norm
    else
        return 0
    end
    rho_rabi = ptrace(ψ, subspaces, size)
    entropy = entropy_vn(rho_rabi)
    return (entropy > 0) * entropy
end


text:
using Distributed
@everywhere using Mamba

## Data
pumps = Dict{Symbol, Any}(
  :y => [5, 1, 5, 14, 3, 19, 1, 1, 4, 22],
  :t => [94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5]
)
pumps[:N] = length(pumps[:y])

## Model Specification
model = Model(

  y = Stochastic(1,
    (theta, t, N) ->
      UnivariateDistribution[(
        lambda = theta[i] * t[i];
        Poisson(lambda)) for i in 1:N
      ],
    false
  ),

  theta = Stochastic(1,
    (alpha, beta) -> Gamma(alpha, 1 / beta),
    true
  ),

  alpha = Stochastic(
    () -> Exponential(1.0)
  ),

  beta = Stochastic(
    () -> Gamma(0.1, 1.0)
  )

)

## Initial Values
inits = [
  Dict(:y => pumps[:y], :alpha => 1.0, :beta => 1.0,
       :theta => rand(Gamma(1.0, 1.0), pumps[:N])),
  Dict(:y => pumps[:y], :alpha => 10.0, :beta => 10.0,
       :theta => rand(Gamma(10.0, 10.0), pumps[:N]))
]

## Sampling Scheme
scheme = [Slice([:alpha, :beta], 1.0, Univariate),
          Slice(:theta, 1.0, Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, pumps, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

## Posterior Predictive Distribution
ppd = predict(sim, :y)
describe(ppd)


text:
#module MapPlot

import .SimpleGraphics
using SparseArrays
#using NMR.Scaling
#using NMR.SimplePlot

export Spy

function Spy(A::AbstractSparseArray{T,Int64};
    g=SimpleGraphics.Group(),
    frame=true,
    axes=false,
    show=true,
    Reverse=[false,true],
    FrameW=1000,
    FrameH=1000) where {T}

  (n,m)=size(A);
  xrange = 1:m;
  yrange = 1:n;

  I,J,V=findnz(A);
  xscale=NiceScale(0.7,m+0.3);
  yscale=NiceScale(0.7,n+0.3);

  revrs=Reverse
  if frame g=Frame(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end
  if axes  g=Axes(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end

  (xreverse,yreverse)=Reverse

  x0=scaled(minimum(xrange),xscale)
  x1=scaled(maximum(xrange),xscale)
  y0=scaled(minimum(yrange),yscale)
  y1=scaled(maximum(yrange),yscale)

  Q=[FrameW*(x1-x0)/(n-1) 0 ; 0 FrameH*(y1-y0)/(m-1)]

  style=SimpleGraphics.GraphicsAttributes(["fill"=>"hsl(,50%,60%)"])
  maxv=maximum(abs.(V));
  for k=1:length(V)
    (i,j,v) = (I[k],J[k],V[k]) ;
    style=SimpleGraphics.GraphicsAttributes(["fill"=>"hsl($(180/π*angle(v)),$(round(100*abs(v)/maxv))%,60%)"])
    if(v != 0)
      circ=SimpleGraphics.Circle(m-j,i-1,0.35);
      circ.attr=style;
      SimpleGraphics.push!(g,SimpleGraphics.translate(SimpleGraphics.transform(circ,Q),(x0)*FrameW,(y0)*FrameH))
    end
  end
  if show Show(g) end
  return(g)
end

function Pcolor(A::Array{T,2};
    g=SimpleGraphics.Group(),
    frame=true,
    axes=false,
    show=true,
    Reverse=[false,true],
    FrameW=1000,
    FrameH=1000) where {T}

  (n,m)=size(A);
  xrange = 1:m;
  yrange = 1:n;

  xscale=NiceScale(0.7,m+0.3);
  yscale=NiceScale(0.7,n+0.3);

  revrs=Reverse
  if frame g=Frame(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end
  if axes  g=Axes(xscale,yscale,FrameW=FrameW,FrameH=FrameH,Reverse=revrs) end

  (xreverse,yreverse)=Reverse

  x0=scaled(minimum(xrange),xscale)
  x1=scaled(maximum(xrange),xscale)
  y0=scaled(minimum(yrange),yscale)
  y1=scaled(maximum(yrange),yscale)

  Q=[FrameW*(x1-x0)/(n-1) 0 ; 0 FrameH*(y1-y0)/(m-1)]

  style=SimpleGraphics.GraphicsAttributes(["fill"=>"hsl(,50%,60%)"])
  maxv=maximum(abs.(A));
  for i=1:m,j=1:n
    v=A[i,j] ;
    style=SimpleGraphics.GraphicsAttributes(["fill"=>"hsl($(180/π*angle(v)),60%,$(round(100*abs(v)/maxv))%)"])
    if(v != 0)
      circ=SimpleGraphics.Circle(m-j,i-1,0.75);
      circ.attr=style;
      SimpleGraphics.push!(g,SimpleGraphics.translate(SimpleGraphics.transform(circ,Q),(x0)*FrameW,(y0)*FrameH))
    end
  end
  if show Show(g) end
  return(g)
end

#end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
using BinaryBuilder

include("../../../fancy_toys.jl")

name = "CUDA_full"
version = v"9.2.148"

sources_linux = [
    FileSource("https://developer.nvidia.com/compute/cuda/9.2/Prod2/local_installers/cuda_9.2.148_396.37_linux",
               "f5454ec2cfdf6e02979ed2b1ebc18480d5dded2ef2279e9ce68a505056da8611", "installer.run")
]
sources_macos = [
    FileSource("https://developer.nvidia.com/compute/cuda/9.2/Prod2/local_installers/cuda_9.2.148_mac",
               "defb095aa002301f01b2f41312c9b1630328847800baa1772fe2bbb811d5fa9f", "installer.dmg")
]
sources_windows = [
    FileSource("https://developer.nvidia.com/compute/cuda/9.2/Prod2/local_installers2/cuda_9.2.148_win10",
               "7d99a6d135587d029c2cf159ade4e71c02fc1a922a5ffd06238b2bde8bedc362", "installer.exe")
]

script = raw"""
cd ${WORKSPACE}/srcdir

# use a temporary directory to avoid running out of tmpfs in srcdir on Travis
temp=${WORKSPACE}/tmpdir
mkdir ${temp}

apk del ninja   # hack
apk add p7zip

mkdir ${prefix}/cuda
if [[ ${target} == x86_64-linux-gnu ]]; then
    sh installer.run --tmpdir="${temp}" --extract="${temp}"
    cd ${temp}
    sh cuda-linux.*.run --noexec --keep
    cd pkg
    find .

    # clean-up
    rm -r libnsight libnvvp nsightee_plugins jre doc

    mv * ${prefix}/cuda
elif [[ ${target} == x86_64-apple-darwin* ]]; then
    7z x installer.dmg -o${temp}
    cd ${temp}
    tar -xzf CUDAMacOSXInstaller/CUDAMacOSXInstaller.app/Contents/Resources/payload/cuda_mac_installer_tk.tar.gz
    cd Developer/NVIDIA/CUDA-*/
    find .

    # clean-up
    rm -r libnsight libnvvp nsightee_plugins jre doc

    mv * ${prefix}/cuda
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    7z x installer.exe -o${temp}
    cd ${temp}
    find .

    mv EULA.txt ${prefix}/cuda

    for project in cuobjdump memcheck nvcc cupti nvdisasm curand cusparse npp cufft \
                   cublas cudart cusolver nvrtc nvgraph nvprof nvprune; do
        [[ -d ${project} ]] || { echo "${project} does not exist!"; exit 1; }
        cp -a ${project}/* ${prefix}/cuda
    done

    # NVIDIA Tools Extension Library
    7z x "CUDAVisualStudioIntegration/NVIDIA NVTX Installer.x86_64".*.msi -o${temp}/nvtx_installer
    find nvtx_installer
    for file in nvtx_installer/*.*_*; do
        mv $file $(echo $file | sed 's/\.\(\w*\)_.*/.\1/')
    done
    mv nvtx_installer/*.dll ${prefix}/cuda/bin
    mv nvtx_installer/*64_*.lib ${prefix}/cuda/lib/x64
    mv nvtx_installer/*32_*.lib ${prefix}/cuda/lib/Win32
    mv nvtx_installer/*.h ${prefix}/cuda/include
    find .

    # fixup
    chmod +x ${prefix}/cuda/bin/*.{exe,dll}

    # clean-up
    rm ${prefix}/cuda/*.nvi
fi

cd ${prefix}/cuda
install_license EULA.txt
"""

products = Product[
    # this JLL isn't meant for use by Julia packages, but only as build dependency
]

dependencies = []

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux, script,
                   [Platform("x86_64", "linux")], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("x86_64-apple-darwin14")
    build_tarballs(non_reg_ARGS, name, version, sources_macos, script,
                   [Platform("x86_64", "macos")], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("x86_64-w64-mingw32")
    build_tarballs(ARGS, name, version, sources_windows, script,
                   [Platform("x86_64", "windows")], products, dependencies;
                   skip_audit=true)
end

max_stars_repo_name: koehlerson/Yggdrasil

max_issues_repo_name: koehlerson/Yggdrasil

max_forks_repo_name: koehlerson/Yggdrasil


text:
using Impero
using Test
using SymbolicUtils
import SymbolicUtils: Chain, Postwalk, Fixpoint

@wrapper a=1 b=2

@testset "symbolic utils swapping" begin
    println("for $a = ", a.data, ", $b = ", b.data)
    println("changing a+b to a*b")
    c = a+b
    check_c = a*b
    r1 = @acrule ~ra + ~rb => ~ra * ~rb
    new_c = Fixpoint(Postwalk(Chain([r1])))(c)
    @test compute(check_c) == compute(new_c)
    println("changing sin(a+b) * sin(a+a+b) to sin(a*b) * sin(a*a*b)")
    c = sin(a+b) * sin(a+a+b)
    check_c = sin(a*b) * sin(a*a*b)
    new_c = Fixpoint(Postwalk(Chain([r1])))(c)
    @test compute(check_c) == compute(new_c)
end


text:
include("./Operation.jl")
using .Operation
using QML
qmlfunction("float_to_frac", Operation.float_to_frac)
qmlfunction("transH", Operation.transH)
qmlfunction("mul", Operation.mul)
qmlfunction("minus", Operation.minus)
qmlfunction("plus", Operation.plus)
qmlfunction("rang", Operation.rang)
qmlfunction("obr", Operation.obr)
qmlfunction("transfake", Operation.transfake)
qmlfunction("det", Operation.det)
qmlfunction("tr", Operation.tr)
loadqml("dev//SolMatrix//src//QMLinterface.qml")
exec_async()

max_stars_repo_name: NikitaKutelev/kurs

max_issues_repo_name: NikitaKutelev/kurs

max_forks_repo_name: NikitaKutelev/kurs


text:
@testset "1.2.4 Improper" begin
    include("1.2.4.2 (d x)^m (a x^q+b x^n+c x^(2 n-q))^p.jl")
end

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.4 Improper/.jl


text:
# FizzBuzzの値を返す関数を定義 --- (*1)
function fizzbuzz(n)
  if n % 3 == 0 && n % 5 == 0
    return "FizzBuzz"
  elseif n % 3 == 0
    return "Fizz"
  elseif n % 5 == 0
    return "Buzz"
  else
    return n
  end
end

# 繰り返し関数を実行 --- (*2)
for i = 1:100
  println(fizzbuzz(i))
end

max_stars_repo_name: KiYOKing/fizzbuzz

max_issues_repo_name: KiYOKing/fizzbuzz

max_forks_repo_name: KiYOKing/fizzbuzz


text:
# @references:
# 1. `https://github.com/evetion/GeoArrays.jl/blob/master/src/geoarray.jl`

# using Stars

abstract type AbstractGeoArray{T,N} <: AbstractArray{T,N} end

# abstract type AbstractGeoArray{T,N,D,A} <: AbstractDimensionalArray{T,N,D,A} end

const RealOrMissing = Union{Missing,Real}
const TypeCRS = Union{WellKnownText,AbstractString}
"""
    GeoArray(A::AbstractArray{<:RealOrMissing}, f::AffineMap = default_affinemap, crs::WellKnownText = crs2wkt(""))
    GeoArray(A::AbstractArray{<:RealOrMissing}, f::AffineMap, crs::AbstractString)
    GeoArray(A::AbstractArray{<:RealOrMissing}, x::AbstractRange, y::AbstractRange, args...)

    GeoArray(fn::AbstractString, bands = nothing)
    GeoArray(ga::AbstractGeoArray; vals::AbstractArray{<:RealOrMissing})

Construct a GeoArray from any Array. Meanwhile, `AffineMap` and `CRS` will be generated.

# Arguments

- `A`: An AbstractArray, at least 2d Array

# Examples
```julia
ga = GeoArray(rand(10,10,1))
GeoArray(rand(10,10))

GeoArray(ga, vals = ga.A)                        # copy `f` and `crs` from `ga`
GeoArray(rand(10,10), 1:10, 2:11)                # define coordinate by `x` and `y`
GeoArray(rand(7, 4), bbox([70, 15, 140, 55]...)) # define coordinate by `bbox`

# 4-d array
ga = GeoArray(rand(10, 10, 2, 3))
GeoArray(ga, vals = ga.A)
```
"""
Base.@kwdef mutable struct GeoArray{T<:RealOrMissing,N} <: AbstractGeoArray{T,N}
    A::AbstractArray{T,N}
    f::AffineMap
    crs::TypeCRS
    names
    time
end

# using Stars
function GeoArray(A::AbstractArray{<:RealOrMissing,N},
    f::AffineMap=default_affinemap,
    crs::TypeCRS=crs2wkt(""),
    names=nothing) where {N}

    # @show N
    if N == 1
        error("`A` should be matrix or array!")
    elseif N == 2
        A = reshape(A, (size(A)..., 1))
    end
    ntime = size(A) |> last
    if names !== nothing
        names = length(names) > ntime ? names[1:ntime] : nothing
    end
    GeoArray(; A=A, f=f, crs=crs2wkt(crs), names=names, time=nothing)
end

# # also suit for high-dimension
function GeoArray(A::AbstractArray{<:RealOrMissing},
    x::AbstractRange, y::AbstractRange, args...)

    if size(A)[1:2] != (length(x), length(y))
        # Note that this function takes *center coordinates*.
        error("Size of `GeoArray` $(size(A)) does not match size of (x,y): $((length(x), length(y))).")
    end
    f = unitrange_to_affine(x, y)
    # b = st_bbox(x, y)
    GeoArray(A, f, args...)
end

function GeoArray(A::AbstractArray{<:RealOrMissing}, b::bbox; proj=4326)
    ga = GeoArray(A)
    st_bbox!(ga, b)
    st_crs!(ga, proj)
    ga
end

GeoArray(fn::AbstractString, bands=nothing) = st_read(fn, bands)

function GeoArray(ga::AbstractGeoArray; vals::AbstractArray{<:RealOrMissing}, names=nothing)
    if names === nothing
        names = ga.names
    end
    GeoArray(vals, ga.f, ga.crs, names)
end

rast = GeoArray
export GeoArray, rast, AbstractGeoArray

include("st_[.jl")
include("st_Ops.jl")


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

"""
    charge_trapflt!(samples::AbstractVector{<:RealOrSIMD{<:AbstractFloat}}, navg::Integer, ngap::Integer)

Apply a trapezoidal FIR filter to a charge signal in `samples`.
"""
function charge_trapflt!(samples::AbstractVector{<:RealOrSIMD{<:AbstractFloat}}, navg::Integer, ngap::Integer)
    @fastmath begin
        T = eltype(samples)
        idxs = eachindex(samples)
        n = length(idxs)

        norm_factor = inv(T(navg))

        navg >= 1 && ngap >= 0 || throw(ArgumentError("Require navg >= 1 and ngap >= 0"))
        (2 * navg + ngap) <= length(idxs) || throw(ArgumentError("filter must not be longer than input"))

        flt_shift = 0

        offs1 = flt_shift
        offs2 = flt_shift + navg
        offs3 = flt_shift + navg + ngap
        offs4 = flt_shift + navg + ngap + navg

        padding_value = samples[last(idxs)]

        acc::T = zero(T)

        @inbounds @simd for j in idxs[offs3+1:offs4] acc += samples[j] end
        @inbounds @simd for j in idxs[offs1+1:offs2] acc -= samples[j] end
        samples[idxs[1]] = acc * norm_factor

        @inbounds @simd for i in idxs[2:n-offs4]
            acc = acc + samples[i + offs1] - samples[i + offs2] - samples[i + offs3] + samples[i + offs4]
            samples[i] = acc * norm_factor
        end

        @inbounds @simd for i in idxs[n-offs4+1:n-offs3]
            acc = acc + samples[i + offs1] - samples[i + offs2] - samples[i + offs3] + padding_value
            samples[i] = acc * norm_factor
        end

        @inbounds @simd for i in idxs[n-offs3+1:n-offs2]
            acc = acc + samples[i + offs1] - samples[i + offs2]
            samples[i] = acc * norm_factor
        end

        @inbounds @simd for i in idxs[n-offs2+1:n-offs1]
            acc = acc + samples[i + offs1] - padding_value
            samples[i] = acc * norm_factor
        end
    end

    samples
end
export charge_trapflt!


text:
# Standard rectlin layer: f(A,x,b)=max(A*x+b,gamma(Ax+b)) where gamma=0.25
# Note that the bias (a column vector) gets expanded here to match the dimension of A*x

gamma=0.25

FkinklinAXplusBias(A,X,b)=begin; a=A*X+b*ones(1,size(X,2)); return (max(a,gamma*a),gamma+(1-gamma)*(a.>0)); end

function FkinklinAXplusBias_inplace(value,aux,A,X,b)
    a=A*X+b*ones(1,size(X,2))
    copy!(value,max(a,gamma*a))
    copy!(aux,gamma+(1-gamma)*(a.>0))
end

function DkinklinAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X,b)
    if derivativeIDX==1
        axpy!(1.0,(grad_c.*faux_c)*X',grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,A'*(grad_c.*faux_c),grad_n)
    elseif derivativeIDX==3
        axpy!(1.0,sum(grad_c.*faux_c,2),grad_n)
    end
end

if PROC=="GPU"

    function FkinklinAXplusBias_inplace(value,aux,A::CudaArray,X::CudaArray,b::CudaArray)
        FAXplusBias_inplace(value,[],A,X,b)
        copy!(aux,value)
        kinklin!(value,value)
    end

    function DkinklinAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray,b::CudaArray)
        tmp=CudaArray(Float64,size(grad_c)); fill!(tmp,0.0)
        A_emult_Bg0!(grad_c,faux_c,tmp)
        scale!((1-gamma),tmp)
        axpy!(gamma,grad_c,tmp)

        if derivativeIDX==1
            gemm!('N','T',1.0,tmp,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,tmp,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float64,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,tmp,ons,1.0,grad_n)
            free(ons)
        end
        free(tmp)
    end

    function DkinklinAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float32},X::CudaArray{Float32},b::CudaArray{Float32})
        tmp=CudaArray(Float32,size(grad_c)); fill!(tmp,0.0)
        A_emult_Bg0!(grad_c,faux_c,tmp)
        scale!((1-gamma),tmp)
        axpy!(gamma,grad_c,tmp)

        if derivativeIDX==1
            gemm!('N','T',Float32(1.0),tmp,X,Float32(1.0),grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',Float32(1.0),A,tmp,Float32(1.0),grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float32,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',Float32(1.0),tmp,ons,Float32(1.0),grad_n)
            free(ons)
        end
        free(tmp)
    end

end

Derivative[FkinklinAXplusBias]=DkinklinAXplusBias
Inplace[FkinklinAXplusBias]=FkinklinAXplusBias_inplace

kinklinAXplusBias(A,X,b)=ADnode(FkinklinAXplusBias,[A X b])
export kinklinAXplusBias

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# Functions for dynamic systems
# Increment of x[ix] through RK4-integration
# dxdh is the array of derivatives, dt is the time increment
function dx_rk4(dxdt, x, ix, dt)
    x = collect(x)
    incr = zeros(length(x))
    k1 = dxdt[ix](x...)
    incr[ix] = dt*k1/2
    k2 = dxdt[ix]((x+incr)...)
    incr[ix] = dt*k2/2
    k3 = dxdt[ix]((x+incr)...)
    incr[ix] = dt*k3
    k4 = dxdt[ix]((x+incr)...)
    dt/6*(k1 + 2k2 + 2k3 + k4)
end
# Integrate dynamical system through n intervals of length dt
function dynamical_system(x0, dxdt, dt, n)
    m = length(x0)
    x = zeros(n,m)
    x[1,:] = collect(x0)
    for t=1:n-1
        dx = zeros(m)
        for ix=1:m
            dx[ix] = dx_rk4(dxdt, x[t,:], ix, dt)
        end
        x[t+1,:] = x[t,:][:] + dx
    end
    x
end
# Rössler system
rossler_eq(a, b, c) = (
    (x,y,z) -> -(y+z),
    (x,y,z) -> x + a*y,
    (x,y,z) -> b + z*(x-c))

a=.25; b=.25; c=4
x0 = zeros(3)
m = zeros(3000,24)
for i=0:11
    x = randn(3)
    rossler_data = dynamical_system(x0, rossler_eq(a,b,c), .05, 4000)
    m[:,2i.+(1:2)] .= rossler_data[1001:end,1:2]
end
writedlm("rossler.txt",m,'\t')


text:
abstract type HeterosynapticConstraint end
abstract type HeterosynapticTarget end
abstract type HeterosynapticMethod end

struct HetUpperLimit <: HeterosynapticConstraint
  wsum_max::Float64
  wmin::Float64
  wmax::Float64
  tolerance::Float64
end
struct HetStrictSum <: HeterosynapticConstraint
  wsum_max::Float64
  wmin::Float64
  wmax::Float64
  tolerance::Float64
end

@inline function hardbounds(x::Float64,hc::HeterosynapticConstraint)
  hardbounds(x,hc.wmin,hc.wmax)
end

function _hetconstraint_vs_sum(sums::Vector{Float64},c::HetUpperLimit)
  _lim  =  c.wsum_max + c.tolerance
  return sums .> _lim
end
function _hetconstraint_vs_sum(sums::Vector{Float64},c::HetStrictSum)
  return .!isapprox.(sums,c.wsum_max;atol=c.tolerance)
end

struct HetIncoming <: HeterosynapticTarget end
struct HetOutgoing <: HeterosynapticTarget end
struct HetBoth <: HeterosynapticTarget end

struct HetAdditive <: HeterosynapticMethod end
struct HetMultiplicative <: HeterosynapticMethod end

struct PlasticityHeterosynapticSpikeTriggered{
    HC<:HeterosynapticConstraint,
    HM<:HeterosynapticMethod,
    HT<:HeterosynapticTarget} <:PlasticityRule
  constraint::HC
  method::HM
  target::HT
  Δt_min_update::Float64
  t_last_update::Vector{Float64}
end
function PlasticityHeterosynapticSpikeTriggered(
    nneus::Integer, Δt::Float64,
    hc::HeterosynapticConstraint,
    hm::HeterosynapticMethod,
    ht::HeterosynapticTarget)
  t_last = fill(-Inf,nneus)
  return PlasticityHeterosynapticSpikeTriggered(hc,hm,ht,Δt,t_last)
end

function reset!(plast::PlasticityHeterosynapticSpikeTriggered)
  fill!(plast.t_last_update,-Inf)
  return nothing
end

function plasticity_update!(t::R,::R,
     pspost::PopulationState,conn::Connection,::PopulationState,
     plast::PlasticityHeterosynapticSpikeTriggered) where R

  if !any(pspost.isfiring) # to optimize a little
    return nothing
  end
  idxs_fire = findall(pspost.isfiring)
  # check time
  to_keep = falses(length(idxs_fire))
  for (i,neu) in enumerate(idxs_fire)
    if (t-plast.t_last_update[neu]) > plast.Δt_min_update
      # keep it, reset timer
      to_keep[i] = true
      plast.t_last_update[neu] = t
    end
  end
  if !any(to_keep)
    return nothing
  end
  # meh, it's not much faster if regularize less neurons
  plast.t_last_update[idxs_fire].= t
  # check if must be normalized, returns sum and n elements
  to_change,sumvals,nvals = _hetplast_check(idxs_fire,conn.weights,
    plast.target,plast.constraint)
  # correct the weights only if needed
  if !any(to_change)
    return nothing
  end
  # apply correction
  _apply_hetplast_spiketriggered!(to_change,idxs_fire,sumvals,nvals,
      conn.weights, plast.target,plast.method,plast.constraint)
  return nothing
end

# incoming : check row value
# faster to go over all elements, than consider few ones :-(
function _hetplast_check(idxs::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetIncoming,constraint::HeterosynapticConstraint) where {I<:Integer}
  weights_nonzeros = nonzeros(weights)
  weights_rows = rowvals(weights)
  N = size(weights,1)
  n_el = fill(0,N)
  sums_all = fill(0.0,N)
  for (_r,row) in enumerate(weights_rows)
    sums_all[row] += weights_nonzeros[_r]
    n_el[row] += 1
  end
  _sums_val = sums_all[idxs]
  _sums_nel = n_el[idxs]
  tochange = _hetconstraint_vs_sum(_sums_val,constraint)
  return tochange,_sums_val,_sums_nel
end

# outgoing : column value
function _hetplast_check(idxs::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetOutgoing,constraint::HeterosynapticConstraint) where {I<:Integer}
  weights_nonzeros = nonzeros(weights)
  nsp=length(idxs)
  _sums_val = fill(NaN,nsp)
  _sums_nel = fill(-1,nsp)
  # now, sums over columns
  for (k,col) in enumerate(idxs)
    _range = nzrange(weights,col)
    lr = length(_range)
    if lr > 0
      _sums_val[k] = sum(view(weights_nonzeros,_range))
      _sums_nel[k] = length(_range)
    end
  end
  tochange = _hetconstraint_vs_sum(_sums_val,constraint)
  return tochange,_sums_val,_sums_nel
end

# both incoming end outgoing. Returns two elements, one for
# outgoing one for incoming
function _hetplast_check(idxs::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetBoth,constraint::HeterosynapticConstraint) where {
      I<:Integer }
  weights_nonzeros = nonzeros(weights)
  weights_rows = rowvals(weights)
  N = size(weights,1)
  n_el = fill(0,N)
  sums_all = fill(0.0,N)
  for (_r,row) in enumerate(weights_rows)
    sums_all[row] += weights_nonzeros[_r]
    n_el[row] += 1
  end
  twonsp=2*length(idxs)
  _sums_val = Vector{Float64}(undef,twonsp)
  _sums_nel = Vector{Int64}(undef,twonsp)
  # sum over rows : odd elements
  _sums_val[1:2:end] = sums_all[idxs]
  _sums_nel[1:2:end] = n_el[idxs]
  # now, sums over columns : even elements
  for (k,col) in enumerate(idxs)
    _range = nzrange(weights,col)
    lr = length(_range)
    if lr > 0
      _sums_val[2k] = sum(view(weights_nonzeros,_range))
      _sums_nel[2k] = length(_range)
    end
  end
  tochange = _hetconstraint_vs_sum(_sums_val,constraint)
  return tochange,_sums_val,_sums_nel
end

# incoming additive (rows)
function _apply_hetplast_spiketriggered!(to_change::BitArray{1},
    neus::Vector{I},
    sumvals::Vector{F},nelvals::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetIncoming,
    ::HetAdditive,constraint::HeterosynapticConstraint) where {
      I<:Integer,F<:Real}
  weights_nonzeros = nonzeros(weights)
  weights_rows = rowvals(weights)
  N = size(weights,1)
  # feels like MATLAB !!!
  fix_vals_all = fill(0.0,N) # much empty!
  for (k,neu) in enumerate(neus)
    if to_change[k]
      fix_vals_all[neu] = (constraint.wsum_max - sumvals[k]) ./ nelvals[k]
    end
  end
  for (_r,row) in enumerate(weights_rows)
    _fix = fix_vals_all[row]
    if !iszero(_fix)
    weights_nonzeros[_r] = hardbounds(
          weights_nonzeros[_r]+_fix,constraint)
    end
  end
  return nothing
end

# outgoing additive (columns)
function _apply_hetplast_spiketriggered!(to_change::BitArray{1},
    neus::Vector{I},
    sumvals::Vector{F},nelvals::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetOutgoing,
    ::HetAdditive,constraint::HeterosynapticConstraint) where {
      I<:Integer,F<:Real}
  # pass over columns
  weights_nonzeros = nonzeros(weights)
  for (k,col) in enumerate(neus)
    if to_change[k]
      fixval_k = (constraint.wsum_max - sumvals[k]) / nelvals[k]
      for _r in nzrange(weights,col)
        weights_nonzeros[_r] = hardbounds(
          weights_nonzeros[_r]+fixval_k,constraint)
      end
    end
  end
  return nothing
end

# both outoing and incoming additive (rows and columns)
function _apply_hetplast_spiketriggered!(to_change::BitArray{1},
    _neus::Vector{I},
    sumvals::Vector{F},nelvals::Vector{I},
    weights::SparseMatrixCSC{Float64,Int64},
    ::HetBoth,
    ::HetAdditive,constraint::HeterosynapticConstraint) where {
      I<:Integer,F<:Real}
  # repeat for rows and columsn
  neus = repeat(_neus;inner=2)
  @assert length(neus) == length(to_change)
  # row part , odd elements
  weights_nonzeros = nonzeros(weights)
  weights_rows = rowvals(weights)
  N = size(weights,1)
  fix_vals_all = fill(0.0,N)
  for k in 1:2:length(neus)
    if to_change[k]
      row = neus[k]
      fix_vals_all[row] = (constraint.wsum_max - sumvals[k]) ./ nelvals[k]
    end
  end
  for (_r,row) in enumerate(weights_rows)
    _fix = fix_vals_all[row]
    if !iszero(_fix)
      weights_nonzeros[_r] = hardbounds(
          weights_nonzeros[_r]+_fix,constraint)
    end
  end
  # now, pass over columns : even elements
  for k in 2:2:length(neus)
    if to_change[k]
      col = neus[k]
      fix_val_k = (constraint.wsum_max - sumvals[k]) / nelvals[k]
      for _r in nzrange(weights,col)
        weights_nonzeros[_r] = hardbounds(
          weights_nonzeros[_r]+fix_val_k,constraint)
      end
    end
  end
  return nothing
end

# rewrite EASY plasticity for strict sum and both outgoing and incoming

struct PlasticityHeterosynapticApprox{
    HC<:HeterosynapticConstraint,
    HM<:HeterosynapticMethod,
    HT<:HeterosynapticTarget} <:PlasticityRule
  constraint::HC
  method::HM
  target::HT
  Δt_update::Float64
  _tcounter::Ref{Float64}
  Nel_rows::Vector{Int64} # pre-allocate memory
  Nel_cols::Vector{Int64}
  allocrows::Vector{Float64}
  alloccols::Vector{Float64}
  function PlasticityHeterosynapticApprox(
      Npost::Int64,Npre::Int64,
      Δt_update::Float64,
      hc::HC,hm::HM,ht::HT) where {
        HC <: HeterosynapticConstraint,
        HM <: HeterosynapticMethod,
        HT <: HeterosynapticTarget }
    _tcounter = Ref(0.0)
    Nel_rows = Vector{Int64}(undef,Npost)
    Nel_cols = Vector{Int64}(undef,Npre)
    allocrows = Vector{Float64}(undef,Npost)
    alloccols = Vector{Float64}(undef,Npre)
    return new{HC,HM,HT}(hc,hm,ht,
      Δt_update,_tcounter,Nel_rows,Nel_cols,alloccols,allocrows)
  end
end

function reset!(plast::PlasticityHeterosynapticApprox)
  plast._tcounter[] = zero(Float64)
  return nothing
end

function plasticity_update!(::R,dt::R,
     pspost::PopulationState,conn::Connection,pspre::PopulationState,
     plast::PlasticityHeterosynapticApprox) where R
  if plast._tcounter[] < plast.Δt_update
    plast._tcounter[] += dt
    return nothing
  end
  # reset timer
  plast._tcounter[] = zero(R)
  # find correction values at col/row , depending on target
  _het_plasticity_fix_rows!(plast.allocrows,plast.Nel_rows,
    conn.weights,plast.constraint,plast.method,plast.target)
  _het_plasticity_fix_cols!(plast.alloccols,plast.Nel_cols,
    conn.weights,plast.constraint,plast.method,plast.target)
  # apply the fix
  _het_plasticity_apply_fix!(
      plast.allocrows,plast.alloccols,
      conn.weights,
      plast.constraint,plast.method,plast.target)
  return nothing
end

function sum_and_count_over_rows!(rowsum::Vector{R},nels::Vector{I},
    M::SparseMatrixCSC{R,I}) where {R,I}
  Mnz = nonzeros(M) # direct access to weights
  fill!(rowsum,zero(R))
  fill!(nels,zero(I))
  @inbounds for (i,r) in enumerate(rowvals(M))
    rowsum[r] += Mnz[i]
    nels[r] += 1
  end
  return nothing
end
function sum_and_count_over_cols!(colsum::Vector{R},nels::Vector{I},
    M::SparseMatrixCSC{R,I}) where {R,I}
  Mnz = nonzeros(M) # direct access to weights
  ncols = size(M,2)
  @inbounds for col in 1:ncols
    rang = nzrange(M,col)
    if isempty(rang)
      nels[col] = 0
      colsum[col] = 0.0
    else
      nels[col] = length(rang)
      colsum[col] = sum(view(Mnz,rang))
    end
  end
  return nothing
end

function _het_plasticity_fix_rows!(alloc::Vector{Float64},nel::Vector{Int64},
    weights::SparseMatrixCSC,
    constraint::HeterosynapticConstraint,::HetAdditive,::Union{HetBoth,HetIncoming})
  sum_and_count_over_rows!(alloc,nel,weights)
  sum_max = constraint.wsum_max
  @. alloc = (sum_max - alloc )/ nel
  return nothing
end
function _het_plasticity_fix_rows!(alloc::Vector{Float64},nel::Vector{Int64},
    ::SparseMatrixCSC,::HeterosynapticConstraint,
    ::HeterosynapticMethod,::HetOutgoing)
  return nothing
end
function _het_plasticity_fix_cols!(alloc::Vector{Float64},nel::Vector{Int64},
    weights::SparseMatrixCSC,
    constraint::HeterosynapticConstraint,::HetAdditive,::Union{HetBoth,HetOutgoing})
  sum_and_count_over_cols!(alloc,nel,weights)
  sum_max = constraint.wsum_max
  @. alloc = (sum_max - alloc )/ nel
  return nothing
end
function _het_plasticity_fix_cols!(alloc::Vector{Float64},nel::Vector{Int64},
    weights::SparseMatrixCSC,
    constraint::HeterosynapticConstraint,::HeterosynapticMethod,::HetIncoming)
  return nothing
end

function _het_plasticity_apply_fix!(
    fixrows::Vector{Float64},fixcols::Vector{Float64},
    weights::SparseMatrixCSC,constraint::HetStrictSum,
    ::HetAdditive,::HetBoth)
  ncols = size(weights,2)
  rows = rowvals(weights)
  weights_nonzeros = nonzeros(weights)
  @inbounds for col in 1:ncols
    rang = nzrange(weights,col)
    fixcol=fixcols[col]
    for r in rang
      row = rows[r]
      fix_val = 0.5*(fixrows[row]+fixcol)
      weights_nonzeros[r] = hardbounds(
          weights_nonzeros[r]+fix_val,constraint)
    end
  end
end

#=
slower :-( , and bugged
function _het_plasticity_apply_fix!(
    fixrows::Vector{Float64},fixcols::Vector{Float64},
    weights::SparseMatrixCSC,constraint::HetStrictSum,
    ::HetAdditive,::HetBoth)
  rows = rowvals(weights)
  weights_nonzeros = nonzeros(weights)
  colptr = SparseArrays.getcolptr(weights)
  @inbounds for i in eachindex(weights_nonzeros)
    col = searchsortedfirst(colptr,i) # faster ?
    row = rows[i]
    fix_val = 0.5*(fixrows[row]+fixcols[col])
    weights_nonzeros[i] = hardbounds(
          weights_nonzeros[i]+fix_val,constraint)
  end
end
=#

# and here is the EASY plasticity, instead of spike triggered
# ... to rewrite to make it more uniform to the above

# struct PlasticityHeterosynapticEasy2{HetMeth<:HeterosynapticPlasticityMethod,
#     HetTarg<:HeterosynapticPlasticityTarget} <: PlasticityRule
#   Δt_update::Float64
#   method::HetMeth
#   target::HetTarg
#   _tcounter::Ref{Float64}
#   function PlasticityHeterosynapticEasy2(Δt_update::Float64,
#       method::M,target::T) where {M<:HeterosynapticPlasticityMethod,T<:HeterosynapticPlasticityTarget}
#     counter = Ref(0.0)
#     new{M,T}(Δt_update,method,target,counter)
#   end
# end

abstract type HeterosynapticPlasticityMethod end
struct HeterosynapticAdditive <: HeterosynapticPlasticityMethod
  wmin::Float64
  wmax::Float64
end
struct HeterosynapticMultiplicative <: HeterosynapticPlasticityMethod
  wmin::Float64
  wmax::Float64
end

abstract type HeterosynapticPlasticityTarget end
struct HeterosynapticIncoming <: HeterosynapticPlasticityTarget
  sum_max::Float64
end
struct HeterosynapticOutgoing <: HeterosynapticPlasticityTarget
  sum_max::Float64
end

struct PlasticityHeterosynapticEasy{HetMeth<:HeterosynapticPlasticityMethod,
    HetTarg<:HeterosynapticPlasticityTarget} <: PlasticityRule
  Δt_update::Float64
  method::HetMeth
  target::HetTarg
  _tcounter::Ref{Float64}
  function PlasticityHeterosynapticEasy(Δt_update::Float64,
      method::M,target::T) where {M<:HeterosynapticPlasticityMethod,T<:HeterosynapticPlasticityTarget}
    counter = Ref(0.0)
    new{M,T}(Δt_update,method,target,counter)
  end
end

function reset!(plast::PlasticityHeterosynapticEasy)
  plast._tcounter[] = zero(Float64)
  return nothing
end

function plasticity_update!(::R,dt::R,
     pspost::PopulationState,conn::Connection,pspre::PopulationState,
     plast::PlasticityHeterosynapticEasy) where R
  if plast._tcounter[] < plast.Δt_update
    plast._tcounter[] += dt
    return nothing
  end
  # reset timer
  plast._tcounter[] = zero(R)
  # apply plasticity at each col/row , depending on target
  _apply_easy_het_plasticity!(conn.weights,plast.method,plast.target)
  return nothing
end

@inline function hardbounds(x::Real,m::HeterosynapticPlasticityMethod)
  return  hardbounds(x,m.wmin,m.wmax)
end

# this should be a little faster than sum(M;dims=2)
function sum_over_rows!(dest::Vector,M::SparseMatrixCSC)
  _rowidx = SparseArrays.rowvals(M)
  Mnz = nonzeros(M) # direct access to weights
  @inbounds for (i,r) in enumerate(_rowidx)
    dest[r] += Mnz[i]
  end
  return dest
end

# how many elements are present along each col
function count_col_elements(M::SparseMatrixCSC)
  ncols = size(M,2)
  return length.(nzrange.(Ref(M),1:ncols))
end

# how many elements are present along each row
function count_row_elements(M::SparseMatrixCSC)
  nrows = size(M,1)
  ret = fill(0,nrows)
  for r in rowvals(M)
    ret[r]+=1
  end
  return ret
end

# additive, incoming connections corresponds to rows
function _apply_easy_het_plasticity!(sparseweights::SparseMatrixCSC,
    method::HeterosynapticAdditive,
    target::HeterosynapticIncoming)
  sum_max = target.sum_max
  nzw = nonzeros(sparseweights)
  wrows = rowvals(sparseweights)
  fixthing = Vector{Float64}(undef,size(sparseweights,1))
  sum_over_rows!(fixthing,sparseweights)
  n_el_row = count_row_elements(sparseweights)
  for i in eachindex(fixthing)
    _sum = fixthing[i]
    if _sum > sum_max
      fixthing[i] = (sum_max-_sum)/n_el_row[i]
    else
      fixthing[i] = 0.0
    end
  end
  for (_r,row) in enumerate(wrows)
    fix = fixthing[row]
    if !iszero(fix) # faster ?
      nzw[_r] = hardbounds(nzw[_r]+fix,method)
    end
  end
  return nothing
end

# additive, outgoing connections corresponds to columns
function _apply_easy_het_plasticity!(sparseweights::SparseMatrixCSC,
    method::HeterosynapticAdditive,
    target::HeterosynapticOutgoing)
  sum_max = target.sum_max
  nzw = nonzeros(sparseweights)
  sumswee = sum(sparseweights;dims=1)
  for (c,_sumr) in enumerate(sumswee)
    if _sumr > sum_max
      idxs = nzrange(sparseweights,c)
      for idx in idxs
        newval = nzw[idx] - (sumswee[c] - sum_max)/length(idxs)
        nzw[idx] = hardbounds(newval,method)
      end
    end
  end
  return nothing
end

# "precise" implementation... BUT...
# this code is just too slow :-(
struct PlasticityHeterosynaptic{HetMeth<:HeterosynapticPlasticityMethod,
    HetTarg<:HeterosynapticPlasticityTarget} <: PlasticityRule
  Δt_update::Float64
  method::HetMeth
  target::HetTarg
  _tcounter::Ref{Float64}
  function PlasticityHeterosynaptic(Δt_update::Float64,
      method::M,target::T) where {M<:HeterosynapticPlasticityMethod,T<:HeterosynapticPlasticityTarget}
    counter = Ref(0.0)
    new{M,T}(Δt_update,method,target,counter)
  end
end

function reset!(plast::PlasticityHeterosynaptic)
  plast._tcounter[] = zero(Float64)
  return nothing
end

function _get_row_idxs(M::SparseMatrixCSC,row_idx::Integer)
  return findall(==(row_idx), SparseArrays.rowvals(M))
end
function _get_col_idxs(M::SparseMatrixCSC,col_idx::Integer)
  cptr = SparseArrays.getcolptr(M)
  return collect(cptr[col_idx]:cptr[col_idx+1]-1)
end

# let's have an iterator to handle those two above, for each row/col
struct HeterosynapticIdxsIterator{HetTarg}
  M::SparseMatrixCSC
  target::HetTarg
end
Base.eltype(::Type{HeterosynapticIdxsIterator}) = Vector{Int64}
Base.getindex(hetit::HeterosynapticIdxsIterator,idx::Integer) = iterate(hetit,idx)[1]

# incoming connections: I am iterating over rows
Base.length(hetit::HeterosynapticIdxsIterator{HeterosynapticIncoming}) = size(hetit.M,1)
function Base.iterate(hetit::HeterosynapticIdxsIterator{HeterosynapticIncoming})
  idxs = _get_row_idxs(hetit.M,1)
  return (idxs,2)
end
function Base.iterate(hetit::HeterosynapticIdxsIterator{HeterosynapticIncoming},
    row::Int64)
  if row > length(hetit)
    return nothing
  else
    idxs = _get_row_idxs(hetit.M,row)
    row += 1
    return (idxs,row)
  end
end
# outgoing connections : iterate over columns
Base.length(hetit::HeterosynapticIdxsIterator{HeterosynapticOutgoing}) = size(hetit.M,2)
function Base.iterate(hetit::HeterosynapticIdxsIterator{HeterosynapticOutgoing})
  idxs = _get_col_idxs(hetit.M,1)
  return (idxs,2)
end
function Base.iterate(hetit::HeterosynapticIdxsIterator{HeterosynapticOutgoing},
    row::Int64)
  if row > length(hetit)
    return nothing
  else
    idxs = _get_col_idxs(hetit.M,row)
    row += 1
    return (idxs,row)
  end
end

# if values less than wmin, re-distributes the correction on the other neurons
function _heterosynaptic_fix!(nzvals::Vector{Float64},idxs::Vector{<:Integer},
    method::HeterosynapticAdditive,target::HeterosynapticPlasticityTarget)
  wmin = method.wmin
  nsyn = length(idxs)
  # fix without wmin
  fix_val = (sum(view(nzvals,idxs))-target.sum_max)/nsyn
  (fix_val < 0.0) && (return nothing)
  while true
    docycle = false
    accu = 0.0
    idx_low = Int64[]
    for idx in idxs
      nzvals[idx] -= fix_val
      if nzvals[idx] < wmin
        accu += wmin - nzvals[idx]
        nzvals[idx] = wmin
        push!(idx_low,idx)
        docycle=true
      end
    end
    if docycle
      filter!(!(in)(idx_low),idxs)
      fix_val = accu / (nsyn-length(idx_low))
      nsyn = length(idxs)
    else
      break
    end
  end
  return nothing
end

function _heterosynaptic_fix!(nzvals::Vector{Float64},idxs::Vector{<:Integer},
    method::HeterosynapticMultiplicative,target::HeterosynapticPlasticityTarget)
  wmin = method.wmin
  sum_targ = target.sum_max
  sum_val = sum(view(nzvals,idxs))
  fix_val = sum_targ/sum_val
  (fix_val > 1.0) && (return nothing)
  while true
    docycle = false
    sum_val = 0.0
    idx_low = Int64[]
    for idx in idxs
      nzvals[idx] *= fix_val # apply fix
      if nzvals[idx] < wmin  # check the fix
        nzvals[idx] = wmin
        push!(idx_low,idx)
        docycle = true
        sum_targ -= wmin  # update target sum to exclude min synapses
      else
        sum_val += nzvals[idx]  # update synapses sum
      end
    end
    if docycle
      filter!(!(in)(idx_low),idxs)
      fix_val = sum_targ/sum_val
    else
      break
    end
  end
  return nothing
end

function plasticity_update!(::Real,dt::Real,
     pspost::PopulationState,conn::Connection,pspre::PopulationState,
     plast::PlasticityHeterosynaptic)
  plast._tcounter[] += dt
  if plast._tcounter[] < plast.Δt_update
    return nothing
  end
  # reset timer
  plast._tcounter[] = zero(Float64)
  # apply plasticity at each col/row , depending on target
  wnzvals = nonzeros(conn.weights)
  for idxs in HeterosynapticIdxsIterator(conn.weights,plast.target)
    _heterosynaptic_fix!(wnzvals,idxs,plast.method,plast.target)
  end
  return nothing
end


text:
using Revise

using JDF

using CSV, DataFrames, Blosc, JLSO, Base.GC

# use 12 threads
Blosc.set_num_threads(6)
# GC.gc()
@time a = CSV.read(
    "C:/data/Performance_All/Performance_2016Q1.txt",
    delim = '|',
    header = false
);

GC.gc()

t = time()
@time metadatas = psavejdf(a, "c:/data/large.jdf");
time() - t
GC.gc()

GC.gc()
@time metadatas = savejdf(a, "c:/data/large.jdf");
GC.gc()

@time JLSO.save("c:/data/large.meta.jlso", metadatas)
GC.gc()

using Revise, JLSO, JDF
@time metadatas = JLSO.load("c:/data/large.meta.jlso")["data"];
GC.gc()

@time a2 = loadjdf("c:/data/large.jdf", metadatas);
GC.gc()

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))

# CSVFiles vs CSV
using Revise
using CSV, DataFrames, JDF, Base.GC
import Base.Threads: @spawn

# a = gf()
fp = [f for f in readdir("c:/data/AirOnTimeCSV/") if endswith(f, ".csv")];

using CSVFiles

@time a= load(joinpath("c:/data/AirOnTimeCSV", fp[1]), type_detect_rows = 2000)

#psavejdf(a, joinpath("c:/data/hehe/", fp[1]))
sort!(fp, by = x->filesize(joinpath("c:/data/AirOnTimeCSV/", x)), rev=true)

fn4(fp) = begin
    #@time a = load(joinpath("C:/data/AirOnTimeCSV", fp), type_detect_rows = 2000) |> DataFrame
    @time a = CSV.read(joinpath("C:/data/AirOnTimeCSV", fp))
    #mkdir(joinpath("c:/data/hehe/", fp))
    res = psavejdf(a, joinpath("c:/data/hehe/", fp))
    a = nothing
    res
end

res = Vector{Any}(undef, length(fp))

@time for (i, f) in enumerate(fp[1:6])
    print(i)
    res[i] = fn4(f)
end

#@time fn4(fp[1])
#@time fn4.(fp)
2+2

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

include("api_types.jl")
include("api_functions.jl")


text:
abstract AResult{T,N} <: AbstractArray{T,N}

immutable SimResult{T,N} <: AResult{T,N}
  sol::Array{T,N}

  t # todo: type this
  b::Array{T, 2}
end

immutable CumsimResult{T,N} <: AResult{T,N}
  sol::Array{T,N}
  t # todo: type this
end

Base.size(r::AResult) = size(r.sol)
Base.getindex(r::AResult, i...) = getindex(r.sol, i...)


text:
using YaoLang: YaoIR, is_pure_quantum
using ZXCalculus
using BenchmarkTools, JSON

function zx_load_qasm(filename::String)
    srcs = readlines("circuits/$(filename)")
    src = prod([srcs[1]; srcs[3:end]])
    # m = @__MODULE__
    # ir = YaoIR(m, src, :qasm_circ)
    # ir.pure_quantum = is_pure_quantum(ir)
    # zxd = ZXDiagram(ir)
    zxd = ZXDiagram(src, Val(:qasm))
    return zxd
end
function run_benchmark()
    filenames = readdir("circuits")
    benchmarks = Dict()
    for circ_name in filenames
        println("benchmarking $circ_name")
        zxd = zx_load_qasm(circ_name)
        b = @benchmark phase_teleportation($zxd)
        tc = tcount(phase_teleportation(zxd))
        # println(circ_name, "\ttime = ", (mean(b).time / 1e9), "(s)\tT-count = ", tc)
        benchmarks[circ_name] = Dict()
        benchmarks[circ_name]["times"] = Dict()
        benchmarks[circ_name]["times"]["mean"] = mean(b).time
        benchmarks[circ_name]["times"]["min"] = minimum(b).time
        benchmarks[circ_name]["times"]["max"] = maximum(b).time
        benchmarks[circ_name]["T-count"] = tc
    end
    return benchmarks
end

benchmarks = run_benchmark()

if !ispath("data")
    mkpath("data")
end

write("data/data.json", JSON.json(benchmarks))


text:
using Documenter, Low-Rank-Tensor-Approximation-Julia

makedocs(;
    modules=[Low-Rank-Tensor-Approximation-Julia],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/PierreSp/Low-Rank-Tensor-Approximation-Julia/blob/{commit}{path}#L{line}",
    sitename="Low-Rank-Tensor-Approximation-Julia.jl",
    authors="Pierre Springer",
    assets=String[],
)

deploydocs(;
    repo="github.com/PierreSp/Low-Rank-Tensor-Approximation-Julia",
)


text:
"""Throw `ArgumentError` when `\$x` is an immutable field."""
ExceptionImmutableField(x) = ArgumentError("`$x` is an immutable field.")

#= Pulse =#

export AbstractPulse
"""
An abstract type for pulses.

See also: [`Pulse`](@ref).
"""
abstract type AbstractPulse end

## Basics
Base.isequal(a::AbstractPulse, b::AbstractPulse) =
  all([isequal(getproperty(a,s), getproperty(b,s)) for s in fieldnames(Pulse)])

export Pulse
"""
A struct for typical MR pulses: `Pulse <: AbstractPulse`.

# Fields:
*Mutable*:
- `rf::TypeND(RF0D, [1,2])` (nT,) or (nT, nCoils).
- `gr::TypeND(GR0D, [2])` (nT, 3), where 3 accounts for x-y-z channels.
- `dt::T0D` (1,), simulation temporal step size, i.e., dwell time.
- `des::String`, an description of the pulse to be constructed.

# Usages:
    Pulse(rf, gr; dt=(4e-6)u"s", des="generic pulse")

See also: [`AbstractPulse`](@ref).
"""
mutable struct Pulse <: AbstractPulse
  rf::TypeND(RF0D, [1,2])
  gr::TypeND(GR0D, [2])
  dt::T0D
  des::String

  function Pulse(rf=missing, gr=missing; dt=4e-6u"s", des="generic pulse")
    rf_miss, gr_miss = ismissing(rf), ismissing(gr)
    rf_miss&&gr_miss && ErrorException("Missing both inputs.")
    if rf_miss rf = zeros(size(gr,1))u"Gauss" end
    if gr_miss gr = zeros(size(rf,1),3)u"Gauss/cm" end
    size(gr,2)==3 || throw(DimensionMismatch)

    if isa(rf, Number) rf = [rf] end
    return new(rf, gr, dt, des)
  end

end

## set and get
Base.setproperty!(p::Pulse, sym::Symbol, x) = begin
  if (sym==:gr) size(x)==(size(p.rf,1),3)||throw(DimensionMismatch) end
  if (sym==:rf) size(x,1)==size(p.gr,1)||throw(DimensionMismatch) end
  setfield!(p, sym, x)
end

#= Spin =#

export AbstractSpinArray
"""
This type keeps the essentials of magnetic spins. Its instance struct must
contain all fields listed listed in the exemplary struct `mSpinArray`.

# Misc
Might make `AbstractSpinArray <: AbstractArray` in a future version

See also: [`mSpinArray`](@ref), [`AbstractSpinCube`](@ref).
"""
abstract type AbstractSpinArray end

## set and get
Base.setproperty!(spa::AbstractSpinArray, s::Symbol, x) = begin
  s ∈ (:dim, :mask) && throw(ExceptionImmutableField(s))

  nM = prod(spa.dim)
  if (s==:M)&&(size(x,1)==1) x=repeat(x, nM, 1) end
  if (s ∈ (:T1,:T2,:γ,:M)) size(x,1)∈(1,nM)||throw(DimensionMismatch) end

  setfield!(spa, s, x)
end

## AbstractArray-like interface
Base.size(spa::AbstractSpinArray) = spa.dim
Base.size(spa::AbstractSpinArray, d) = (d ≤ length(spa.dim)) ? spa.dim[d] : 1

Base.isequal(a::AbstractSpinArray, b::AbstractSpinArray) =
  all([isequal(getproperty.((a,b),s)...) for s in fieldnames(mSpinArray)])

## Concrete mSpinArray
export mSpinArray
"""
    mSpinArray(dim::Dims; T1=1.47u"s", T2=0.07u"s", γ=γ¹H, M=[0. 0. 1.])
    mSpinArray(mask::BitArray; T1=1.47u"s", T2=0.07u"s", γ=γ¹H, M=[0. 0. 1.])
An exemplary struct instantiating `AbstractSpinArray`.

# Fields:
*Immutable*:
- `dim::Dims` (nd,): `nM ← prod(dim)`, dimension of the object.
- `mask::BitArray` (nx,(ny,(nz))): Mask for `M`, `dim == (nx,(ny,(nz)))`
*Mutable*:
- `T1::TypeND(T0D, [0,1])` (1,) or (nM,): Longitudinal relaxation coeff.
- `T2::TypeND(T0D, [0,1])` (1,) or (nM,): Transversal relaxation coeff.
- `γ::TypeND(Γ0D, [0,1])`  (1,) or (nM,): Gyromagnetic ratio.
- `M::TypeND(Real, [2])`   (`count(mask)`, 3):  Magnetic spins, (𝑀x,𝑀y,𝑀z).

# Notes:
off-resonance, `Δf`, and locations, `loc`, are intentionally unincluded, as they
are not intrinsic to spins, and can change over time. Unincluding them allows
extensional subtypes specialized for, e.g., arterial spin labelling.

See also: [`AbstractSpinArray`](@ref).
"""
mutable struct mSpinArray <: AbstractSpinArray
  # *Immutable*:
  dim ::Dims
  mask::BitArray
  # *Mutable*:
  T1 ::TypeND(T0D, [0,1])
  T2 ::TypeND(T0D, [0,1])
  γ  ::TypeND(Γ0D, [0,1])
  M  ::TypeND(AbstractFloat, [2])

  function mSpinArray(mask::BitArray;
                      T1=1.47u"s", T2=0.07u"s", γ=γ¹H, M=[0. 0. 1.])
    dim = size(mask)
    nM = prod(dim)
    M = eltype(M)<:AbstractFloat ? M : float(M)
    if size(M,1)==1 M=repeat(M, nM, 1) end
    (all(map(x->(size(x,1) ∈ (1,nM)), [T1,T2,γ,M]))) || throw(DimensionMismatch)

    return new(dim, mask, T1, T2, γ, M)
  end

  mSpinArray(dim::Dims=(1,); kw...) = mSpinArray(trues(dim); kw...)

end

#= Cube =#

export AbstractSpinCube
"""
`AbstractSpinCube <: AbstractSpinArray`.
This type inherits `AbstractSpinArray` as a field. Its instance struct must
contain all fields listed in the exemplary struct `mSpinCube`.

See also: [`AbstractSpinArray`](@ref), [`mSpinCube`](@ref).
"""
abstract type AbstractSpinCube <: AbstractSpinArray end

## set and get
Base.setproperty!(cb::AbstractSpinCube, s::Symbol, x) = begin
  s ∈ (:spinarray,:fov,:ofst,:loc) && throw(ExceptionImmutableField(s))
  s ∈ fieldnames(typeof(cb)) ? setfield!(cb, s,x) : setfield!(cb.spinarray, s,x)
end

Base.getproperty(cb::AbstractSpinCube, s::Symbol) =
  s ∈ fieldnames(typeof(cb)) ? getfield(cb, s) : getfield(cb.spinarray, s)

## AbstractArray-like interface
Base.size(cb::AbstractSpinCube, a...) = Base.size(cb.spinarray, a...)

Base.isequal(a::AbstractSpinCube, b::AbstractSpinCube) =
  all([isequal(getproperty.((a,b),s)...) for s in fieldnames(mSpinCube)])

## Concrete mSpinCube
export mSpinCube
"""
    spincube = mSpinCube(dim::Dims{3}, fov; ofst, Δf, T1, T2, γ)
    spincube = mSpinCube(mask::BitArray{3}, fov; ofst, Δf, T1, T2, γ)
`dim`, `mask`, `T1`, `T2`, and `γ` are passed to `mSpinArray` constructors.

An exemplary struct instantiating `AbstractSpinCube`, designed to model a set of
regularly spaced spins, e.g., a volume.

# Fields:
*Immutable*:
- `spinarray::AbstractSpinArray` (1,): inherited `AbstractSpinArray` struct
- `fov ::TypeND(L0D, [2])` (1, 3): field of view.
- `ofst::TypeND(L0D, [2])` (1, 3): fov offset from magnetic field iso-center.
- `loc ::TypeND(L0D, [2])` (nM, 3): location of spins.
*Mutable*:
- `Δf::TypeND(F0D, [0,1])` (1,) or (nM,): off-resonance map.

See also: [`AbstractSpinCube`](@ref).
"""
mutable struct mSpinCube <: AbstractSpinCube
  # *Immutable*:
  spinarray::AbstractSpinArray
  fov ::TypeND(L0D, [2])
  ofst::TypeND(L0D, [2])
  loc ::TypeND(L0D, [2])
  # *Mutable*:
  Δf  ::TypeND(F0D, [0,1])

  function mSpinCube(mask::BitArray{3}, fov::TypeND(L0D, [2]);
                     ofst::TypeND(L0D, [2])=[0 0 0]u"cm", Δf=0u"Hz",
                     T1=1.47u"s", T2=0.07u"s", γ=γ¹H)
    size(fov)==size(ofst)==(1,3) || throw(DimensionMismatch)
    spa = mSpinArray(mask; T1=T1, T2=T2, γ=γ)
    loc = CartesianLocations(spa.dim)./(reshape([spa.dim...], 1,:)./fov) .+ ofst
    return new(spa, fov, ofst, loc, Δf)
  end

  mSpinCube(dim::Dims{3}, a...; kw...) = mSpinCube(trues(dim), a...; kw...)

end

#= Bolus (*Under Construction*) =#
# TODO
# export AbstractSpinBolus
"""
*UNDER CONSTRUCTION*

`AbstractSpinBolus <: AbstractSpinArray`.
This type inherits `AbstractSpinArray` as a field. Its instance struct must
contain all fields listed in the exemplary struct `mSpinBolus`.

See also: [`AbstractSpinArray`](@ref), [`mSpinBolus`](@ref).
"""
abstract type AbstractSpinBolus <: AbstractSpinArray end

## set and get

Base.getproperty(bl::AbstractSpinBolus, s::Symbol) =
  s ∈ fieldnames(typeof(bl)) ? getfield(bl, s) : getfield(bl.spinarray, s)

## AbstractArray-like interface
Base.size(bl::AbstractSpinBolus, a...) = Base.size(bl.spinarray, a...)

Base.isequal(a::AbstractSpinBolus, b::AbstractSpinBolus) =
  all([isequal(getproperty.((a,b),s)...) for s in fieldnames(mSpinBolus)])

## Concrete mSpinBolus
# export mSpinBolus
"""
*UNDER CONSTRUCTION*

An exemplary struct instantiating `AbstractSpinBolus`, designed to model a set
of moving spins, e.g., a blood bolus in ASL context.

See also: [`AbstractSpinBolus`](@ref).
"""
mutable struct mSpinBolus <: AbstractSpinBolus
  # *Immutable*:
  # *Mutable*:
end

#= mObjects utils =#

export Pulse2B
"""
    B = Pulse2B(pulse::Pulse, loc; Δf, b1Map, γ)
Create effective magnetic field, 𝐵, from input `pulse`.

See also: [`rfgr2B`](@ref), [`B2UΦ`](@ref), [`blochSim`](@ref).
"""
Pulse2B(p::Pulse, loc; kw...) = rfgr2B(p.rf, p.gr, loc; kw...)

"""
    B = Pulse2B(pulse::Pulse, spa::AbstractSpinArray, loc; Δf, b1Map)
...with `γ=spa.γ`.
"""
Pulse2B(p::Pulse, spa::AbstractSpinArray, loc; kw...) =
  Pulse2B(p, loc; γ=spa.γ, kw...)

"""
    B = Pulse2B(pulse::Pulse, cb::AbstractSpinCube; b1Map)
...with `loc, Δf, γ = cb.loc, cb.Δf, cb.γ`.
"""
Pulse2B(p::Pulse, cb::AbstractSpinCube; kw...) =
  Pulse2B(p, cb.loc, cb.Δf; γ=cb.γ, kw...)

export applyPulse, applyPulse!
"""
    applyPulse(spa::AbstractSpinArray, p::Pulse, loc; Δf, b1Map, doHist)
Turn `p` into 𝐵-effective and apply it on `spa.M`, using its own `M, T1, T2, γ`.

See also: [`blochSim`](@ref), [`freePrec`](@ref).
"""
applyPulse(spa::AbstractSpinArray, p::Pulse, loc; doHist=false, kw...) =
  blochSim(spa.M, Pulse2B(p, spa, loc; kw...);
           T1=spa.T1, T2=spa.T2, γ=spa.γ, dt=p.dt, doHist=doHist)

"""
    applyPulse!(spa::AbstractSpinArray, p::Pulse, loc; Δf, b1Map, doHist)
Update `spa.M` before return.
"""
applyPulse!(spa::AbstractSpinArray, p::Pulse, loc; doHist=false, kw...) = begin
  _, Mhst = blochSim!(spa.M, Pulse2B(p, spa, loc; kw...);
                      T1=spa.T1, T2=spa.T2, γ=spa.γ, dt=p.dt, doHist=doHist)
  return (M=spa.M, Mhst=Mhst)
end

"""
    applyPulse(cb::AbstractSpinCube, p::Pulse; b1Map, doHist)
Turn `p` into 𝐵-effective and apply it on `cb.M`, using its own `M, T1, T2, γ`.
"""
applyPulse(cb::AbstractSpinCube, p::Pulse; b1Map=1, doHist=false) =
  blochSim(cb.M, Pulse2B(p, cb; b1Map=b1Map);
           T1=cb.T1, T2=cb.T2, γ=cb.γ, dt=p.dt, doHist=doHist)

"""
    applyPulse!(cb::AbstractSpinCube, p::Pulse; b1Map, doHist)
Update `cb.M` before return.
"""
applyPulse!(cb::AbstractSpinCube, p::Pulse; b1Map=1, doHist=false) = begin
  _, Mhst = blochSim!(cb.M, Pulse2B(p, cb; b1Map=b1Map);
                      T1=cb.T1, T2=cb.T2, γ=cb.γ, dt=p.dt, doHist=doHist)
  return (M=cb.M, Mhst=Mhst)
end

export freePrec!, freePrec
"""
    freePrec(spa::AbstractSpinArray, t; Δf)
`spa::AbstractSpinArray` free precess by `t`. `spa.M` will not be updated.

See also: [`applyPulse`](@ref), [`freePrec!`](@ref).
"""
freePrec(spa::AbstractSpinArray, t::T0D; Δf::TypeND(F0D,[0,1])=0u"Hz") =
  freePrec(spa.M, t; Δf=Δf, T1=spa.T1, T2=spa.T2)

"""
    freePrec!(spa::AbstractSpinArray, t; Δf)
...`spa.M` will updated by the results.
"""
freePrec!(spa::AbstractSpinArray, t::T0D; Δf::TypeND(F0D,[0,1])=0u"Hz") =
  freePrec!(spa.M, t; Δf=Δf, T1=spa.T1, T2=spa.T2)

"""
    freePrec(cb::AbstractSpinCube, t)
`cb::AbstractSpinCube` free precess by `t`. `cb.M` will not be updated.

See also: [`applyPulse`](@ref), [`freePrec`](@ref).
"""
freePrec(cb::AbstractSpinCube, t::T0D) =
  freePrec(cb.M, t; Δf=cb.Δf, T1=cb.T1, T2=cb.T2)

"""
    freePrec!(cb::AbstractSpinCube, t)
...`cb.M` will be updated by the results.
"""
freePrec!(cb::AbstractSpinCube, t::T0D) =
  freePrec!(cb.M, t; Δf=cb.Δf, T1=cb.T1, T2=cb.T2)


text:
"""
    Genser.@fieldtype(type, fieldkey, gensertype)

Override the genser type for a struct field.

```@example
using Genser

struct Person
    id::Base.UUID
    name::String
end

# Always convert the UUID to a string
@fieldtype Persion :id GenserString
```
"""
macro fieldtype(type, field, gensertype)
    type = esc(type)
    # gensertype = esc(gensertype)
    return quote
        Genser.fieldtype(::Type{$type}, ::Type{Val{$field}}) = $gensertype
    end
end

"""
    Genser.@fieldencoding(type, fieldkey, encodingkey)

Set the string encoding for a binary value.

```@example
using Genser

struct DataContainer
    contents::Vector{UInt8}
end

# Encode as base64 when serialized as a string
@fieldencoding DataContainer :contents :base64
```
"""
macro fieldencoding(type, field, encoding)
    type = esc(type)
    return quote
        Genser.fieldencoding(::Type{$type}, ::Type{Val{$field}}) = Encoding{$encoding}
    end
end


text:
module LongestCommonSubsequence

# ----------------------------------------
# EXPORTED INTERFACE
# ----------------------------------------

export lcs

# ----------------------------------------
# EXPORTED FUNCTIONS
# ----------------------------------------

function lcs(left, right, comparator)
    splits = split(walk(table(left, right, comparator), left, right, comparator))

    grouper   = ( x, y ) -> y[1] == x[1] + 1 && x[2] == y[2]
    collector = group -> ( group[1][2], range(group[1][1], group[end][1], step = 1) )

    return map(groups -> map(collector, groups), map(list -> group_by(list, grouper), splits))
end

# ----------------------------------------
# PRIVATE FUNCTIONS
# ----------------------------------------

# Group up the sequence pairs ( boolean, index ) into a sequence that uses ranges ... ( boolean, [start ... end] ).
#
function group_by(ungrouped, comparator)
    function worker(head, tail, accumulator, accumulators)
        push!(accumulator, head)
        if isempty(tail)
            push!(accumulators, accumulator)

        else
            if !comparator(head, tail[1])
                push!(accumulators, accumulator)
                accumulator = [ ]
            end

            worker(tail[1], tail[2:end], accumulator, accumulators)
        end
    end

    accumulators = [ ]

    if !isempty(ungrouped)
        worker(ungrouped[1], ungrouped[2:end], [ ], accumulators)
    end

    return accumulators
end

# This builds the longest-common-subsequence table of costs.
#
function table(left, right, comparator)
    rows    = length(left) + 1
    columns = length(right) + 1
    costs   = zeros(UInt64, rows, columns)

    for column = 2:columns
        for row = 2:rows
            if comparator(left[row - 1], right[column - 1])
                costs[row, column] = costs[row - 1, column - 1] + 1

            else
                costs[row, column] = max(costs[row, column - 1], costs[row - 1, column])
            end
        end
    end

    return costs
end

# Partition the sets into left and right groupings.
#
function split(boths)
    return ( map(left -> ( left[2], left[1] == "L" ), filter(both -> both[1] == "L" || both[1] == "=", boths)),
             map(right -> ( right[3], right[1] == "R" ), filter(both -> both[1] == "R" || both[1] == "=", boths)) )
end

# Walk the LCS costs table to generate diff op-codes.
#
function walk(costs, left, right, comparator)
    function worker(i, j, accumulator::Array{Any, 1})
        if i > 1 && j > 1 && comparator(left[i - 1], right[j - 1])
            worker(i - 1, j - 1, accumulator)
            push!(accumulator, ( "=", i - 1, j - 1, left[i - 1] ))

        elseif j > 1 && (i == 1 || costs[i, j - 1] >= costs[i - 1, j])
            worker(i, j - 1, accumulator)
            push!(accumulator, ( "R", i - 1, j - 1, right[j - 1] ))

        elseif i > 1 && (j == 1 || costs[i, j - 1] < costs[i - 1, j])
            worker(i - 1, j, accumulator)
            push!(accumulator, ( "L", i - 1, j - 1, left[i - 1] ))
        end
    end

    accumulator = [ ]
    worker(size(costs)[1], size(costs)[2], accumulator)

    return accumulator
end

end

max_stars_repo_name: jlegare/Longest-Common-Subsequence

max_issues_repo_name: jlegare/Longest-Common-Subsequence

max_forks_repo_name: jlegare/Longest-Common-Subsequence


text:
@testset "example" begin

my_graph = "graph {
  1 -- 2
  2 -- 3
  3 -- 1
}
"

root = parse_dot(my_graph)[1]

for node in nodes(root)
    println("node $(node)")
end
for (node1, node2) in edges(root)
    println("edge $(node1) - $(node2)")
end

end


text:
using Test
using SetOps

A = Set(1:5)
B = Set(3:6)

@test A+B == Set(1:6)
@test A|B == A+B
@test A&B == Set(3:5)
@test B-A == Set(6)
@test A$B == Set([1,2,6])
@test length(A*B) == length(A)*length(B)
@test length(powerset(B)) == 16


text: # This is a comment in Julia


text:
function cal_flux_residual(loc_globaldata, loc_ghost_holder, power, limiter_flag, vl_const, gamma)
    loc_globaldata = loc_globaldata[:L]
    loc_ghost_holder = loc_ghost_holder[:L]

    # main_store = zeros(MVector{4})
    phi_i = zeros(MVector{4})
	phi_k = zeros(MVector{4})
	G_i = zeros(MVector{4})
    G_k = zeros(MVector{4})
	result = zeros(MVector{4})
	qtilde_i = zeros(MVector{4})
	qtilde_k = zeros(MVector{4})
	Gxp = zeros(MVector{4})
	Gxn = zeros(MVector{4})
	Gyp = zeros(MVector{4})
	Gyn = zeros(MVector{4})
    ∑_Δx_Δf = zeros(MVector{4})
    ∑_Δy_Δf = zeros(MVector{4})

    dist_length = length(loc_globaldata)

	for idx in 1:dist_length
		if loc_globaldata[idx].flag_1 == 0
			wallindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const )
		elseif loc_globaldata[idx].flag_1 == 2
			outerindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const )
		elseif loc_globaldata[idx].flag_1 == 1
			interiorindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const)
		end
	end
	return nothing
end

function wallindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const )
		wall_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxp )
		wall_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxn )
		wall_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gyn )
		loc_globaldata[idx].flux_res = SVector{4}((Gxp + Gxn + Gyn) * 2)

	return nothing
end

function outerindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const )
	outer_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxp )
	outer_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxn )
	outer_dGy_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gyp )
	loc_globaldata[idx].flux_res = SVector{4}(Gxp + Gxn + Gyp)
	return nothing
end

function interiorindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, gamma, idx, Gxp, Gxn, Gyp, Gyn, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const)
	interior_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxp)
	interior_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gxn)
	interior_dGy_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gyp)
	interior_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, gamma, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k, ∑_Δx_Δf, ∑_Δy_Δf, power, limiter_flag, vl_const, Gyn)
	loc_globaldata[idx].flux_res = SVector{4}(Gxp + Gxn + Gyp + Gyn)
	return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
promote_rule(::Type{S}, ::Type{T}) where S<:SafeSigned   where T<:SafeSigned =
    safeint(promote_type(baseint(S), baseint(T)))
promote_rule(::Type{S}, ::Type{T}) where S<:SafeSigned   where T<:SafeUnsigned =
    safeint(promote_type(baseint(S), baseint(T)))
promote_rule(::Type{S}, ::Type{T}) where S<:SafeUnsigned where T<:SafeUnsigned =
    safeint(promote_type(baseint(S), baseint(T)))

promote_rule(::Type{S}, ::Type{T}) where S<:SafeSigned   where T<:Signed   = S
promote_rule(::Type{S}, ::Type{T}) where S<:SafeSigned   where T<:Unsigned = S
promote_rule(::Type{S}, ::Type{T}) where S<:SafeUnsigned where T<:Signed   = S
promote_rule(::Type{S}, ::Type{T}) where S<:SafeUnsigned where T<:Unsigned = S

promote_rule(::Type{S}, ::Type{T}) where S<:SafeSigned   where T<:Base.IEEEFloat = T
promote_rule(::Type{S}, ::Type{T}) where S<:SafeUnsigned where T<:Base.IEEEFloat = T

promote_rule(::Type{S}, ::Type{BigFloat}) where S<:SafeSigned  = BigFloat
promote_rule(::Type{S}, ::Type{BigFloat}) where S<:SafeUnsigned = BigFloat