text:
#' plots.jl: make plots for the paper
#' ================================================

#' To generate a document from the code output, run the following at the Julia terminal:
#' `using Weave` then
#' `weave("scripts/plots.jl", doctype = "github", out_path = "scripts")`

#= ------------ for loading a simulation run from a remote server
local_filename   = "out_figure4.jld"
remote_filename  = "out_figure4.jld"
remote_machine   =  "poisson" #
@eval run(`
    scp
    anderes@$remote_machine.ucdavis.edu:~/PrimordialB.jl/scripts/$remote_filename
    /Users/ethananderes/Dropbox/PrimordialB/scripts/$local_filename
`)
=# #---------------------------------------------

#   .oooooo.                   o8o            oooo                        oooo                .
#  d8P'  `Y8b                  `"'            `888                        `888              .o8
# 888      888    oooo  oooo  oooo   .ooooo.   888  oooo       oo.ooooo.   888   .ooooo.  .o888oo
# 888      888    `888  `888  `888  d88' `"Y8  888 .8P'         888' `88b  888  d88' `88b   888
# 888      888     888   888   888  888        888888.          888   888  888  888   888   888
# `88b    d88b     888   888   888  888   .o8  888 `88b.        888   888  888  888   888   888 .
#  `Y8bood8P'Ybd'  `V88V"V8P' o888o `Y8bod8P' o888o o888o       888bod8P' o888o `Y8bod8P'   "888"
#                                                               888
#                                                              o888o
#
#

#= Quick plot for viewing results

using JLD, PyPlot
loadpath = "scripts/saved_run_0x3f77869c073f2346/out_figure2c.jld"
rng_true = load(loadpath, "rng_true")
rest_B_1 = load(loadpath, "rest_B")
rest_1   = load(loadpath, "rest_modes")
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)
standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))

fig, ax = plt[:subplots](1,1, figsize=(5,3.7), sharex=true,  sharey=true)
fontsize = 10

ax[:plot](rng_true, rng_true./standard_error_B_1, label = "quadratic delenser")
ax[:plot](rng_true, rng_true./standard_error_L_1, label = "likelihood approx")
ax[:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:set_title](L"\mu K arcmin = 1.0", fontsize=fontsize)

fig, ax = plt[:subplots](1,1, figsize=(5,3.7), sharex=true,  sharey=true)
fontsize = 10

ax[:plot](rng_true, bias_B_1./rng_true, label = "quadratic delenser")
ax[:plot](rng_true, bias_L_1./rng_true, label = "likelihood approx")
ax[:set_ylabel](L"$\frac{Bias(r)}{r}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:set_title](L"\mu K arcmin = 1.0", fontsize=fontsize)
ax[:axhline](y=.0, color = "k", linestyle="dashed")

=#

#  .o88o.  o8o                                                   .o
#  888 `"  `"'                                                 o888
# o888oo  oooo   .oooooooo oooo  oooo  oooo d8b  .ooooo.        888
#  888    `888  888' `88b  `888  `888  `888""8P d88' `88b       888
#  888     888  888   888   888   888   888     888ooo888       888
#  888     888  `88bod8P'   888   888   888     888    .o       888
# o888o   o888o `8oooooo.   `V88V"V8P' d888b    `Y8bod8P'      o888o
#               d"     YD
#               "Y88888P'
#
#
#
# Plot of phi estimation noise used in the simulations

using PyPlot, JLD
savefigs = true
savepath = "scripts/saved_run_0x3f77869c073f2346/"

#' f1.pdf : phi N0ℓ for CMB-S3 vrs CMB-S4
#' ================================================
include("src/PrimordialB.jl")
kbins, N0ℓ_a, cϕϕ = let
    nside     = load(savepath*"out_figure2a.jld", "nside")
    period    = load(savepath*"out_figure2a.jld", "period")
    μKarcminT = load(savepath*"out_figure2a.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2a.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2a.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end
kbins, N0ℓ_b, cϕϕ = let
    nside     = load(savepath*"out_figure2b.jld", "nside")
    period    = load(savepath*"out_figure2b.jld", "period")
    μKarcminT = load(savepath*"out_figure2b.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2b.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2b.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end
kbins, N0ℓ_c, cϕϕ = let
    nside     = load(savepath*"out_figure2c.jld", "nside")
    period    = load(savepath*"out_figure2c.jld", "period")
    μKarcminT = load(savepath*"out_figure2c.jld", "μKarcminT")
    pixel_side_length = load(savepath*"out_figure2c.jld", "pixel_side_length")

    phiest_pixel_side_length = 2 # fixed
    phiest_nside          = ceil(Int64, pixel_side_length*nside/phiest_pixel_side_length) #nside * 5
    phiest_μKarcminT      = μKarcminT
    beamFWHM              = deg2rad(pixel_side_length/60)
    phiest_beamFWHM       = deg2rad(phiest_pixel_side_length/60)

    N0ℓ       = load(savepath*"out_figure2c.jld", "N0ℓ")
    g = PrimordialB.FFTgrid(PrimordialB.Pix{period, nside}())
    cls, mCls, cls_no_B    = PrimordialB.get_spectra(μKarcminT, beamFWHM, period, nside) # for B survey

    kbins, N0ℓ_tmp  = PrimordialB.radial_ave((g.r.^4) .* N0ℓ, 1, g)
    _, cϕϕ    = PrimordialB.radial_ave((g.r.^4) .* mCls.cϕϕk, 1, g)
    kbins, N0ℓ_tmp, cϕϕ
end

#' Plot
fig, ax = plt[:subplots](1,1, figsize=(4.5,3.5))
fontsize = 10
ax[:plot](kbins, cϕϕ, label = L"$\ell^4 C_\ell^{\phi\phi}$")
ax[:plot](kbins, N0ℓ_a, label = L"\mu K arcmin = 10")
ax[:plot](kbins, N0ℓ_b, label = L"\mu K arcmin = 1")
ax[:plot](kbins, N0ℓ_c, label = L"\mu K arcmin = 0.5")
ax[:set_xscale]("log")
ax[:set_yscale]("log")
ax[:set_title](L"$\phi$ uncertainty", fontsize=fontsize)
ax[:set_xlabel](L"multiple $\ell$", fontsize=fontsize)
ax[:legend](loc="best", fontsize=fontsize-2)
#ax[:set_ylim]([1e-8, 1.5e-6])
ax[:grid](true)

#' Now save
savename  = "f1.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

# using PyPlot
# loglog(cls[:ell], cls[:ell].^2 .* cls[:tt])
# loglog(cls[:ell], cls[:ell].^2 .* cls[:ee])
# loglog(cls[:ell], cls[:ell].^2 .* cls[:bb])
# #loglog(cls[:ell], cls[:ell].^2 .* cls_no_B[:bb])
# loglog(cls[:ell], cls[:ell].^2 .* cls_no_B[:ln_bb])
# ksl = g.k[2][1:25,1]
# loglog(ksl, ksl.^2 .* mCls.cBBnoisek[1:25,1])
# loglog(ksl, ksl.^2 .* mCls.cBBk[1:25,1])
# loglog(ksl, ksl.^2 .* mCls.cEEk[1:25,1])
# loglog(ksl, ksl.^2 .* mCls.cTTk[1:25,1])

#  .o88o.  o8o                                                   .oooo.
#  888 `"  `"'                                                 .dP""Y88b
# o888oo  oooo   .oooooooo oooo  oooo  oooo d8b  .ooooo.             ]8P'
#  888    `888  888' `88b  `888  `888  `888""8P d88' `88b          .d8P'
#  888     888  888   888   888   888   888     888ooo888        .dP'
#  888     888  `88bod8P'   888   888   888     888    .o      .oP     .o
# o888o   o888o `8oooooo.   `V88V"V8P' d888b    `Y8bod8P'      8888888888
#               d"     YD
#               "Y88888P'
#
#
#
#

using PyPlot, JLD
savefigs = true
savepath = "scripts/saved_run_0x3f77869c073f2346/"

#' f2a.pdf, f2b.pdf, f2c.pdf:A comparison of
#' ================================================

rng_true = load(savepath*"out_figure2a.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure2a.jld", "rest_B")
rest_1   = load(savepath*"out_figure2a.jld", "rest_modes")
rest_B_2 = load(savepath*"out_figure2b.jld", "rest_B")
rest_2   = load(savepath*"out_figure2b.jld", "rest_modes")
rest_B_3 = load(savepath*"out_figure2c.jld", "rest_B")
rest_3   = load(savepath*"out_figure2c.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_B_2  = mean((rest_B_2 .- rng_true),2)
bias_B_3  = mean((rest_B_3 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)
bias_L_2  = mean((rest_2 .- rng_true),2)
bias_L_3  = mean((rest_3 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_B_2 = sqrt.(mean((rest_B_2 .- rng_true).^2, 2))
standard_error_B_3 = sqrt.(mean((rest_B_3 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))
standard_error_L_2 = sqrt.(mean((rest_2 .- rng_true).^2, 2))
standard_error_L_3 = sqrt.(mean((rest_3 .- rng_true).^2, 2))

maxylab = maximum(vcat(
    rng_true./standard_error_B_1,
    rng_true./standard_error_B_2,
    rng_true./standard_error_B_3,
    rng_true./standard_error_L_1,
    rng_true./standard_error_L_2,
    rng_true./standard_error_L_3
    ))

minbias, maxbias = vcat(
         bias_B_1./standard_error_B_1,
         bias_B_2./standard_error_B_2,
         bias_B_3./standard_error_B_3,
         bias_L_1./standard_error_L_1,
         bias_L_2./standard_error_L_2,
         bias_L_3./standard_error_L_3,
         ) |> x-> (minimum(x), maximum(x))

#' f2a.pdf : plot r / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,3, figsize=(10,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[1][:set_title](L"\mu K arcmin = 10", fontsize=fontsize)
ax[1][:plot](rng_true, rng_true./standard_error_B_1, label = "quadratic delenser")
ax[1][:plot](rng_true, rng_true./standard_error_L_1, label = "likelihood approx")
ax[1][:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[1][:set_xlabel](L"$r$", fontsize=fontsize)
ax[1][:set_xticks]([.01, .03, .05, .07, .09])
ax[1][:tick_params](axis="y", labelsize = fontsize-2)
ax[1][:tick_params](axis="x", labelsize = fontsize-2)
ax[1][:legend](loc="best", fontsize=fontsize)
ax[1][:grid](true)

#  CMB-S4 for ϕ est
ax[2][:set_title](L"\mu K arcmin = 1", fontsize=fontsize)
ax[2][:plot](rng_true, rng_true./standard_error_B_2, label = "quadratic delenser")
ax[2][:plot](rng_true, rng_true./standard_error_L_2, label = "likelihood approx")
ax[2][:set_xlabel](L"$r$", fontsize=fontsize)
ax[2][:tick_params](axis="y", left="off", right="off")#labelsize = fontsize-2, pad = 5, labelright="on")
ax[2][:tick_params](axis="x", labelsize = fontsize-2)
ax[2][:axis]("tight")
ax[2][:grid](true)

#  CMB-S4 for ϕ est
ax[3][:set_title](L"\mu K arcmin = 0.5", fontsize=fontsize)
ax[3][:plot](rng_true, rng_true./standard_error_B_3, label = "quadratic delenser")
ax[3][:plot](rng_true, rng_true./standard_error_L_3, label = "likelihood approx")
ax[3][:set_xlabel](L"$r$", fontsize=fontsize)
ax[3][:tick_params](axis="y", labelsize = fontsize-2, labelright="on")
ax[3][:tick_params](axis="x", labelsize = fontsize-2, pad = 5)
ax[3][:axis]("tight")
ax[3][:grid](true)

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f2a.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f2b.pdf : plot  Bias(r) / r
#' ================================================

fig, ax = plt[:subplots](1,3, figsize=(10,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[1][:set_title](L"\mu K arcmin = 10", fontsize=fontsize)
ax[1][:plot](rng_true, bias_B_1./rng_true, label = "quadratic delenser")
ax[1][:plot](rng_true, bias_L_1./rng_true, label = "likelihood approx")
ax[1][:set_ylabel](L"$\frac{Bias(r)}{r}$        ", fontsize=fontsize+4,  rotation="horizontal")
ax[1][:set_xlabel](L"$r$", fontsize=fontsize)
ax[1][:set_xticks]([.01, .03, .05, .07, .09])
ax[1][:tick_params](axis="y", labelsize = fontsize-2)
ax[1][:tick_params](axis="x", labelsize = fontsize-2)
ax[1][:legend](loc="best", fontsize=fontsize)
ax[1][:grid](true)
ax[1][:axhline](y=.0, color = "k", linestyle="dashed")

#  CMB-S4 for ϕ est
ax[2][:set_title](L"\mu K arcmin = 1", fontsize=fontsize)
ax[2][:plot](rng_true, bias_B_2./rng_true, label = "quadratic delenser")
ax[2][:plot](rng_true, bias_L_2./rng_true, label = "likelihood approx")
ax[2][:set_xlabel](L"$r$", fontsize=fontsize)
ax[2][:tick_params](axis="y", left="off", right="off")
ax[2][:tick_params](axis="x", labelsize = fontsize-2)
ax[2][:axis]("tight")
ax[2][:grid](true)
ax[2][:axhline](y=.0, color = "k", linestyle="dashed")

#  CMB-S4 for ϕ est
ax[3][:set_title](L"\mu K arcmin = 0.5", fontsize=fontsize)
ax[3][:plot](rng_true, bias_B_3./rng_true, label = "quadratic delenser")
ax[3][:plot](rng_true, bias_L_3./rng_true, label = "likelihood approx")
ax[3][:set_xlabel](L"$r$", fontsize=fontsize)
ax[3][:tick_params](axis="y", labelsize = fontsize-2, labelright="on")
ax[3][:tick_params](axis="x", labelsize = fontsize-2)
ax[3][:axis]("tight")
ax[3][:grid](true)
ax[3][:axhline](y=.0, color = "k", linestyle="dashed")

#fig[:suptitle]("    Bias", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f2b.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#  .o88o.  o8o                                              .oooo.
#  888 `"  `"'                                            .dP""Y88b
# o888oo  oooo   .oooooooo oooo  oooo  oooo d8b  .ooooo.        ]8P'
#  888    `888  888' `88b  `888  `888  `888""8P d88' `88b     <88b.
#  888     888  888   888   888   888   888     888ooo888      `88b.
#  888     888  `88bod8P'   888   888   888     888    .o o.   .88P
# o888o   o888o `8oooooo.   `V88V"V8P' d888b    `Y8bod8P' `8bd88P'
#               d"     YD
#               "Y88888P'
#
#
#
# Explore higher order bias impact on the quadraitc delenser

using PyPlot, JLD
savefigs = true
savepath = "scripts/saved_run_0x3f77869c073f2346/"

rng_true = load(savepath*"out_figure3a.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure3a.jld", "rest_B")
rest_1   = load(savepath*"out_figure3a.jld", "rest_modes")
rest_B_2 = load(savepath*"out_figure3b.jld", "rest_B")
rest_2   = load(savepath*"out_figure3b.jld", "rest_modes")

rest_B_3 = load(savepath*"out_figure3c.jld", "rest_B")
rest_3   = load(savepath*"out_figure3c.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_B_2  = mean((rest_B_2 .- rng_true),2)
bias_B_3  = mean((rest_B_3 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)
bias_L_2  = mean((rest_2 .- rng_true),2)
bias_L_3  = mean((rest_3 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_B_2 = sqrt.(mean((rest_B_2 .- rng_true).^2, 2))
standard_error_B_3 = sqrt.(mean((rest_B_3 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))
standard_error_L_2 = sqrt.(mean((rest_2 .- rng_true).^2, 2))
standard_error_L_3 = sqrt.(mean((rest_3 .- rng_true).^2, 2))

maxylab = maximum(vcat(
    rng_true./standard_error_B_1,
    rng_true./standard_error_B_2,
    rng_true./standard_error_B_3,
    rng_true./standard_error_L_1,
    rng_true./standard_error_L_2,
    rng_true./standard_error_L_3
    ))

minbias, maxbias = vcat(
         bias_B_1./standard_error_B_1,
         bias_B_2./standard_error_B_2,
         bias_B_3./standard_error_B_3,
         bias_L_1./standard_error_L_1,
         bias_L_2./standard_error_L_2,
         bias_L_3./standard_error_L_3,
         ) |> x-> (minimum(x), maximum(x))

#' f2a.pdf : plot r / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[:set_title](L"Quadratic delenser ($\mu K arcmin = 1$) ", fontsize=fontsize-1)
ax[:plot](rng_true, rng_true./standard_error_B_1, "C0", label = "Linear order (null test)")
ax[:plot](rng_true, rng_true./standard_error_B_2, "C3.", label = "All order E from E")
ax[:plot](rng_true, rng_true./standard_error_B_3, "C4", label = "All order E from B")
ax[:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f3a.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f2b.pdf : plot  Bias(r) / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[:set_title](L"Quadratic delenser ($\mu K arcmin = 1$) ", fontsize=fontsize-1)
ax[:plot](rng_true, bias_B_1./rng_true, "C0", label = "Linear order (null test)")
ax[:plot](rng_true, bias_B_2./rng_true, "C3.", label = "All order E from E")
ax[:plot](rng_true, bias_B_3./rng_true, "C4", label = "All order E from B")
ax[:set_ylabel](L"$\frac{Bias(r)}{r}$        ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:axhline](y=.0, color = "k", linestyle="dashed")

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f3b.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#  .o88o.  o8o                                                       .o
#  888 `"  `"'                                                     .d88
# o888oo  oooo   .oooooooo oooo  oooo  oooo d8b  .ooooo.         .d'888
#  888    `888  888' `88b  `888  `888  `888""8P d88' `88b      .d'  888
#  888     888  888   888   888   888   888     888ooo888      88ooo888oo
#  888     888  `88bod8P'   888   888   888     888    .o           888
# o888o   o888o `8oooooo.   `V88V"V8P' d888b    `Y8bod8P'          o888o
#               d"     YD
#               "Y88888P'
#
#
#
#
#
#
#
#
#
#
#

using PyPlot, JLD
savefigs = true
savepath = "scripts/saved_run_0x3f77869c073f2346/"

rng_true = load(savepath*"out_figure4.jld", "rng_true")

# low resolution survey for ϕ est
rest_B_1 = load(savepath*"out_figure4.jld", "rest_B")
rest_1   = load(savepath*"out_figure4.jld", "rest_modes")

# compute bias and standard error
bias_B_1  = mean((rest_B_1 .- rng_true),2)
bias_L_1  = mean((rest_1 .- rng_true),2)

standard_error_B_1 = sqrt.(mean((rest_B_1 .- rng_true).^2, 2))
standard_error_L_1 = sqrt.(mean((rest_1 .- rng_true).^2, 2))

maxylab = maximum(vcat(
    rng_true./standard_error_B_1,
    rng_true./standard_error_L_1,
    ))

minbias, maxbias = vcat(
         bias_B_1./standard_error_B_1,
         bias_L_1./standard_error_L_1,
         ) |> x-> (minimum(x), maximum(x))

#' f4a.pdf : plot r / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[:set_title](L"Nonstationary noise (baseline $\mu K arcmin = 1$)", fontsize=fontsize-1)
ax[:plot](rng_true, rng_true./standard_error_B_1, label = "quadratic delenser")
ax[:plot](rng_true, rng_true./standard_error_L_1, label = "likelihood approx")
ax[:set_ylabel](L"$\frac{r}{\sigma(r)}$       ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)

#fig[:suptitle]("    Detection", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f4a.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f4b.pdf : plot  Bias(r) / σ(r)
#' ================================================

fig, ax = plt[:subplots](1,1, figsize=(4.2,3.7), sharex=true,  sharey=true)
fontsize = 10

#  CMB-S3 for ϕ est
ax[:set_title](L"Nonstationary noise (baseline $\mu K arcmin = 1$)", fontsize=fontsize-1)
ax[:plot](rng_true, bias_B_1./rng_true, label = "quadratic delenser")
ax[:plot](rng_true, bias_L_1./rng_true, label = "likelihood approx")
ax[:set_ylabel](L"$\frac{Bias(r)}{r}$        ", fontsize=fontsize+4,  rotation="horizontal")
ax[:set_xlabel](L"$r$", fontsize=fontsize)
ax[:set_xticks]([.01, .03, .05, .07, .09])
ax[:tick_params](axis="y", labelsize = fontsize-2)
ax[:tick_params](axis="x", labelsize = fontsize-2)
ax[:legend](loc="best", fontsize=fontsize)
ax[:grid](true)
ax[:axhline](y=.0, color = "k", linestyle="dashed")

#fig[:suptitle]("    Bias", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f4b.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end

#' f4c.pdf :
#' ================================================

fig = figure(figsize=(4.2,3.7))
fontsize = 10
σmod       = load(savepath*"out_figure4.jld", "σmod")
μKarcminT  = load(savepath*"out_figure4.jld", "μKarcminT")
#imshow(fftshift(σmod .* √2 .* μKarcminT))
imshow(fftshift(σmod))
colorbar()
title("Non-stationary noise modulation", fontsize=fontsize)
axis("tight")
axis("off")

#fig[:suptitle]("    Bias", fontsize=fontsize+3)
fig[:tight_layout]()

#' Now save
savename  = "f4c.pdf"
if savefigs
     fig[:savefig](joinpath(savepath, savename),bbox_inches="tight")
end


text:
# Coded by Adrien I, 12/27/2020

using Statistics
using Random

mutable struct Self
    features::Array{Any, 2} #change if you have differents features
    targets::Array{Any, 1} #change if you have differents targets, but in general, this size is good
    features_for_forward::Array{Any,1}
    weights::Array{Any,1}
    bias::Array{Any,1}
    eta::Float64
    result::Array
end

function activation(z)
    1 ./ (1 .+ exp.(-z))
end

function derivative_activation(z)
    sigmoid(z) .* (1 .- sigmoid(z))
end

function for_prop(self::Self)
    self.features_for_forward = [self.features]
    for i in 1:length(self.weights)
        z = self.features_for_forward[i] * self.weights[i] .+ self.bias[i] #pre_activation of our neuron
        push!(self.features_for_forward, sigmoid(z)) #activation of the neuron
    end
    self.result = self.features_for_forward[end]
    return self
end

function back_prop(self::Self)
    nabla = (self.targets - self.result) .* sigmoid_derivative(self.result)

    self.weights[end] += self.eta * transpose(self.features_for_forward[end-1]) * nabla
    self.bias[end] .+= self.eta * mean(nabla)

    for i in 1:length(self.weights)-1 #update our gradient
        J = length(self.weights) - i

        cost = nabla * transpose(self.weights[J+1])
        nabla = cost .* sigmoid_derivative(self.features_for_forward[J+1])
        self.weights[J] += self.eta * transpose(self.features_for_forward[J]) * nabla
        self.bias[J] .+= self.eta * mean(nabla)
    end

    return self
end

function train(self::Self, epochs=1000)
    for i in 1:epochs
        for_prop(self) #forward propagation
        back_prop(self) #back prop
    end
end

function predict(self::Self)
    forward(self) #simple forward througout our NN to predict
end

function init_var(input_layer::Int64, hidden_layer, output_layer::Int64, eta)
    sizes = [input_layer hidden_layer output_layer]

    weights = []
    bias = []
    for i in 1:length(sizes)
        push!(weights, randn((sizes[i-1], sizes[i]))) #create weights
        push!(bias, zeros((1, sizes[i]))) #create bias
    end
   return self([], [], [], weights, bias, eta, [])
end

#=

Example of code to generate features and targets, you just have to add it into your code
Base on an example sick/healthy (basic example to test a regression or a NN)

function gen_features(self::Self)
    input = 2
    hidden = [10, 20, 30, 20, 15, 10]
    output = 1
    eta = 0.005

    size_per_class = 200

    sick = randn((size_per_class, 2)) .+ 2         #sick > 0
    healthy = randn((size_per_class, 2)) .- 2      #healthy > 0

    features = vcat(sick, healthy)
    targets = vcat(zeros(size_per_class), ones(size_per_class))

    self.features = features
    self.targets = targets

    init_var(input, hidden, output, eta)

    return self
end

=#

ext: jl


text:
# # Flock model

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../flocking.mp4" type="video/mp4">
# </video>
# ```

# The flock model illustrates how flocking behavior can emerge when each bird follows three simple rules:
#
# * maintain a minimum distance from other birds to avoid collision
# * fly towards the average position of neighbors
# * fly in the average direction of neighbors

# It is also available from the `Models` module as [`Models.flocking`](@ref).

# ## Defining the core structures

# We begin by calling the required packages and defining an agent type representing a bird.

using Agents, LinearAlgebra
using Random # hide

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    separate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

# The fields `id` and `pos` are required for every agent.
# The field `vel` is required for using [`move_agent!`](@ref) in `ContinuousSpace`.
# `speed` defines how far the bird travels in the direction defined by `vel` per `step`.
# `seperation` defines the minimum distance a bird must maintain from its neighbors.
# `visual_distance` refers to the distance a bird can see and defines a radius of neighboring birds.
# The contribution of each rule defined above recieves an importance weight: `cohere_factor`
# is the importance of maintaining the average position of neighbors,
# `match_factor` is the importance of matching the average trajectory of neighboring birds,
# and `separate_factor` is the importance of maining the minimum
# distance from neighboring birds.

# The function `initialize_model` generates birds and returns a model object using default values.
function initialize_model(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    extent = (100, 100),
    spacing = visual_distance / 1.5,
)
    space2d = ContinuousSpace(extent, spacing)
    model = ABM(Bird, space2d, scheduler = Schedulers.randomly)
    for _ in 1:n_birds
        vel = Tuple(rand(model.rng, 2) * 2 .- 1)
        add_agent!(
            model,
            vel,
            speed,
            cohere_factor,
            separation,
            separate_factor,
            match_factor,
            visual_distance,
        )
    end
    return model
end
nothing # hide

# ## Defining the agent_step!
# `agent_step!` is the primary function called for each step and computes velocity
# according to the three rules defined above.
function agent_step!(bird, model)
    ## Obtain the ids of neighbors within the bird's visual distance
    neighbor_ids = nearby_ids(bird, model, bird.visual_distance)
    N = 0
    match = separate = cohere = (0.0, 0.0)
    ## Calculate behaviour properties based on neighbors
    for id in neighbor_ids
        N += 1
        neighbor = model[id].pos
        heading = neighbor .- bird.pos

        ## `cohere` computes the average position of neighboring birds
        cohere = cohere .+ heading
        if edistance(bird.pos, neighbor, model) < bird.separation
            ## `separate` repels the bird away from neighboring birds
            separate = separate .- heading
        end
        ## `match` computes the average trajectory of neighboring birds
        match = match .+ model[id].vel
    end
    N = max(N, 1)
    ## Normalise results based on model input and neighbor count
    cohere = cohere ./ N .* bird.cohere_factor
    separate = separate ./ N .* bird.separate_factor
    match = match ./ N .* bird.match_factor
    ## Compute velocity based on rules defined above
    bird.vel = (bird.vel .+ cohere .+ separate .+ match) ./ 2
    bird.vel = bird.vel ./ norm(bird.vel)
    ## Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

# ## Plotting the flock
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

# The great thing about [`abm_plot`](@ref) is its flexibility. We can incorporate the
# direction of the birds when plotting them, by making the "marker" function `am`
# create a `Polygon`: a triangle with same orientation as the bird's velocity.
# It is as simple as defining the following function:

const bird_polygon = Polygon(Point2f0[(-0.5, -0.5), (1, 0), (-0.5, 0.5)])
function bird_marker(b::Bird)
    φ = atan(b.vel[2], b.vel[1]) #+ π/2 + π
    scale(rotate2D(bird_polygon, φ), 2)
end

# Where we have used the utility functions `scale` and `rotate2D` to act on a
# predefined polygon. We now give `bird_marker` to `abm_plot`, and notice how
# the `as` keyword is meaningless when using polygons as markers.

model = initialize_model()
figure, = abm_plot(model; am = bird_marker)
figure

# And let's also do a nice little video for it:
abm_video(
    "flocking.mp4", model, agent_step!;
    am = bird_marker,
    framerate = 20, frames = 100,
    title = "Flocking"
)

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../flocking.mp4" type="video/mp4">
# </video>
# ```


text:
# Timings for the skeleton BATSRUS for porting to GPU with OpenACC.
# 4 subplots, each of block size 8^3, 16^3, 32^3, and 64^3.
# Each subplot contains a bar plot, with nBlock as x axis, timing as y axis,
# and different columns represent nWork.
# Column have two colors, blue represent CPU, and orange represent GPU.
#
# Hongyang Zhou, hyzhou@umich.edu 04/22/2020

using PyPlot, CSV, DataFrames

#
df = DataFrame(CSV.File("results/timing_gpu.txt"))
t1 = df[df.nCell.==8,:]
t2 = df[df.nCell.==16,:]
t3 = df[df.nCell.==32,:]
t4 = df[df.nCell.==64,:]

t1_cpu = t1[t1.type.=="CPU",:]

# set width of bar
const barWidth = 0.3

# Set position of bar on X axis
r1 = collect(0:size(t1,1)-1)
r2 = [x + barWidth for x in r1]
r3 = [x + barWidth for x in r2]

fig, ax = subplots(2,2,figsize=(14,6))
#axis("tight")

# Make the plot
#a1 = ax[1,1].bar(r1, t1[:,:time], width=barWidth, edgecolor="white", label="1 thread * 256 MPI")
#a2 = ax[1,2].bar(r2, t2[:,:time], width=barWidth, edgecolor="white", label="8 threads * 32 MPI")
#a3 = ax[2,1].bar(r3, t3[:,:time], width=barWidth, edgecolor="white", label="32 threads * 8 MPI")
#a3 = ax[2,2].bar(r3, t4[:,:time], width=barWidth, edgecolor="white", label="32 threads * 8 MPI")
#=
# Add xticks on the middle of the group bars
xlabel("module", fontweight=:bold, fontsize=18)
xticks([r + barWidth for r in 0:length(t1)-1],
   ["Advance", "FaceValue", "FaceFlux", "Update", "UpdateCheck", "MessagePass"],
   fontsize=18)
plt.minorticks_on()
ax.tick_params(axis=:x, which=:minor, bottom=false)
ax.yaxis.grid(true, linestyle="--")
#ax.yaxis.grid(true, which="minor")
ylabel("time [s]", fontsize=18)
ax.set_title("Timings per time step for 256 cores, explicit 2nd Order Linde + mc3 limiter",fontsize=18)

# Create legend & Show graphic
legend(fontsize=18)

function autolabel(rects, t_percent)
   """Attach a text label above each bar in *rects*, displaying its height."""
   for (i, rect) in enumerate(rects)
      i == 1 && continue
      height = rect.get_height()
      percent = t_percent[i]
      ax.annotate(string(percent),
         xy=(rect.get_x() + rect.get_width() / 2, height),
         xytext=(0, 3),  # 3 points vertical offset
         textcoords="offset points",
         ha=:center, va=:bottom)
   end
end

autolabel(a1, t1_percent)
autolabel(a2, t2_percent)
autolabel(a3, t3_percent)
=#


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 2, Zadanie 3.
#

# imports
using LinearAlgebra
include("hilb.jl")
include("matcond.jl")
using Printf

# solve Ax = b using two methods
# in this particular case x is just a column of ones, but we’re trying
# to recalculate it to see how bad is the given matrix in terms of its
# condition number `cond(A)`
function solve(A)
    # prepare
    x = ones(size(A)[1])
    b = A * x
    # calc using two methods
    gauss = A \ b
    inverted = inv(A) * b

    return gauss, inverted
end

# this function calculates the relative error between the original
# and the recalculated value
function calc_relative_error(x, x̃)
    return (norm(x - x̃)/norm(x))
end

# scientific notation formatting function
function scn(n)
    return @sprintf("%.0e", n)
    # return string(n)
end

# generates and print out the results
function printout(A, n, c, latex_table=false)

    gauss, inverted = solve(A)

    if !latex_table
        print(
            "n = ", n, "\n",
            (c != 0 ? string("c = ", scn(c), "\n") : ""),
            "  cond(A)                   = ", cond(A), "\n",
            "  rank(A)                   = ", rank(A), "\n",
            "  relative error (gauss)    = ", calc_relative_error(ones(n), gauss), "\n",
            "  relative error (inverted) = ", calc_relative_error(ones(n), inverted), "\n\n"
        )
    else
        println(
            n, " & ",
            (c != 0 ? string(scn(c), " & ") : ""),
            cond(A), " & ",
            rank(A), " & ",
            calc_relative_error(ones(n), gauss), " & ",
            calc_relative_error(ones(n), inverted), "\\\\"
        )
    end
end

# 1. Hilbert’s matrix Hₙ for n > 1
for n in 2:24

    A = hilb(n)
    printout(A, n, 0, false)

end

# 2. Random matrix Rₙ for n = 5,10,20 and c = 1, 10, 10³, 10⁷, 10^12, 10^16
for n in [5,10,20]

    for c in [1, 10, 10^3, 10^7, 10^12, 10^16]

        A = matcond(n, Float64(c))
        printout(A, n, c, false)

    end

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
using PkgBenchmark
using CategoricalArrays

@benchgroup "isequal(A, v::String)" begin
    function sumequals(A::AbstractArray, v::Any)
        n = 0
        @inbounds for x in A
            n += isequal(x, v)
        end
        n
    end

    ca = CategoricalArray(repeat(string.('A':'J'), outer=1000))
    nca = NullableCategoricalArray(repeat([Nullable(); string.('A':'J')], outer=1000))
    @bench "CategoricalArray" sumequals(ca, "D")
    @bench "NullableCategoricalArray" sumequals(nca, Nullable("D"))
end

@benchgroup "isequal(A, v::CategoricalValue)" begin
    function sumequals(A::AbstractArray, v::CategoricalValue)
        n = 0
        @inbounds for x in A
            n += isequal(x, v)
        end
        n
    end

    ca = CategoricalArray(repeat(string.('A':'J'), outer=1000))
    nca = NullableCategoricalArray(repeat([Nullable(); string.('A':'J')], outer=1000))
    @bench "CategoricalArray" sumequals(ca, ca[1])
    @bench "NullableCategoricalArray" sumequals(nca, nca[1])
end


text:
###############################################################################
#
# update the passwords for data portals
#
###############################################################################
"""
    update_MODIS_password!()

Update global user name and password for LP DAAC, if eithe of them is empty
"""
function update_MODIS_password!()
    # input MODIS_USER_ID and MODIS_USER_PWD
    global MODIS_USER_ID, MODIS_USER_PWD;
    if MODIS_USER_ID == "" || MODIS_USER_PWD == ""
        @warn "Do not share your user name and password with others.";
        @info "Please indicate your user name for LP DAAC data portal:";
        MODIS_USER_ID  = readline();
        @info "Please indicate your password for LP DAAC Data portal:";
        MODIS_USER_PWD = read(Base.getpass("Password (invisble)"), String);
    end

    return nothing
end

"""
    update_CDSAPI_password!()

Update global user name and password for CDSAPI, note that the `uid:api-key`
    should not contain `{}`.
"""
function update_CDSAPI_password!()
    global CDSAPI_CLIENT, CDSAPI_KEY, CDSAPI_PORTAL;
    if CDSAPI_KEY=="" || CDSAPI_PORTAL==""
        # Install cdsapi if not existing already
        try
            pyimport("cdsapi");
        catch e
            install_cdsapi!();
        end

        # update password and client
        if isfile("$(homedir())/.cdsapirc")
            @info "Load cdsapi information from $(homedir())/.cdsapirc";
            cdsapi        = pyimport("cdsapi");
            CDSAPI_CLIENT = cdsapi.Client();
            CDSAPI_KEY    = CDSAPI_CLIENT.key;
            CDSAPI_PORTAL = CDSAPI_CLIENT.url;
        else
            @warn "Do not share your user name and password with others.";
            @info "Please indicate the cdsapi data portal (url):";
            CDSAPI_PORTAL = readline();
            @info "Please indicate the cdsapi key combo (uid:api-key):";
            CDSAPI_KEY    = read(Base.getpass("Password (invisble)"), String);
            cdsapi        = pyimport("cdsapi");
            CDSAPI_CLIENT = cdsapi.Client(url=CDSAPI_PORTAL, key=CDSAPI_KEY);
        end
    end

    return nothing
end


text:
using CSV, StateSpaceModels, Plots, Statistics, Dates

# Load the AirPassengers dataset
AP = CSV.read("AirPassengers.csv")

# Take the log of the series
logAP = log.(Vector{Float64}(AP[:Passengers]))

# Plot the data
p1 = plot(AP[:Date], logAP, label = "Log-airline passengers", legend = :topleft, color = :black)

# Specify the state-space model
model = structural(logAP, 12)

# Estimate the state-space model
ss = statespace(model)

# Analyze its decomposition in trend and seasonal
p2 = plot(AP[:Date], [ss.smoother.alpha[:, 1] ss.smoother.alpha[:, 3]], layout = (2, 1),
            label = ["Trend component" "Seasonal component"], legend = :topleft)

# Forecast 24 months ahead
N = 24
pred, dist = forecast(ss, N)

# Define forecasting dates
firstdate = AP[:Date][end] + Month(1)
newdates = collect(firstdate:Month(1):firstdate + Month(N - 1))

p3 = plot!(p1, newdates, pred, label = "Forecast")


text:
function mulαβ!(C::Matrix{T}, A::Matrix{T}, adjB::Adjoint{T,<:Matrix{T}},
        α=true, β=false) where {T<:BlasFloat}
    BLAS.gemm!('N', 'C', T(α), A, adjB.parent, T(β), C)
end

function mulαβ!(C::Matrix{T}, A::Matrix{T}, adjB::Diagonal{T,S},
        α=true, β=false) where {T<:BlasFloat,S}
    # adapted from LinearAlgebra/src/diagonal.jl in 1.3
    C .= (A .* permutedims(adjB.diag)) .* α .+ C .* β
end

function mulαβ!(C::SparseMatrixCSC{T}, A::SparseMatrixCSC{T}, adjB::Diagonal{T,S},
        α=true, β=false) where {T<:BlasFloat,S}
    # adapted from LinearAlgebra/src/diagonal.jl in 1.3
    C .= (A .* permutedims(adjB.diag)) .* α .+ C .* β
end

function mulαβ!(C::Matrix{T}, A::SparseMatrixCSC{T}, adjB::Adjoint{T,<:SparseMatrixCSC{T}},
        α=true, β=false) where T <: Number
    B = adjB.parent
    B.m == size(C, 2) && A.m == size(C, 1) && A.n == B.n || throw(DimensionMismatch(""))
    anz = nonzeros(A)
    arv = rowvals(A)
    bnz = nonzeros(B)
    brv = rowvals(B)
    isone(β) || rmul!(C, β)
    @inbounds for j = 1:A.n
        for ib in nzrange(B, j)
            αbnz = α * bnz[ib]
            jj = brv[ib]
            for ia in nzrange(A, j)
                C[arv[ia], jj] += anz[ia] * αbnz
            end
        end
    end
    C
end

mulαβ!(C::Matrix{T}, A::BlockedSparse{T}, adjB::Adjoint{T,<:BlockedSparse{T}}, α=true,
       β=false) where {T} = mulαβ!(C, A.cscmat, adjB.parent.cscmat', α, β)

function mulαβ!(C::Matrix{T}, A::SparseMatrixCSC{T}, adjB::Adjoint{T,Matrix{T}},
        α=true, β=false) where {T}
    B = adjB.parent
    A.n == size(B, 2) || throw(DimensionMismatch())
    A.m == size(C, 1) || throw(DimensionMismatch())
    size(B, 1) == size(C, 2) || throw(DimensionMismatch())
    nzv = A.nzval
    rv = A.rowval
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        @inbounds for col = 1:A.n
            αxj = α*B[k,col]
            for j = nzrange(A, col)
                C[rv[j], k] += nzv[j]*αxj
            end
        end
    end
    C
end

mulαβ!(C::Matrix{T}, A::BlockedSparse{T}, adjB::Adjoint{T,<:Matrix{T}}, α=true, β=false) where {T} =
    mulαβ!(C, A.cscmat, adjB, α, β)

function mulαβ!(C::SparseMatrixCSC{T}, A::SparseMatrixCSC{T}, adjB::Adjoint{T,<:SparseMatrixCSC{T}},
        α=true, β=false) where {T}
    B = adjB.parent
    C.m == A.m && C.n == B.m && A.n == B.n || throw(DimensionMismatch(""))
    Anz = nonzeros(A)
    Bnz = nonzeros(B)
    Cnz = nonzeros(C)
    isone(β) || rmul!(Cnz, β)
    Arv = rowvals(A)
    Brv = rowvals(B)
    Crv = rowvals(C)
    for j in 1:A.n
        for K in nzrange(B, j)
            k = Brv[K]
            alphabjk = α * Bnz[K]
            colkfirstr = Int(C.colptr[k])
            colklastr = Int(C.colptr[k + 1] - 1)
            for I in nzrange(A, j)
                i = Arv[I]
                searchk = searchsortedfirst(Crv, i, colkfirstr, colklastr, Base.Order.Forward)
                if searchk <= colklastr && Crv[searchk] == i
                    Cnz[searchk] += alphabjk * Anz[I]
                else
                    throw(ArgumentError("C does not have the nonzero pattern of A*B'"))
                end
            end
        end
    end
    C
end

mulαβ!(C::BlockedSparse{T}, A::BlockedSparse{T}, adjB::Adjoint{T,<:BlockedSparse{T}}, α=true, β=false) where {T} =
    mulαβ!(C.cscmat, A.cscmat, adjB.parent.cscmat', α, β)

function mulαβ!(C::StridedVecOrMat{T}, A::StridedVecOrMat{T}, adjB::Adjoint{T,<:SparseMatrixCSC{T}},
        α=true, β=false) where T
    B = adjB.parent
    m, n = size(A)
    p, q = size(B)
    r, s = size(C)
    r == m && s == p && n == q || throw(DimensionMismatch(""))
    isone(β) || rmul!(C, β)
    nz = nonzeros(B)
    rv = rowvals(B)
    @inbounds for j in 1:q, k in nzrange(B, j)
        rvk = rv[k]
        anzk = α * nz[k]
        for jj in 1:r
            C[jj, rvk] += A[jj, j] * anzk
        end
    end
    C
end

mulαβ!(C::StridedVecOrMat{T}, A::StridedVecOrMat{T}, adjB::Adjoint{T,<:BlockedSparse{T}},
    α=true, β=false) where {T} = mulαβ!(C, A, adjB.parent.cscmat', α, β)

mulαβ!(C::StridedVector{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, B::StridedVector{T},
    α=true, β=false) where {T<:BlasFloat} = BLAS.gemv!('C', T(α), adjA.parent, B, T(β), C)

mulαβ!(C::StridedVector{T}, adjA::Adjoint{T,<:SparseMatrixCSC{T}}, B::StridedVector{T},
    α=true, β=false) where {T} = mul!(C, adjA, B, T(α), T(β))

mulαβ!(C::StridedVector{T}, adjA::Adjoint{T,<:BlockedSparse{T}}, B::StridedVector{T},
    α=true, β=false) where {T} = mulαβ!(C, adjA.parent.cscmat', B, α, β)

function LinearAlgebra.ldiv!(adjA::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
        B::StridedVector{T}) where {T}
    A = adjA.parent
    length(B) == size(A, 2) || throw(DimensionMismatch(""))
    m, n, k = size(A.data.data)
    fv = A.data.facevec
    bb = reshape(B, (n, k))
    for j in 1:k
        ldiv!(adjoint(LowerTriangular(fv[j])), view(bb, :, j))
    end
    B
end

function LinearAlgebra.rdiv!(A::Matrix{T},
        adjB::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}}) where T
    Bd = adjB.parent.data
    m, n, k = size(Bd.data)
    size(A, 2) == size(Bd, 1) && m == n || throw(DimensionMismatch(""))
    inds = 1:m
    for (i, f) in enumerate(Bd.facevec)
        BLAS.trsm!('R', 'L', 'T', 'N', one(T), f, view(A, :, inds .+ m * (i-1)))
    end
    A
end

function LinearAlgebra.rdiv!(A::BlockedSparse{T,S,P},
        B::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}}) where {T,S,P}
    Bpd = B.parent.data
    j,k,l = size(Bpd.data)
    cbpt = A.colblkptr
    nzv = A.cscmat.nzval
    P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
    for (j,f) in enumerate(Bpd.facevec)
        rdiv!(reshape(view(nzv, cbpt[j]:(cbpt[j + 1] - 1)), :, P), adjoint(LowerTriangular(f)))
    end
    A
end


text:
using Documenter
using PotentialCalculation, PotentialCalculation.psi4

makedocs(sitename="PotentialCalculation.jl",
         pages=["Home" => "index.md",
                "Install" => "install.md",
                "Usage" => "use.md",
                "References" => "reference.md"]

)

deploydocs(
    repo = "github.com/MatrixLabTools/PotentialCalculation.jl.git",
)


text:
@testset "faces" begin

    @testset "constructors" begin
        f1 = Face(1,2,3)
        @test Face{3, ZeroIndex{Int}}(f1) == Face{3, Int}(1,2,3)
        # round trip
        i = 2
        tmp = Face{3, OffsetInteger{3, Int32}}(Face{3}(ZeroIndex{Int}(i)))
        @test Int(tmp[1]) == i
    end

    @testset "getindex" begin
        a = [1,2,3,4]
        @test a[Face{3, Int}(1,2,3)] == (1,2,3)
        @test a[Face{3, ZeroIndex{Int}}(1,2,3)] == (1,2,3)
    end

    @testset "setindex" begin
        a = [1,2,3,4]
        a[Face(1,2,3)] = [7,6,5]
        @test a == [7,6,5,4]
    end

    @testset "views" begin
        x = GLNormalMesh(AABB(Vec3f0(0), Vec3f0(1)))
        v = view(vertices(x), faces(x))
        vs = vertices(x)
        vf = [vs[face] for face in faces(x)]
        @test v == vf
    end

    @testset "negative offests and unsigned integers" begin
        oi = OffsetInteger{-1, UInt64}(UInt64(1))
        show(IOBuffer(), oi)

        oi = OffsetInteger{1, UInt64}(UInt64(1))
        show(IOBuffer(), oi)
    end
end


text:
@testset "MCDM functions" begin
    @testset "Grey Relational Analysis" begin

        tol = 0.0001

        df = DataFrame(
            :K1 => [105000.0, 120000, 150000, 115000, 135000],
            :K2 => [105.0, 110, 120, 105, 115],
            :K3 => [10.0, 15, 12, 20, 15],
            :K4 => [4.0, 4, 3, 4, 5],
            :K5 => [300.0, 500, 550, 600, 400],
            :K6 => [10.0, 8, 12, 9, 9],
        )
        functionlist = [minimum, maximum, minimum, maximum, maximum, minimum]

        w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40]

        result = grey(df, w, functionlist)

        @test isa(result, GreyResult)

        @test isapprox(
            result.scores,
            [
                0.525,
                0.7007142857142857,
                0.5464285714285715,
                0.5762820512820512,
                0.650952380952381,
            ],
            atol = tol,
        )

        @test result.bestIndex == 2

        setting = MCDMSetting(df, w, functionlist)
        result2 = grey(setting)
        @test result2 isa GreyResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, GreyMethod())
        @test result3 isa GreyResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "TOPSIS" begin
        tol = 0.00001
        df = DataFrame()
        df[:, :x] = Float64[9, 8, 7]
        df[:, :y] = Float64[7, 7, 8]
        df[:, :z] = Float64[6, 9, 6]
        df[:, :q] = Float64[7, 6, 6]
        w = Float64[4, 2, 6, 8]

        fns = makeminmax([maximum, maximum, maximum, maximum])

        result = topsis(df, w, fns)

        @test isa(result, TopsisResult)
        @test result.bestIndex == 2
        @test isapprox(result.scores, [0.3876870, 0.6503238, 0.0834767], atol = tol)

        setting = MCDMSetting(df, w, fns)
        result2 = topsis(setting)
        @test isa(result2, TopsisResult)
        @test result2.bestIndex == result.bestIndex
        @test result.scores == result2.scores

        result3 = mcdm(setting, TopsisMethod())
        @test isa(result3, TopsisResult)
        @test result3.bestIndex == result.bestIndex
        @test result3.scores == result.scores
    end

    @testset "VIKOR" begin
        tol = 0.00001
        w = [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
        Amat = [
            100 92 10 2 80 70 95 80
            80 70 8 4 100 80 80 90
            90 85 5 0 75 95 70 70
            70 88 20 18 60 90 95 85
        ]
        dmat = makeDecisionMatrix(Amat)
        fns = makeminmax([
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
        ])

        result = vikor(dmat, w, fns)

        @test isa(result, VikorResult)
        @test result.bestIndex == 4

        @test isapprox(result.scores[1], 0.1975012087551764, atol = tol)
        @test isapprox(result.scores[2], 0.2194064473270817, atol = tol)
        @test isapprox(result.scores[3], 0.3507643203516215, atol = tol)
        @test isapprox(result.scores[4], -0.16727341435277993, atol = tol)

        setting = MCDMSetting(dmat, w, fns)
        result2 = vikor(setting)
        @test result.scores == result2.scores
        @test result.bestIndex == result2.bestIndex

        result3 = mcdm(setting, VikorMethod())
        @test result.scores == result3.scores
        @test result.bestIndex == result3.bestIndex
    end

    @testset "ELECTRE" begin
        tol = 0.00001
        w = [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
        Amat = [
            100 92 10 2 80 70 95 80
            80 70 8 4 100 80 80 90
            90 85 5 0 75 95 70 70
            70 88 20 18 60 90 95 85
        ]
        dmat = makeDecisionMatrix(Amat)
        fns = makeminmax([maximum for i = 1:8])
        result = electre(dmat, w, fns)

        @test isa(result, ElectreResult)
        @test isa(result.bestIndex, Tuple)
        @test result.bestIndex[1] == 4

        @test isapprox(
            result.C,
            [0.36936937, 0.01501502, -2.47347347, 2.08908909],
            atol = tol,
        )
        @test isapprox(result.D, [0.1914244, -0.1903929, 2.8843076, -2.8853391], atol = tol)

        setting = MCDMSetting(dmat, w, fns)
        result2 = electre(setting)
        @test result isa ElectreResult
        @test result.bestIndex == result2.bestIndex

        result3 = mcdm(setting, ElectreMethod())
        @test result3 isa ElectreResult
        @test result3.bestIndex == result.bestIndex
    end

    @testset "MOORA" begin
        tol = 0.00001
        w = [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
        Amat = [
            100 92 10 2 80 70 95 80
            80 70 8 4 100 80 80 90
            90 85 5 0 75 95 70 70
            70 88 20 18 60 90 95 85
        ]
        dmat = makeDecisionMatrix(Amat)
        fns = makeminmax([
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
            maximum,
        ])
        result = moora(dmat, w, fns)

        @test isa(result, MooraResult)
        @test isa(result.bestIndex, Int64)
        @test result.bestIndex == 4

        @test isapprox(
            result.scores,
            [0.33159387, 0.29014464, 0.37304311, 0.01926526],
            atol = tol,
        )

        setting = MCDMSetting(dmat, w, fns)
        result2 = moora(setting)
        @test result2 isa MooraResult
        @test result2.bestIndex == result.bestIndex
        @test result2.scores == result.scores

        result3 = mcdm(setting, MooraMethod())
        @test result3 isa MooraResult
        @test result3.bestIndex == result.bestIndex
        @test result3.scores == result.scores
    end

    @testset "DEMATEL" begin
        tol = 0.00001
        K = [
            0 3 0 2 0 0 0 0 3 0
            3 0 0 0 0 0 0 0 0 2
            4 1 0 2 1 3 1 2 3 2
            4 1 4 0 1 2 0 1 0 0
            3 2 3 1 0 3 0 2 0 0
            4 1 4 4 0 0 0 1 1 3
            3 0 0 0 0 2 0 0 0 0
            3 0 4 3 2 3 1 0 0 0
            4 3 2 0 0 1 0 0 0 2
            2 1 0 0 0 0 0 0 3 0
        ]

        dmat = makeDecisionMatrix(K)

        result = dematel(dmat)

        @test isapprox(result.threshold, 0.062945, atol = tol)

        @test isapprox(
            result.c,
            [
                0.3991458,
                0.2261648,
                1.0204318,
                0.7538625,
                0.8096760,
                0.9780926,
                0.2717874,
                0.9455390,
                0.5960514,
                0.2937537,
            ],
            atol = tol,
        )

        @test isapprox(
            result.r,
            [
                1.5527024,
                0.7251791,
                0.8551461,
                0.6895615,
                0.2059141,
                0.6790404,
                0.1057168,
                0.3163574,
                0.6484014,
                0.5164858,
            ],
            atol = tol,
        )

        @test isapprox(
            result.influenceMatrix,
            [
                0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0
                1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
                1.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0
                1.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
                1.0 1.0 1.0 1.0 0.0 1.0 0.0 1.0 0.0 0.0
                1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0
                1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
                1.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
                1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
                1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
            ],
            atol = tol,
        )

        @test isapprox(
            result.weights,
            [
                0.1686568559124561,
                0.07991375718719543,
                0.14006200243438863,
                0.10748052790517183,
                0.08789022388276985,
                0.12526272598854982,
                0.03067915023486491,
                0.10489168834828348,
                0.092654758940811,
                0.06250830916550884,
            ],
            atol = tol,
        )

    end

    @testset "AHP - Consistency" begin

        tol = 0.00001

        K = [
            1 7 (1/5) (1/8) (1/2) (1/3) (1/5) 1
            (1/7) 1 (1/8) (1/9) (1/4) (1/5) (1/9) (1/8)
            5 8 1 (1/3) 4 2 1 1
            8 9 3 1 7 5 3 3
            2 4 (1/4) (1/7) 1 (1/2) (1/5) (1/5)
            3 5 (1/2) (1/5) 2 1 (1/3) (1/3)
            5 9 1 (1/3) 5 3 1 1
            1 8 1 (1/3) 5 3 1 1
        ]

        dmat = makeDecisionMatrix(K)
        result::AHPConsistencyResult = ahp_consistency(dmat)

        @test isa(result, AHPConsistencyResult)

        @test result.isConsistent

        @test isapprox(result.CR, 0.07359957154133831, atol = tol)

        @test isapprox(result.CI, 0.10377539587328702, atol = tol)

        @test isapprox(result.lambda_max, 8.72642777111301, atol = tol)

        @test isapprox(
            result.pc,
            [
                8.40982182534543
                8.227923797501001
                8.95201017080425
                8.848075127995905
                8.860427038870013
                8.941497805932498
                8.946071053879708
                8.62559534857526
            ],
            atol = tol,
        )
    end

    @testset "AHP" begin

        tol = 0.00001

        K = [
            1 7 1/5 1/8 1/2 1/3 1/5 1
            1/7 1 1/8 1/9 1/4 1/5 1/9 1/8
            5 8 1 1/3 4 2 1 1
            8 9 3 1 7 5 3 3
            2 4 1/4 1/7 1 1/2 1/5 1/5
            3 5 1/2 1/5 2 1 1/3 1/3
            5 9 1 1/3 5 3 1 1
            1 8 1 1/3 5 3 1 1
        ]

        A1 = [
            1 3 1/5 2
            1/3 1 1/7 1/3
            5 7 1 4
            1/2 3 1/4 1
        ]
        A2 = [
            1 1/2 4 5
            2 1 6 7
            1/4 1/6 1 3
            1/5 1/7 1/3 1
        ]
        A3 = [
            1 1/2 1/6 3
            2 1 1/4 5
            6 4 1 9
            1/3 1/5 1/9 1
        ]
        A4 = [
            1 7 1/4 2
            1/7 1 1/9 1/5
            4 9 1 5
            1/2 5 1/5 1
        ]
        A5 = [
            1 6 2 3
            1/6 1 1/4 1/3
            1/2 4 1 2
            1/3 3 1/2 1
        ]
        A6 = [
            1 1/4 1/2 1/7
            4 1 2 1/3
            2 1/2 1 1/5
            7 3 5 1
        ]
        A7 = [
            1 3 7 1
            1/3 1 4 1/3
            1/7 1/4 1 1/7
            1 3 7 1
        ]
        A8 = [
            1 2 5 8
            1/2 1 3 6
            1/5 1/3 1 3
            1/8 1/6 1/3 1
        ]

        km = makeDecisionMatrix(K)
        as = map(makeDecisionMatrix, [A1, A2, A3, A4, A5, A6, A7, A8])

        result = ahp(as, km)

        @test isa(result, AHPResult)

        @test result.bestIndex == 3

        @test isapprox(
            result.scores,
            [0.2801050, 0.1482273, 0.3813036, 0.1903641],
            atol = tol,
        )
    end

    @testset "NDS" begin
        cases = [
            1.0 2.0 3.0
            2.0 1.0 3.0
            1.0 3.0 2.0
            4.0 5.0 6.0
        ]

        nd = makeDecisionMatrix(cases)

        result = nds(nd)

        @test isa(result, NDSResult)

        @test result.bestIndex == 4

        @test result.ranks == [0, 0, 0, 3]

    end

    @testset "SAW" begin

        @testset "Example 1: 4 criteria × 4 alternatives" begin
            tol = 0.0001
            df = DataFrame(
                :c1 => [25.0, 21, 19, 22],
                :c2 => [65.0, 78, 53, 25],
                :c3 => [7.0, 6, 5, 2],
                :c4 => [20.0, 24, 33, 31],
            )
            weights = [0.25, 0.25, 0.25, 0.25]
            fns = [maximum, maximum, minimum, maximum]
            result = saw(df, weights, fns)

            @test result isa SawResult

            @test isapprox(
                result.scores,
                [0.681277, 0.725151, 0.709871, 0.784976],
                atol = tol,
            )
        end

        @testset "Example 2: 7 criteria × 5 alternatives " begin
            tol = 0.0001
            decmat = [
                4.0 7 3 2 2 2 2
                4.0 4 6 4 4 3 7
                7.0 6 4 2 5 5 3
                3.0 2 5 3 3 2 5
                4.0 2 2 5 5 3 6
            ]

            df = makeDecisionMatrix(decmat)
            weights = [0.283, 0.162, 0.162, 0.07, 0.085, 0.162, 0.076]
            fns = makeminmax([maximum for i = 1:7])
            result = saw(df, weights, fns)

            @test result isa SawResult
            @test isapprox(
                result.scores,
                [0.553228, 0.713485, 0.837428, 0.514657, 0.579342],
                atol = tol,
            )
            @test result.bestIndex == 3
            @test result.ranking == [3, 2, 5, 1, 4]

            setting = MCDMSetting(df, weights, fns)
            result2 = saw(setting)
            @test result2 isa SawResult
            @test result2.scores == result.scores
            @test result2.bestIndex == result.bestIndex

            result3 = mcdm(setting, SawMethod())
            @test result3 isa SawResult
            @test result3.scores == result.scores
            @test result3.bestIndex == result.bestIndex
        end
    end

    @testset "ARAS Additive Ratio Assessment" begin

        tol = 0.0001

        df = DataFrame(
            :K1 => [105000.0, 120000, 150000, 115000, 135000],
            :K2 => [105.0, 110, 120, 105, 115],
            :K3 => [10.0, 15, 12, 20, 15],
            :K4 => [4.0, 4, 3, 4, 5],
            :K5 => [300.0, 500, 550, 600, 400],
            :K6 => [10.0, 8, 12, 9, 9],
        )
        functionlist = [minimum, maximum, minimum, maximum, maximum, minimum]

        w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40]

        result = aras(df, w, functionlist)
        @test isa(result, ARASResult)
        @test isapprox(
            result.scores,
            [0.81424068, 0.89288620, 0.76415790, 0.84225462, 0.86540635],
            atol = tol,
        )
        @test result.bestIndex == 2

        setting = MCDMSetting(df, w, functionlist)
        result2 = aras(setting)
        @test result2 isa ARASResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, ArasMethod())
        @test result3 isa ARASResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "WPM" begin
        tol = 0.0001
        decmat = [
            3 12.5 2 120 14 3
            5 15 3 110 38 4
            3 13 2 120 19 3
            4 14 2 100 31 4
            3 15 1.5 125 40 4
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201]

        fns = [maximum, minimum, minimum, maximum, minimum, maximum]

        result = wpm(df, weights, fns)
        @test result isa WPMResult
        @test isapprox(
            result.scores,
            [
                0.7975224331331,
                0.7532541470585,
                0.7647463553356,
                0.7873956894791,
                0.7674278741782,
            ],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = wpm(setting)
        @test result2 isa WPMResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, WPMMethod())
        @test result3 isa MCDMResult
        @test result3 isa WPMResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "WASPAS" begin
        tol = 0.0001
        decmat = [
            3 12.5 2 120 14 3
            5 15 3 110 38 4
            3 13 2 120 19 3
            4 14 2 100 31 4
            3 15 1.5 125 40 4
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201]

        fns = [maximum, minimum, minimum, maximum, minimum, maximum]

        result = waspas(df, weights, fns)
        @test result isa WASPASResult
        @test isapprox(
            result.scores,
            [0.805021, 0.775060, 0.770181, 0.796424, 0.788239],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = waspas(setting)
        @test result2 isa WASPASResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, WaspasMethod())
        @test result3 isa WASPASResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "EDAS" begin
        tol = 0.0001
        decmat = [
            5000 5 5300 450
            4500 5 5000 400
            4500 4 4700 400
            4000 4 4200 400
            5000 4 7100 500
            5000 5 5400 450
            5500 5 6200 500
            5000 4 5800 450
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.25, 0.25, 0.25, 0.25]

        fns = [maximum, maximum, minimum, minimum]

        result = edas(df, weights, fns)
        @test result isa EDASResult
        @test isapprox(
            result.scores,
            [
                0.759594,
                0.886016,
                0.697472,
                0.739658,
                0.059083,
                0.731833,
                0.641691,
                0.385194,
            ],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = edas(setting)
        @test result2 isa EDASResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, EdasMethod())
        @test result3 isa EDASResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "MARCOS" begin
        tol = 0.0001
        decmat = [
            8.675 8.433 8.000 7.800 8.025 8.043
            8.825 8.600 7.420 7.463 7.825 8.229
            8.325 7.600 8.040 7.700 7.925 7.600
            8.525 8.667 7.180 7.375 7.750 8.071
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.19019, 0.15915, 0.19819, 0.19019, 0.15115, 0.11111]

        fns = [maximum, maximum, maximum, maximum, maximum, maximum]

        Fns = convert(Array{Function,1}, fns)

        result = marcos(df, weights, Fns)
        @test result isa MARCOSResult
        @test isapprox(
            result.scores,
            [0.684865943528, 0.672767106696, 0.662596906139, 0.661103207660],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, Fns)
        result2 = marcos(setting)
        @test result2 isa MARCOSResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, MarcosMethod())
        @test result3 isa MARCOSResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "MABAC" begin

        tol = 0.0001
        decmat = [
            2 1 4 7 6 6 7 3000
            4 1 5 6 7 7 6 3500
            3 2 6 6 5 6 8 4000
            5 1 5 7 6 7 7 3000
            4 2 5 6 7 7 6 3000
            3 2 6 6 6 6 6 3500
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.293, 0.427, 0.067, 0.027, 0.053, 0.027, 0.053, 0.053]

        fns = [maximum, maximum, maximum, maximum, maximum, maximum, maximum, minimum]

        result = mabac(df, weights, fns)
        @test result isa MABACResult
        @test isapprox(
            result.scores,
            [-0.31132, -0.10898, 0.20035, 0.04218, 0.34452, 0.20035],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = mabac(setting)
        @test result2 isa MABACResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, MabacMethod())
        @test result3 isa MABACResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "MAIRCA" begin

        tol = 0.0001
        decmat = [
            6.952 8.000 6.649 7.268 8.000 7.652 6.316
            7.319 7.319 6.604 7.319 8.000 7.652 5.313
            7.000 7.319 7.652 6.952 7.652 6.952 4.642
            7.319 6.952 6.649 7.319 7.652 6.649 5.000
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.172, 0.165, 0.159, 0.129, 0.112, 0.122, 0.140]

        fns = [maximum, maximum, maximum, maximum, maximum, maximum, minimum]

        result = mairca(df, weights, fns)
        @test result isa MAIRCAResult
        @test isapprox(
            result.scores,
            [0.1206454, 0.0806646, 0.1458627, 0.1454237],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = mairca(setting)
        @test result2 isa MAIRCAResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, MaircaMethod())
        @test result3 isa MAIRCAResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "COPRAS" begin

        tol = 0.0001
        decmat = [
            2.50 240 57 45 1.10 0.333333
            2.50 285 60 75 4.00 0.428571
            4.50 320 100 65 7.50 1.111111
            4.50 365 100 90 7.50 1.111111
            5.00 400 100 90 11.00 1.111111
            2.50 225 60 45 1.10 0.333333
            2.50 270 57 60 4.00 0.428571
            4.50 330 100 70 7.50 1.111111
            4.50 365 100 80 7.50 1.111111
            5.00 380 110 65 8.00 1.111111
            2.50 285 65 80 4.00 0.400000
            4.00 280 75 65 4.00 0.400000
            4.50 365 102 95 7.50 1.111111
            4.50 400 102 95 7.50 1.111111
            6.00 450 110 95 11.00 1.176471
            6.00 510 110 105 11.00 1.176471
            6.00 330 140 110 18.50 1.395349
            2.50 240 65 80 4.00 0.400000
            4.00 280 75 75 4.00 0.400000
            4.50 355 102 95 7.50 1.111111
            4.50 385 102 90 7.50 1.111111
            5.00 385 114 95 7.50 1.000000
            6.00 400 110 90 11.00 1.000000
            6.00 480 110 95 15.00 1.000000
            6.00 440 140 100 18.50 1.200000
            6.00 500 140 100 18.50 1.200000
            5.00 450 125 100 15.00 1.714286
            6.00 500 150 125 18.50 1.714286
            6.00 515 180 140 22.00 2.307692
            7.00 550 200 150 30.00 2.307692
            6.00 500 180 140 15.00 2.307692
            6.00 500 180 140 18.50 2.307692
            6.00 500 180 140 22.00 2.307692
            7.00 500 180 140 30.00 2.307692
            7.00 500 200 140 37.00 2.307692
            7.00 500 200 140 45.00 2.307692
            7.00 500 200 140 55.00 2.307692
            7.00 500 200 140 75.00 2.307692
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.1667, 0.1667, 0.1667, 0.1667, 0.1667, 0.1667]

        fns = [maximum, maximum, maximum, maximum, maximum, minimum]

        result = copras(df, weights, fns)
        @test result isa COPRASResult
        @test isapprox(
            result.scores,
            [
                0.44194,
                0.44395,
                0.41042,
                0.44403,
                0.48177,
                0.44074,
                0.42430,
                0.41737,
                0.43474,
                0.44382,
                0.46625,
                0.48602,
                0.45019,
                0.45825,
                0.51953,
                0.54265,
                0.56134,
                0.45588,
                0.49532,
                0.44788,
                0.45014,
                0.48126,
                0.51586,
                0.56243,
                0.58709,
                0.60091,
                0.51850,
                0.61085,
                0.65888,
                0.75650,
                0.61430,
                0.63486,
                0.65542,
                0.72065,
                0.77680,
                0.82379,
                0.88253,
                1.00000,
            ],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = copras(setting)
        @test result2 isa COPRASResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, CoprasMethod())
        @test result3 isa COPRASResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "PROMETHEE" begin
        tol = 0.005
        decmat = [
            42.0 35 43 51
            89 72 92 85
            14 85 17 40
            57 60 45 80
            48 32 43 40
            71 45 60 85
            69 40 72 55
            64 35 70 60
        ]
        df = makeDecisionMatrix(decmat)
        qs = [49, nothing, 45, 30]
        ps = [100, 98, 95, 80]
        weights = [0.25, 0.35, 0.22, 0.18]
        fns = makeminmax([maximum, maximum, maximum, maximum])
        prefs = convert(
            Array{Function,1},
            [prometLinear, prometVShape, prometLinear, prometLinear],
        )

        result = promethee(df, weights, fns, prefs, qs, ps)

        @test result isa PrometheeResult
        @test isapprox(
            result.scores,
            [0.07, -0.15, -0.06, -0.05, 0.10, 0.0, 0.03, 0.06],
            atol = tol,
        )
        @test result.bestIndex == 5

        setting = MCDMSetting(df, weights, fns)
        result2 = promethee(setting, prefs, qs, ps)
        @test result2 isa PrometheeResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, PrometheeMethod(prefs, qs, ps))
        @test result3 isa PrometheeResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex

    end

    @testset "CoCoSo" begin
        tol = 0.0001
        decmat = [
            60.00 0.40 2540.00 500.00 990.00
            6.35 0.15 1016.00 3000.00 1041.00
            6.80 0.10 1727.20 1500.00 1676.00
            10.00 0.20 1000.00 2000.00 965.00
            2.50 0.10 560.00 500.00 915.00
            4.50 0.08 1016.00 350.00 508.00
            3.00 0.10 1778.00 1000.00 920.00
        ]

        df = makeDecisionMatrix(decmat)

        weights = [0.036, 0.192, 0.326, 0.326, 0.120]

        lambda = 0.5

        fns = [maximum, minimum, maximum, maximum, maximum]

        result = cocoso(df, weights, fns, lambda = lambda)
        @test result isa CoCoSoResult
        @test isapprox(
            result.scores,
            [
                2.0413128390265998,
                2.787989783418825,
                2.8823497955972495,
                2.4160457689259287,
                1.2986918936013303,
                1.4431429073391682,
                2.519094173200623,
            ],
            atol = tol,
        )

        setting = MCDMSetting(df, weights, fns)
        result2 = cocoso(setting)
        @test result2 isa CoCoSoResult
        @test result2.scores == result.scores
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, CocosoMethod())
        @test result3 isa CoCoSoResult
        @test result3.scores == result.scores
        @test result3.bestIndex == result.bestIndex
    end

    @testset "CRITIC" begin
        tol = 0.0001
        decmat = [
            12.9918 0.7264 -1.1009 1.598139592
            4.1201 5.8824 3.4483 1.021563567
            4.1039 0.0000 -0.5076 0.984469444
        ]

        df = makeDecisionMatrix(decmat)

        fns = [maximum, maximum, minimum, maximum]

        result = critic(df, fns)
        @test result isa CRITICResult
        @test isapprox(
            result.w,
            [0.16883925, 0.418444976, 0.249124763, 0.163591012],
            atol = tol,
        )

        setting = MCDMSetting(df, zeros(4), fns)
        result2 = critic(setting)
        @test result2 isa CRITICResult
        @test result2.ranking == result.ranking
        @test result2.bestIndex == result.bestIndex

        result3 = mcdm(setting, CriticMethod())
        @test result3 isa CRITICResult
        @test result3.ranking == result.ranking
        @test result3.bestIndex == result.bestIndex
    end

    @testset "Entropy" begin
        tol = 0.0001
        df = DataFrame(
            C1 = [2, 4, 3, 5, 4, 3],
            C2 = [1, 1, 2, 1, 2, 2],
            C3 = [4, 5, 6, 5, 5, 6],
            C4 = [7, 6, 6, 7, 6, 6],
            C5 = [6, 7, 5, 6, 7, 6],
            C6 = [6, 7, 6, 7, 7, 6],
            C7 = [7, 6, 8, 7, 6, 6],
            C8 = [3000, 3500, 4000, 3000, 3000, 3500],
        )

        result = entropy(df)

        @test result isa EntropyResult

        @test isapprox(
            result.w,
            [
                0.29967360960,
                0.44136733892,
                0.07009088720,
                0.02123823711,
                0.04902292895,
                0.02308037885,
                0.04776330969,
                0.04776330969,
            ],
            atol = tol,
        )

    end

    @testset "CODAS" begin
        tol = 0.0001
        decmat = [
            60.000 0.400 2540 500 990
            6.350 0.150 1016 3000 1041
            6.800 0.100 1727.2 1500 1676
            10.000 0.200 1000 2000 965
            2.500 0.100 560 500 915
            4.500 0.080 1016 350 508
            3.000 0.100 1778 1000 920
        ]

        df = makeDecisionMatrix(decmat)

        w = [0.036, 0.192, 0.326, 0.326, 0.12]
        fns = [maximum, minimum, maximum, maximum, maximum]

        result = codas(df, w, fns)
        @test result isa CODASResult
        @test isapprox(
            result.scores,
            [
                0.512176491,
                1.463300035,
                1.07153259,
                -0.212467998,
                -1.851520552,
                -1.17167677,
                0.188656204,
            ],
            atol = tol,
        )

        setting = MCDMSetting(df, w, fns)
        result2 = codas(setting)
        @test result2 isa CODASResult
        @test result2.scores == result.scores

        result3 = mcdm(setting, CodasMethod())
        @test result3 isa CODASResult
        @test result3.scores == result.scores
    end
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
module ReverseDiff

using Base: RefValue
using Random
using LinearAlgebra
using Statistics

using FunctionWrappers: FunctionWrapper

using DiffResults
using DiffResults: DiffResult
using DiffRules, SpecialFunctions, NaNMath

using ForwardDiff
using ForwardDiff: Dual, Partials
using StaticArrays

# Not all operations will be valid over all of these types, but that's okay; such cases
# will simply error when they hit the original operation in the overloaded definition.
const ARRAY_TYPES = (:AbstractArray, :AbstractVector, :AbstractMatrix, :Array, :Vector, :Matrix)
const REAL_TYPES = (:Bool, :Integer, :(Irrational{:e}), :(Irrational{:π}), :Rational, :BigFloat, :BigInt, :AbstractFloat, :Real, :Dual)

const SKIPPED_UNARY_SCALAR_FUNCS  = Symbol[:isinf, :isnan, :isfinite, :iseven, :isodd, :isreal, :isinteger]
const SKIPPED_BINARY_SCALAR_FUNCS = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]

include("tape.jl")
include("tracked.jl")
include("macros.jl")
include("derivatives/propagation.jl")
include("derivatives/scalars.jl")
include("derivatives/elementwise.jl")
include("derivatives/linalg/arithmetic.jl")
include("derivatives/linalg/reductions.jl")
include("derivatives/linalg/special.jl")
include("api/utils.jl")
include("api/Config.jl")
include("api/tape.jl")
include("api/gradients.jl")
include("api/jacobians.jl")
include("api/hessians.jl")

export DiffResults

end # module


text:
function kmr_markov_matrix_sequential{T<:Real}(n::Integer, p::T, ε::T)
    """
    Generate the MarkovChain with the associated transition matrix from the KMR model with *sequential* move

    n: number of players
    p: level of p-dominance for action 1
       = the value of p such that action 1 is the BR for (1-q, q) for any q > p,
         where q (1-q, resp.) is the prob that the opponent plays action 1 (0, resp.)
    ε: mutation probability

    References:
        KMRMarkovMatrixSequential is contributed from https://github.com/oyamad
    """
    x = zeros(T, n+1, n+1)

    x[1, 1], x[1, 2] = 1 - ε/2, ε/2
    @inbounds for i = 1:n-1
        x[i+1, i] = (i/n) * (ε/2 + (1 - ε) *
                             (((i-1)/(n-1) < p) + ((i-1)/(n-1) == p)/2))
        x[i+1, i+2] = ((n-i)/n) * (ε/2 + (1 - ε) *
                                 ((i/(n-1) > p) + (i/(n-1) == p)/2))
        x[i+1, i+1] = 1 - x[i+1, i] - x[i+1, i+2]
    end
    x[end, end-1], x[end, end] = ε/2, 1 - ε/2
    return MarkovChain(x)
end

@testset "Testing mc_tools.jl" begin

    # these matrices come from RMT4 section 2.2.1
    mc1 = [1 0 0; .2 .5 .3; 0 0 1]
    mc2 = [.7 .3 0; 0 .5 .5; 0 .9 .1]
    mc3 = [0.4 0.6; 0.2 0.8]
    mc4 = eye(2)
    mc5 = [
         0 1 0 0 0 0
         1 0 0 0 0 0
         1//2 0 0 1//2 0 0
         0 0 0 0 1 0
         0 0 0 0 0 1
         0 0 0 1 0 0
         ]
    mc5_stationary = zeros(Rational,6,2)
    mc5_stationary[[1,2]] = 1//2; mc5_stationary[[10,11,12]] = 1//3

    mc6 = [2//3 1//3; 1//4 3//4]  # Rational elements
    mc6_stationary = [3//7, 4//7]

    mc7 = [1 0; 0 1]
    mc7_stationary = [1 0;0 1]

    # Reducible mc with a unique recurrent class,
    # where n=2 is a transient state
    mc10 = [1. 0; 1. 0]
    mc10_stationary = [1., 0]

    mc1 = MarkovChain(mc1)
    mc2 = MarkovChain(mc2)
    mc3 = MarkovChain(mc3)
    mc4 = MarkovChain(mc4)
    mc5 = MarkovChain(mc5)
    mc6 = MarkovChain(mc6)
    mc7 = MarkovChain(mc7)
    mc10 = MarkovChain(mc10)

    # examples from
    # Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

    fig1_p = zeros(Rational{Int64}, 5, 5)
    fig1_p[[3, 4, 9, 10, 11, 13, 18, 19, 22, 24]] =
        [1//2, 2//5, 1//10, 1, 1, 1//5, 3//10, 1//5, 1, 3//10]
    fig2_p = zeros(Rational{Int64}, 5, 5)
    fig2_p[[3, 10, 11, 13, 14, 17, 18, 19, 22]] =
        [1//3, 1, 1, 1//3, 1//2, 1//2, 1//3, 1//2, 1//2]

    fig1 = MarkovChain(convert(Matrix{Float64}, fig1_p))
    fig1_rat = MarkovChain(fig1_p)

    fig2 = MarkovChain(convert(Matrix{Float64}, fig2_p))
    fig2_rat = MarkovChain(fig2_p)

    mc8 = kmr_markov_matrix_sequential(27, 1/3, 1e-2)
    mc9 = kmr_markov_matrix_sequential(3, 1/3, 1e-14)

    tol = 1e-15

    @testset "test mc_compute_stationary using exact solutions" begin
        @test mc_compute_stationary(mc1) == eye(3)[:, [1, 3]]
        @test isapprox(mc_compute_stationary(mc2), [0, 9/14, 5/14])
        @test isapprox(mc_compute_stationary(mc3), [1/4, 3/4])
        @test mc_compute_stationary(mc4) == eye(2)
        @test mc_compute_stationary(mc5) == mc5_stationary
        @test mc_compute_stationary(mc6) == mc6_stationary
        @test mc_compute_stationary(mc7) == mc7_stationary
        @test mc_compute_stationary(mc10) == mc10_stationary
    end

    @testset "test gth_solve with KMR matrices" begin
        for d in [mc8,mc9]
            x = mc_compute_stationary(d)

            # Check elements sum to one
            @test isapprox(sum(x), 1; atol=tol)

            # Check elements are nonnegative
            for i in 1:length(x)
                @test x[i] >= -tol
            end

            # Check x is a left eigenvector of P
            @test isapprox(vec(x'*d.p), x; atol=tol)
        end
    end

    @testset "test MarkovChain throws errors" begin
        # not square
        @test_throws DimensionMismatch MarkovChain(rand(4, 5))

        # state_values to long
        @test_throws DimensionMismatch MarkovChain([0.5 0.5; 0.2 0.8], rand(3))

        # first row doesn't sum to 1
        @test_throws ArgumentError MarkovChain([0.0 0.5; 0.2 0.8])

        # negative element, but sums to 1
        @test_throws ArgumentError MarkovChain([-1 1; 0.2 0.8])
    end

    @testset "test graph theoretic algorithms" begin
        for fig in [fig1, fig1_rat]
            @test recurrent_classes(fig) == Vector[[2, 5]]
            @test communication_classes(fig) == Vector[[2, 5], [1, 3, 4]]
            @test is_irreducible(fig) == false
            # @test period(fig) == 2
            # @test is_aperiodic(fig) == false
        end

        for fig in [fig2, fig2_rat]
            @test recurrent_classes(fig) == Vector[[1, 3, 4]]
            @test communication_classes(fig) == Vector[[1, 3, 4], [2, 5]]
            @test is_irreducible(fig) == false
            # @test period(fig) == 1
            # @test is_aperiodic(fig) == true
        end
    end

    @testset "test simulate shape" begin

        for mc in (mc3, MarkovChain(sparse(mc3.p)))
            ts_length = 10
            init = [1, 2]
            nums_reps = [3, 1]

            @test size(simulate(mc, ts_length)) == (ts_length,1)
            @test size(simulate(mc, ts_length, init)) ==
                (ts_length, length(init))
            num_reps = nums_reps[1]
            @test size(simulate(mc, ts_length, init, num_reps=num_reps)) ==
                (ts_length, length(init)*num_reps)
            for num_reps in nums_reps
            @test size(simulate(mc, ts_length, num_reps=num_reps)) ==
                    (ts_length, num_reps)
            end
        end
    end  # testset

    @testset "test simulate init array" begin
        for mc in (mc3, MarkovChain(sparse(mc3.p)))
            ts_length = 10
            init = [1, 2]
            num_reps = 3

            X = simulate(mc, ts_length, init, num_reps=num_reps)
            @test vec(X[1, :]) == repmat(init, num_reps)
        end
    end  # testset

    @testset "Sparse Matrix" begin
        # use fig1_p from above because we already checked the graph theory
        # algorithms for this stochastic matrix.
        p = fig1_p
        p_s = sparse(p)

        @testset "constructors" begin
            mc_s = MarkovChain(p_s)
            @test isa(mc_s, MarkovChain{eltype(p), typeof(p_s)})
        end

        mc_s = MarkovChain(p_s)
        mc = MarkovChain(p)

        @testset "basic correspondence with dense version" begin
            @test n_states(mc) == n_states(mc_s)
            @test maxabs(mc.p - mc_s.p) == 0.0
            @test recurrent_classes(mc) == recurrent_classes(mc_s)
            @test communication_classes(mc) == communication_classes(mc_s)
            @test is_irreducible(mc) == is_irreducible(mc_s)
            @test is_aperiodic(mc) == is_aperiodic(mc_s)
            @test period(mc) == period(mc_s)
            @test maxabs(gth_solve(mc_s.p) - gth_solve(mc.p)) < 1e-15
        end
    end

    @testset "simulate_values" begin
        for _mc in (mc3, MarkovChain(sparse(mc3.p)))
            for T in [Float16, Float32, Float64, Int8, Int16, Int32, Int64]
                mc = MarkovChain(_mc.p, rand(T, size(_mc.p, 1)))
                ts_length = 10
                init = [1, 2]
                nums_reps = [3, 1]

                @test size(@inferred(simulate_values(mc, ts_length))) == (ts_length,1)
                @test size(@inferred(simulate_values(mc, ts_length, init))) ==
                    (ts_length, length(init))
                num_reps = nums_reps[1]
                @test size(simulate(mc, ts_length, init, num_reps=num_reps)) ==
                    (ts_length, length(init)*num_reps)
                for num_reps in nums_reps
                    @test size(simulate(mc, ts_length, num_reps=num_reps)) ==
                        (ts_length, num_reps)
                end
            end  # state_values eltypes
        end

        @testset "test simulate_values init array" begin
            for mc in (mc3, MarkovChain(sparse(mc3.p)))
                ts_length = 10
                init = [1, 2]
                num_reps = 3

                X = simulate(mc, ts_length, init, num_reps=num_reps)
                @test vec(X[1, :]) == repmat(init, num_reps)
            end
        end  # testset
    end

    @testset "simulation and value_simulation" begin
        @test size(@inferred(simulation(mc3, 10, 1))) == (10,)
        @test size(@inferred(value_simulation(mc3, 10, 1))) == (10,)
    end

end  # testset


text:
using SHA

include("config.jl")

# BinaryBuilder => Rust target triple name matching
RUST_TRIPLES = Dict(
    "i686-linux-gnu" => "i686-unknown-linux-gnu",
    "x86_64-linux-gnu" => "x86_64-unknown-linux-gnu",
    "aarch64-linux-gnu" => "aarch64-unknown-linux-gnu",
    "armv7l-linux-gnueabihf" => "armv7-unknown-linux-gnueabihf",
    "powerpc64le-linux-gnu" => "powerpc64le-unknown-linux-gnu",
    "i686-linux-musl" => "i686-unknown-linux-musl",
    "x86_64-linux-musl" => "x86_64-unknown-linux-musl",
    "aarch64-linux-musl" => "aarch64-unknown-linux-musl" ,
    "armv7l-linux-musleabihf" => "armv7-unknown-linux-musleabihf",
    "x86_64-apple-darwin" => "x86_64-apple-darwin",
    "aarch64-apple-darwin" => "aarch64-apple-darwin",
    "x86_64-unknown-freebsd" => "x86_64-unknown-freebsd",
    "i686-w64-mingw32" => "i686-pc-windows-gnu",
    "x86_64-w64-mingw32" => "x86_64-pc-windows-gnu",
)

prebuilt = []
prefix = "chemfiles-static.v$version."
suffix = ".tar.gz"
for path in readdir("products/")
    if !startswith(path, "chemfiles-static")
        continue
    end

    julia_triple = path[length(prefix) + 1:end - length(suffix)]

    fd = open(joinpath("products/", path))
    sha = bytes2hex(sha256(fd))
    close(fd)

    push!(prebuilt, (RUST_TRIPLES[julia_triple], julia_triple, sha))
end

open("prebuilt.rs", "w") do fd
    println(fd, "/// Get the julia triple & sha256 corresponding to the prebuilt libchemfiles")
    println(fd, "/// for a given rust triple, if it exists")
    println(fd, "pub fn get_prebuilt_info(target: &str) -> Option<(&'static str, &'static str)> {")
    println(fd, "    match target {")
    for (rust_triple, julia_triple, sha) in prebuilt
        println(fd, "        \"$rust_triple\" => Some((\"$julia_triple\", \"$sha\")),")
    end
    println(fd, "        _ => None,")
    println(fd, "    }")
    println(fd, "}")
end

println("====================================================\n")
println("Done generating prebuilt.rs")

max_stars_repo_name: chemfiles/chemfiles-prebuilt

max_issues_repo_name: chemfiles/chemfiles-prebuilt

max_forks_repo_name: chemfiles/chemfiles-prebuilt


text: @test_deprecated IC.debug(NumberLimit(2))


text:
module Computer
export Program, read_program, run_program!, step_program!, reset!

struct Memory
    vect :: Vector{Int}
    dict :: Dict{Int, Int}
end

Memory(v::Vector{Int}) = Memory(v, Dict())

function Base.getindex(m::Memory, i::Int)
    if i <= length(m.vect)
        @inbounds m.vect[i]
    else
        get(m.dict, i, 0)
    end
end

function Base.setindex!(m::Memory, val::Int, i::Int)
    if i <= length(m.vect)
        @inbounds m.vect[i] = val
    else
        m.dict[i] = val
    end
end

mutable struct Program
    mem      :: Memory
    pointer  :: Int
    code     :: Int
    rel_base :: Int

    input    :: Int
    output   :: Int
end

Program(p) = Program(Memory(p), 0,0,0,0,0) |> reset!

read_program(fname) = open(fname) do f
    parse.(Int, split(readline(f), ",")) |> Program
end

debug() = false

macro read_args(N)
    ret = Expr(:block)
    for i in 1:N
        r = Symbol("r_", i)
        x = Symbol("x_", i)
        param = Symbol("param_", i)
        block = quote
            $x = prog.mem[prog.pointer]
            if     $param == 0
                $r = $x
                $x = prog.mem[1+$r]
            elseif $param == 1
                $r = -1
            elseif $param == 2
                $r = $x + prog.rel_base
                $x = prog.mem[1+$r]
            end
            prog.pointer += 1
        end |> esc
        push!(ret.args, block)
    end
    ret
end

function run_program!(p::Vector{Int}, inputs=(0))
    p = Program(p)
    run_program!(p, inputs)
    p
end

function run_program!(prog::Program, inputs=(0))
    outputs = Int[]
    for input in inputs
        prog.input = input
        while true
            step_program!(prog)
            if     prog.code == 3
                break # feed next input
            elseif prog.code == 4
                push!(outputs, prog.output)
            elseif prog.code == 99
                return outputs
            else
                error("unexpected program stopping code: $(prog.code)")
            end
        end
    end
end

function step_program!(prog::Program)
    input_available = true

    while true
        pointer_reset = prog.pointer
        instr = prog.mem[prog.pointer]
        (instr, prog.code) = divrem(instr, 100)
        Base.Cartesian.@nexprs 3 j -> begin
            (instr, param_j) = divrem(instr, 10)
        end

        if debug()
            @show prog.mem
            @show prog.code, param_1, param_2, param_3
            @show prog.pointer
            println()
        end

        prog.pointer += 1

        # Add
        if prog.code == 1
            @read_args 3
            @assert r_3 != -1
            prog.mem[1+r_3] = x_1 + x_2
            continue
        end

        # Mul
        if prog.code == 2
            @read_args 3
            @assert r_3 != -1
            prog.mem[1+r_3] = x_1 * x_2
            continue
        end

        # Input
        if prog.code == 3
            @read_args 1
            @assert r_1 != -1
            if input_available
                prog.mem[1+r_1] = prog.input
                input_available = false
            else
                prog.pointer = pointer_reset
                return prog
            end
            continue
        end

        # Output
        if prog.code == 4
            @read_args 1
            prog.output = x_1
            return prog
        end

        # Jump non-zero
        if prog.code == 5
            @read_args 2
            if x_1 != 0
                prog.pointer = 1+x_2
            end
            continue
        end

        # Jump zero
        if prog.code == 6
            @read_args 2
            if x_1 == 0
                prog.pointer = 1+x_2
            end
            continue
        end

        # Less than
        if prog.code == 7
            @read_args 3
            @assert r_3 != -1
            prog.mem[1+r_3] = x_1<x_2 ? 1 : 0
            continue
        end

        # Equals
        if prog.code == 8
            @read_args 3
            @assert r_3 != -1
            prog.mem[1+r_3] = x_1==x_2 ? 1 : 0
            continue
        end

        # Adjust relative base
        if prog.code == 9
            @read_args 1
            prog.rel_base += x_1
            continue
        end

        # Exit
        @assert prog.code == 99
        return prog
    end

    nothing
end

function reset!(dest::Program, src::Program)
    copy!(dest.mem.vect, src.mem.vect)
    copy!(dest.mem.dict, src.mem.dict)
    reset!(dest)
end

function reset!(dest::Program)
    dest.pointer = 1
    dest.code = 0

    dest.input  = typemin(Int)
    dest.output = typemax(Int)
    dest
end

end


text:
using ACME
using Test

include("..\\src\\AudioToaster.jl")

circuitfrommacro = @circuit begin
    r1 = resistor(1e4)
    c1 = capacitor(4.7e-9), [1] ⟷ r1[2], [2] ⟷ gnd
    src = voltagesource(), [+] ⟷ r1[1], [-] ⟷ gnd
    prb = voltageprobe(), [+] ⟷ c1[1], [-] ⟷ gnd
end
circuitfromfile = AudioToaster.loadfile("examples\\simplest.diy")
modelfrommacro = DiscreteModel(circuitfrommacro, 1/44100)
modelfromfile = DiscreteModel(circuitfromfile, 1/44100)
testdata = reshape(rand(1000) .* 2.0 .- 1.0, 1, :)
expected = run!(modelfrommacro, testdata)
actual = run!(modelfromfile, testdata)

for i in 1:length(expected)
    @test expected[i] ≈ actual[i]
end

max_stars_repo_name: stellartux/diyloader

max_issues_repo_name: stellartux/diyloader

max_forks_repo_name: stellartux/diyloader


text:
export RealTimeDynamicPlanner, RTDP

"Planner that uses Real Time Dynamic Programming (RTDP)."
@kwdef mutable struct RealTimeDynamicPlanner <: Planner
    heuristic::Heuristic = GoalCountHeuristic()
    n_rollouts::Int = 50
    max_depth::Int = 50
    rollout_noise::Float64 = 0.0
    action_noise::Float64 = 0.0
end

const RTDP = RealTimeDynamicPlanner

function solve(planner::RealTimeDynamicPlanner,
               domain::Domain, state::State, spec::Specification)
    # Initialize then refine solution
    sol = PolicyValue()
    sol.V[hash(state)] = -planner.heuristic(domain, state, spec)
    sol = solve!(planner, sol, domain, state, spec)
    # Wrap in Boltzmann policy if needed
    return planner.action_noise == 0 ?
        sol : BoltzmannPolicy(sol, planner.action_noise)
end

function solve!(planner::RealTimeDynamicPlanner, sol::PolicyValue,
                domain::Domain, state::State, spec::Specification)
    @unpack heuristic, action_noise = planner
    @unpack n_rollouts, max_depth, rollout_noise = planner
    # Precompute heuristic information
    precompute!(heuristic, domain, state, spec)
    # Perform rollouts from initial state
    initial_state = state
    ro_policy = rollout_noise == 0 ? sol : BoltzmannPolicy(sol, rollout_noise)
    visited = State[]
    for n in 1:n_rollouts
        state = initial_state
        # Rollout until maximum depth
        for t in 1:max_depth
            push!(visited, state)
            if is_goal(spec, domain, state) break end
            update_values!(planner, sol, domain, state, spec)
            act = get_action(ro_policy, state)
            state = transition(domain, state, act)
        end
        # Post-rollout update
        while length(visited) > 0
            state = pop!(visited)
            update_values!(planner, sol, domain, state, spec)
        end
    end
    return sol
end

function solve!(planner::RealTimeDynamicPlanner,
                sol::BoltzmannPolicy{PolicyValue},
                domain::Domain, state::State, spec::Specification)
    sol = solve!(planner, sol.policy, domain, state, spec)
    return BoltzmannPolicy(sol, planner.action_noise)
end

function update_values!(planner::RealTimeDynamicPlanner, sol::PolicyValue,
                        domain::Domain, state::State, spec::Specification)
    actions = available(state, domain)
    state_id = hash(state)
    if is_goal(spec, domain, state)
        qs = zeros(length(actions))
        sol.Q[state_id] = Dict{Term,Float64}(zip(actions, qs))
        sol.V[state_id] = 0.0
        return
    end
    qs = map(actions) do act
        next_state = transition(domain, state, act)
        r = get_reward(spec, domain, state, act, next_state)
        h_val = planner.heuristic(domain, next_state, spec)
        return get_discount(spec) * get!(sol.V, hash(next_state), -h_val) + r
    end
    sol.Q[state_id] = Dict{Term,Float64}(zip(actions, qs))
    sol.V[state_id] = planner.action_noise == 0 ?
        maximum(qs) : sum(softmax(qs ./ planner.action_noise) .* qs)
end


text:
using PyPlot, SetPyPlot
using Base.Test

# write your own tests here

setrcparam()

figure(1)

x = collect(linspace(0, 2π))
y1 = sin(x)
plot(x, sin(x), label = "sin(x)")
plot(x, cos(x), label = "cos(x)")
xlabel("x values")
ylabel("y values")
legend()
tight_layout()

gcf() # to plot in vscode

# savefig("publication_test.png")

@test 1==1