text:
# generate the typical outcome plots found in
# https://arxiv.org/abs/2201.05967

using DyadicKDE
include("./plot_helpers.jl")

# specify parameters
n_data = 100
n_evals = 100
degeneracies = ["total", "partial", "none"]
kernel_name = "epanechnikov_order_2"
evals = collect(range(-2.0, stop=2.0, length=n_evals))
sdp_solver = "cosmo"
n_resample = 10000
significance_level = 0.05
results = []
ps = Dict(
    "total" => [0.5, 0.0, 0.5],
    "partial" => [0.25, 0.0, 0.75],
    "none" => [0.2, 0.2, 0.6],
)

println("Running experiments")
for degen in degeneracies

    println("Degeneracy: $degen")
    data = make_data(n_data, ps[degen])
    h_ROT = estimate_ROT_bandwidth(data, "epanechnikov_order_2")

    est = DyadicKernelDensityEstimator(
        kernel_name, h_ROT, significance_level,
        n_resample, sdp_solver, evals, data,
        Dict("degen" => degen, "p" => ps[degen]))

    fit(est)
    push!(results, est)

end

println("Making plots")
linewidth = 1.0
y_lim = [-0.01, 0.45]
e = 0.22
pci_marker = PyPlot.matplotlib.path.Path([[-1,0], [-1,e], [-1,-e], [-1,0], [1,0], [1,e], [1,-e], [1,0]])
handle_f = PyPlot.matplotlib.lines.Line2D([0], [0], color="k", lw=linewidth, label="\$f_W(w)\$")
handle_fhat = PyPlot.matplotlib.lines.Line2D([0], [0], color="k", lw=linewidth, linestyle=(0, (1,1)), label="\$\\widehat f_W(w)\$")
handle_ucb = PyPlot.matplotlib.patches.Patch(facecolor="lightgray", edgecolor="lightgray", label="UCB")
handle_pci = PyPlot.matplotlib.lines.Line2D([0], [0], color="black", lw=0, markersize=17, marker=pci_marker, label="PCI")
handles = [handle_f, handle_fhat, handle_ucb, handle_pci]

for est in results

    degen = est.meta["degen"]
    p = est.meta["p"]
    f = generate_f(est.evals, p)

    # plot f
    fig, ax = plt.subplots(figsize=(4,4))
    ax.plot(est.evals, f, color = "black", linewidth=linewidth)

    # plot fhat
    ax.plot(est.evals, est.fhat,
        color = "black", linewidth=linewidth, linestyle=(0, (1,1)))

    # plot ucb
    ax.fill_between(est.evals, est.ucb[1,:], est.ucb[2,:],
        color="lightgray", linewidth=0.0)

    # plot pci
    plot_confidence_intervals(ax, est.evals, est.pci[1,:],
        est.pci[2,:], 10, "black", linewidth, "-")

    # save plot
    PyPlot.xlabel("Evaluation point, \$w\$")
    plt.ylim(y_lim)
    x_min = minimum(est.evals)
    x_max = maximum(est.evals)
    plt.xlim((x_min, x_max))
    plt.yticks(range(0.0, stop=0.4, step=0.1))
    legend(handles=handles, loc="upper left")
    plt.ylabel("Density", labelpad=4.0)
    plt.tight_layout()
    PyPlot.savefig("outcome_plot_$degen.pdf")
    close("all")

end


text:
struct eigOptions
    # structure to encapsulate all the different options for the scf iterations
    # the main idea is to define a Hamiltonian, and then run scf on it; however,
    # the Hamiltoinian shouldn't have the scf options embeded
    eigstol::Float64
    eigsiter::Int64
    eigmethod::AbstractString

    function eigOptions()
        new(1e-8, 100, "eigs")
    end

    function eigOptions(eigstol::Float64, eigsiter::Int64, eigmethod::String)
        new(eigstol, eigsiter, eigmethod)
    end
end

# defining abstract type for all the different mixing options
abstract type mixingOptions end

mutable struct simpleMixOptions <: mixingOptions
    betamix::Float64
    function simpleMixOptions(betamix)
        new(betamix)
    end
    # default constructor
    function simpleMixOptions()
        new(0.5)
    end
end

mutable struct andersonMixOptions <: mixingOptions
    ymat::Array{Float64,2}
    smat::Array{Float64,2}
    betamix::Float64
    mixdim::Int64
    iter::Int64
    function andersonMixOptions(Ns, betamix, mixdim )
        ymat = zeros(Ns, mixdim);
        smat = zeros(Ns, mixdim);
        new(ymat,smat,betamix[1],mixdim,1)
    end
end

mutable struct andersonPrecMixOptions <: mixingOptions
    ymat::Array{Float64,2}
    smat::Array{Float64,2}
    betamix::Float64
    mixdim::Int64
    iter::Int64
    prec::Function
    precArgs
    function andersonPrecMixOptions(Ns, betamix, mixdim, prec, precArgs)
        ymat = zeros(Ns, scfOpts.mixdim);
        smat = zeros(Ns, scfOpts.mixdim);
        new(ymat,smat,betamix[1],mixdim,1, prec, precArgs)
    end
end

mutable struct kerkerMixOptions <: mixingOptions
    betamix::Float64
    KerkerB::Float64
    kx::Array{Float64,2}
    YukawaK::Float64
    epsil0::Float64
    function kerkerMixOptions(betamix,KerkerB, kx, YukawaK, epsil0)
        new(betamix,KerkerB, kx, YukawaK, epsil0)
    end
end

function updateMix!( mixOpts::andersonMixOptions, ii )
    mixOpts.iter = ii;
end

struct scfOptions
    # structure to encapsulate all the different options for the scf iterations
    # the main idea is to define a Hamiltonian, and then run scf on it; however,
    # the Hamiltoinian shouldn't have the scf options embeded
    SCFtol::Float64
    scfiter::Int64
    eigOpts::eigOptions
    mixOpts::mixingOptions

    # TODO: this Initializing is ugly, we need to use kwargs to properly initialize
    function scfOptions()
        eigOpts = eigOptions();
        mixOpts = simpleMixOptions();
        new(1e-7,100, eigOpts, mixOpts);
    end
    function scfOptions(SCFtol::Float64,scfiter::Int64,
                        eigOpts::eigOptions,
                        mixOpts::mixingOptions)
        new(SCFtol,scfiter, eigOpts, mixOpts)
    end
    function scfOptions(eigOpts::eigOptions, mixOpts::mixingOptions)
        new(1e-7,100,  eigOpts, mixOpts)
    end

end

ext: jl

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text: @time f(Cosine(2), x, y)


text:
"""
$(README)
"""
module MATDaemon

import DaemonMode
import MAT
import MacroTools
import Pkg

using DocStringExtensions: README, TYPEDFIELDS, TYPEDSIGNATURES
using Downloads: download

export download_jlcall

# Options file for communicating with MATLAB
const JL_OPTIONS = "jlcall_opts.mat"

# Containers which are matlabified by iterating over `pairs`
const KeyValueContainer = Union{<:AbstractDict, <:NamedTuple, <:Base.Iterators.Pairs}

"""
    $(TYPEDSIGNATURES)

Download or copy the `jlcall.m` MATLAB API function to the file `filename`.

By default, `jlcall.m` is copied into the current working directory from the `api` subfolder of the installed `MATDaemon.jl` source tree.
The latest version of `jlcall.m` can instead be downloaded from GitHub by passing `latest = true`.
The destination `filename` can be overwritten if it exists by passing `force = true`.
"""
function download_jlcall(filename::String = joinpath(pwd(), "jlcall.m"); latest::Bool = false, force::Bool = false)
    if latest
        jlcall_github = "https://raw.githubusercontent.com/jondeuce/MATDaemon.jl/master/api/jlcall.m"
        jlcall_local = download(jlcall_github)
    else
        jlcall_local = normpath(@__DIR__, "../api/jlcall.m")
    end
    cp(jlcall_local, filename; force = force)
end

"""
    matlabify(x)

Convert Julia value `x` to equivalent MATLAB representation.
"""
matlabify(x) = x # default
matlabify(::Nothing) = zeros(Float64, 0, 0) # represent `Nothing` as MATLAB's `[]`
matlabify(::Missing) = zeros(Float64, 0, 0) # represent `Missing` as MATLAB's `[]`
matlabify(x::Symbol) = string(x)
matlabify(xs::Tuple) = matlabify_iterable(xs)
matlabify(xs::KeyValueContainer) = matlabify_pairs(xs)

matlabify_iterable(xs) = Any[matlabify(x) for x in xs]
matlabify_pairs(xs) = Dict{String, Any}(string(k) => matlabify(v) for (k, v) in pairs(xs))

matlabify_output(x) = Any[matlabify(x)] # default
matlabify_output(::Nothing) = Any[]
matlabify_output(x::Tuple) = Any[map(matlabify, x)...]

juliafy_args(xs::Array{Any}) = vec(xs)
juliafy_kwargs(xs::Dict{String, Any}) = Pair{Symbol, Any}[Symbol(k) => v for (k, v) in xs]

"""
    JLCallOptions(; kwargs...)

Julia struct for storing [`jlcall.m`](https://github.com/jondeuce/MATDaemon.jl/blob/master/api/jlcall.m) input parser results.

Struct fields/keyword arguments for constructor:

$(TYPEDFIELDS)
"""
Base.@kwdef struct JLCallOptions
    "User function to be parsed and evaluated"
    f::String               = "(args...; kwargs...) -> nothing"
    "MATLAB `.mat` file containing positional and keyword arguments for calling `f`"
    infile::String          = tempname(; cleanup = true) * ".mat"
    "MATLAB `.mat` file for writing outputs of `f` into"
    outfile::String         = tempname(; cleanup = true) * ".mat"
    "Julia runtime binary location"
    runtime::String         = joinpath(Base.Sys.BINDIR, "julia")
    "Julia project to activate before calling `f`"
    project::String         = ""
    "Number of Julia threads"
    threads::Int            = Threads.nthreads()
    "Julia setup script to include before defining and calling the user function"
    setup::String           = ""
    "Treat `f` as a generic Julia expression, not a function: evaluate `f` and return `nothing`"
    nofun::Bool             = false
    "Julia modules to import before defining and calling the user function"
    modules::Vector{Any}    = Any[]
    "Current working directory. Change path to this directory before loading code"
    cwd::String             = pwd()
    "MATDaemon workspace. Local Julia project and temporary files for communication with MATLAB are stored here"
    workspace::String       = mktempdir(; prefix = ".jlcall_", cleanup = true)
    "Start Julia instance on a local server using `DaemonMode.jl`"
    server::Bool            = true
    "Port to start Julia server on"
    port::Int               = 3000
    "Julia code is loaded into a persistent server environment if true. Otherwise, load code in unique namespace"
    shared::Bool            = true
    "Restart the Julia server before loading code"
    restart::Bool           = false
    "Shut down the julia server and return"
    shutdown::Bool          = false
    "Garbage collect temporary files after each call"
    gc::Bool                = true
    "Print debugging information"
    debug::Bool             = false
end

matlabify(opts::JLCallOptions) = Dict{String, Any}(string(k) => matlabify(getproperty(opts, k)) for k in fieldnames(JLCallOptions))

"""
    $(TYPEDSIGNATURES)

Start Julia server.
"""
function start(port::Int; shared::Bool, verbose::Bool = false)
    DaemonMode.serve(port, shared; print_stack = true, async = true, threaded = false)
    verbose && println("* Julia server started\n")
    return nothing
end

"""
    $(TYPEDSIGNATURES)

Kill Julia server. If server is already killed, do nothing.
"""
function kill(port::Int; verbose::Bool = false)
    try
        DaemonMode.sendExitCode(port)
        verbose && println("* Julia server killed\n")
    catch e
        if (e isa Base.IOError)
            #TODO: Check for proper error code:
            #   on linux: abs(e.code) == abs(Libc.ECONNREFUSED)
            #   on windows: ?
            verbose && println("* Julia server inactive; nothing to kill\n")
        else
            rethrow()
        end
    end
    return nothing
end

"""
    $(TYPEDSIGNATURES)

Load [`jlcall.m`](https://github.com/jondeuce/MATDaemon.jl/blob/master/api/jlcall.m) input parser results from `workspace`.
"""
function load_options(workspace::String)
    clean_value(k, v) =
        k == "modules" ? vec(v) : # MATLAB vectors are passed as column matrices
        k == "threads" && v == "auto" ? Threads.nthreads() : # Replace --threads=auto with threads from this session
        v
    mxopts = MAT.matread(joinpath(workspace, JL_OPTIONS))
    kwargs = Dict{Symbol, Any}(Symbol(k) => clean_value(k, v) for (k, v) in mxopts)
    opts = JLCallOptions(; kwargs..., workspace = workspace)
    return opts
end

"""
    $(TYPEDSIGNATURES)

Initialize [`jlcall`](@ref) environment.
"""
function init_environment(opts::JLCallOptions)
    # Change to current MATLAB working directory
    cd(opts.cwd)

    # Activate user project
    if !isempty(opts.project)
        proj_file = abspath(expanduser(opts.project))
        if isdir(proj_file)
            proj_file = joinpath(proj_file, "Project.toml")
        end
        if Base.active_project() != proj_file
            # Passed project is not active; activate it
            Pkg.activate(proj_file; io = devnull)
        end
    end

    return nothing
end

"""
    $(TYPEDSIGNATURES)

Save [`jlcall`](@ref) output results into workspace.
"""
function save_output(output, opts::JLCallOptions)
    # Try writing output to `opts.outfile`
    try
        MAT.matwrite(opts.outfile, Dict{String, Any}("output" => output))
    catch e
        println("* ERROR: Unable to write Julia output to .mat:\n*   ", opts.outfile, "\n")
        rm(opts.outfile; force = true)
        rethrow()
    end

    return nothing
end

"""
    $(TYPEDSIGNATURES)

Run Julia function `f` using [`jlcall.m`](https://github.com/jondeuce/MATDaemon.jl/blob/master/api/jlcall.m) input parser results `opts`.
"""
function jlcall(f::F, opts::JLCallOptions) where {F}
    # Since `f` is dynamically defined in a global scope, try to force specialization on `f` (may help performance)
    f_args = MAT.matread(opts.infile)
    args = juliafy_args(f_args["args"])
    kwargs = juliafy_kwargs(f_args["kwargs"])
    out = f(args...; kwargs...)

    # Save results to workspace
    save_output(matlabify_output(out), opts)
end

"""
    $(TYPEDSIGNATURES)

Location of script for loading code, importing modules, and evaluating the function expression passed from [`jlcall.m`](https://github.com/jondeuce/MATDaemon.jl/blob/master/api/jlcall.m).
"""
jlcall_script() = joinpath(@__DIR__, "..", "api", "jlcall.jl")

end # module MATDaemon


text:
# Hinge loss function

"""
  HingeLoss(b::Array{Real}, mu::Real=1.0)

Returns the function `g(x) = mu * sum(max(0, 1 - b_i * x_i), i=1,...,n )`.
"""

HingeLoss{T <: AbstractArray, R <: Real}(b::T, mu::R=1.0) = Postcompose(Precompose(SumPositive(), -b, 1.0), mu)


text:
module TestSDPA

using Test
using MathOptInterface
import SDPA

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_solver_name()
    @test MOI.get(SDPA.Optimizer(), MOI.SolverName()) == "SDPA"
end

function test_options()
    param = MOI.RawOptimizerAttribute("bad_option")
    err = MOI.UnsupportedAttribute(param)
    @test_throws err MOI.set(
        SDPA.Optimizer(),
        MOI.RawOptimizerAttribute("bad_option"),
        0,
    )
end

function test_runtests()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        MOI.instantiate(SDPA.Optimizer, with_bridge_type = Float64),
    )
    # `Variable.ZerosBridge` makes dual needed by some tests fail.
    MOI.Bridges.remove_bridge(
        model.optimizer,
        MathOptInterface.Bridges.Variable.ZerosBridge{Float64},
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            rtol = 1e-3,
            atol = 1e-3,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
                MOI.ObjectiveBound,
                MOI.SolverVersion,
            ],
        ),
        exclude = String[
            # Unable to bridge RotatedSecondOrderCone to PSD because the dimension is too small: got 2, expected >= 3.
            "test_conic_SecondOrderCone_INFEASIBLE",
            "test_constraint_PrimalStart_DualStart_SecondOrderCone",
            # Expression: MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
            #  Evaluated: MathOptInterface.INFEASIBLE_OR_UNBOUNDED == MathOptInterface.INFEASIBLE
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
            # Incorrect objective
            # See https://github.com/jump-dev/MathOptInterface.jl/issues/1759
            "test_unbounded_MIN_SENSE",
            "test_unbounded_MIN_SENSE_offset",
            "test_unbounded_MAX_SENSE",
            "test_unbounded_MAX_SENSE_offset",
            "test_infeasible_MAX_SENSE",
            "test_infeasible_MAX_SENSE_offset",
            "test_infeasible_MIN_SENSE",
            "test_infeasible_MIN_SENSE_offset",
            "test_infeasible_affine_MAX_SENSE",
            "test_infeasible_affine_MAX_SENSE_offset",
            "test_infeasible_affine_MIN_SENSE",
            "test_infeasible_affine_MIN_SENSE_offset",
            # TODO remove when PR merged
            # See https://github.com/jump-dev/MathOptInterface.jl/pull/1769
            "test_objective_ObjectiveFunction_blank",
            # FIXME investigate
            #  Expression: isapprox(MOI.get(model, MOI.ObjectiveValue()), T(2), config)
            #   Evaluated: isapprox(5.999999984012059, 2.0, ...
            "test_modification_delete_variables_in_a_batch",
            # FIXME investigate
            #  Expression: isapprox(MOI.get(model, MOI.ObjectiveValue()), objective_value, config)
            #   Evaluated: isapprox(-2.1881334077988868e-7, 5.0, ...
            "test_objective_qp_ObjectiveFunction_edge_case",
            # FIXME investigate
            #  Expression: isapprox(MOI.get(model, MOI.ObjectiveValue()), objective_value, config)
            #   Evaluated: isapprox(-2.1881334077988868e-7, 5.0, ...
            "test_objective_qp_ObjectiveFunction_zero_ofdiag",
            # FIXME investigate
            #  Expression: isapprox(MOI.get(model, MOI.ConstraintPrimal(), index), solution_value, config)
            #   Evaluated: isapprox(2.5058846553349667e-8, 1.0, ...
            "test_variable_solve_with_lowerbound",
        ],
    )
    return
end

end  # module

TestSDPA.runtests()


text:
function _JLine(p1, p2, linewidth, color)
    sethue(color)
    setline(linewidth)
    line(p1, p2, :stroke)
    return p2
end

"""
    1. JLine(pt1::Point, pt2::Point; kwargs...)
    2. JLine(pt2::Point; kwargs...)
        - `pt1` is set as the origin `O`

# Keywords for all
- `color` = "black"
- `linewidth` = 1

Draw a line between the points pt1 and pt2.
Returns the final point of the line
"""
JLine(pt1::Point, pt2::Point; linewidth = 1, color = "black") =
    (args...; color = color, linewidth = linewidth, pt1 = pt1, pt2 = pt2) ->
        _JLine(pt1, pt2, linewidth, color)

JLine(pt::Point; kwargs...) = JLine(O, pt; kwargs...)


text:
Dict(
     "T2" => "T2 - en",
     "T3" => "T3 - en",
     "T4" => "T4 - en",
     "T5" => "T5 - en",
    )


text:
include(joinpath(@__DIR__, "config.jl"))
include(joinpath(@__DIR__, "..", "subspace_util.jl"))

using OneClassActiveLearning, SVDD, JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

data_dirs = ["HeartDisease", "Stamps", "Pima", "Cardiotocography", "SpamBase", "Annthyroid", "PageBlocks"]
data_info = "04-subsvdd-smallsubspaces"

classify_precision = SVDD.OPT_PRECISION

### learning scenario ###
initial_pool_strategy = [("Pu", Dict())]
split_strategy = [("Sf", Dict())]

#### models ####
models = [Dict(:type => :SubSVDD,
               :param => Dict{Symbol, Any}(:weight_update_strategy => SVDD.FixedWeightStrategy(10.0, 0.01)))]

init_strategies = [SimpleSubspaceStrategy(SVDD.WangGammaStrategy(solver, range(0.1, stop=20, length=10), 1.0), FixedCStrategy(0.45), gamma_scope=Val(:Subspace))]

subspace_generators = [generate_randsub_from_data(10; min_ss_size=2, max_ss_size=2),
                       generate_randsub_from_data(20; min_ss_size=2, max_ss_size=2),
                       generate_randsub_from_data(10; min_ss_size=2, max_ss_size=4),
                       generate_randsub_from_data(20; min_ss_size=2, max_ss_size=4)]

#### oracle ####
oracle_param = Dict{Symbol, Any}(
    :type => PoolOracle,
    :param => Dict{Symbol, Any}()
)

#### query strategies ####
query_strategies = [
    Dict(:type => SubspaceQs{DecisionBoundaryPQs},
           :param => Dict{Symbol, Any}(:scale_fct => min_max_normalize,
                                       :combination_fct => + )),
    Dict(:type => SubspaceQs{HighConfidencePQs},
        :param => Dict{Symbol, Any}(:scale_fct => min_max_normalize,
                                      :combination_fct => + )),
    Dict(:type => SubspaceQs{RandomOutlierPQs},
         :param => Dict{Symbol, Any}(:scale_fct => min_max_normalize,
                                     :combination_fct => + ))]

num_al_iterations = 50
num_resamples_initial_pool = 1
num_resamples_subspaces = 3
data_output_root = joinpath(data_root, "output", "evaluation-part1")

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using Distributions
using Turing
using Plots
# using Statistics

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 4.64
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))

# %% 4.65, 6.66
d.weight_s = (d.weight .- mean(d.weight)) / std(d.weight)

f_parabola(weight_s, a, b1, b2) = a + b1 * weight_s + b2 * weight_s^2

@model function parabola(weight_s, heights)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    σ ~ Uniform(0, 50)
    μ = f_parabola.(weight_s, a, b1, b2)
    heights ~ MvNormal(μ, σ)
end

q4_5 = quap(parabola(d.weight_s, d.height), NelderMead())

# precis(m4_5)

# %% 4.67
weight_seq = range(-2.2, 2, length = 30)
post = DataFrame(rand(q4_5.distr, 1_000)', q4_5.params)
mu = f_parabola.(weight_seq, post.a', post.b1', post.b2')
sim = rand.(Normal.(mu, post.σ'))

mu = meanlowerupper(mu)
sim = meanlowerupper(sim)

# %% 4.68
scatter(d.weight_s, d.height, ms = 3, alpha = 0.7, legend = false)
plot!(weight_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(weight_seq, sim.lower, fillrange = sim.upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %% 4.69
f_cube(weight_s, a, b1, b2, b3) = a + b1 * weight_s + b2 * weight_s^2 + b3 * weight_s^3

@model function cube(weight_s, heights)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 10)
    b3 ~ Normal(0, 10)
    σ ~ Uniform(0, 50)
    μ = f_cube.(weight_s, a, b1, b2, b3)
    heights ~ MvNormal(μ, σ)
end

q4_6 = quap(cube(d.weight_s, d.height), NelderMead())

# %% 4.70, 4.71
weight_seq = range(-2.2, 2, length = 30)
post = DataFrame(rand(q4_6.distr, 1_000)', q4_6.params)
mu = f_cube.(weight_seq, post.a', post.b1', post.b2', post.b3') |> meanlowerupper
sim = rand.(Normal.(mu.raw, post.σ')) |> meanlowerupper

weight_seq_rescaled = weight_seq .* std(d.weight) .+ mean(d.weight)
scatter(d.weight, d.height, ms = 3, alpha = 0.7, legend = false)
plot!(weight_seq_rescaled, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(weight_seq_rescaled, sim.lower, fillrange = sim.upper, alpha = 0.3, la = 0.0, c = 2)

# %% 4.72
d = DataFrame!(CSV.File(datadir("exp_raw/cherry_blossoms.csv"), missingstring = "NA"))

# precis(d)

scatter(d.year, d.doy)

# %% 4.73
d2 = dropmissing(d, :doy)

num_knots = 15
knot_list = quantile(d2.year, range(0, 1, length = num_knots))

# %% 4.74, 4.75
using BSplines: BSplineBasis, basismatrix

basis = BSplineBasis(4, knot_list)
B = basismatrix(basis, d2.year)

plot(legend = false, xlabel = "year", ylabel = "basis value")
for y in eachcol(B)
    plot!(d2.year, y)
end
plot!()

# %% 4.76
@model function spline(D, B = B)
    α ~ Normal(100, 10)
    w ~ filldist(Normal(0, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = α .+ B * w
    D ~ MvNormal(μ, σ)
    return μ
end

q4_7 = quap(spline(d2.doy))

# %% 4.77
w_str = ["w[$i]" for i in 1:length(basis)]
post = DataFrame(rand(q4_7.distr, 1000)', ["α"; w_str; "σ"])

w = mean.(eachcol(post[:, w_str]))              # either
w = [mean(post[:, col]) for col in w_str]       # or

plot(legend = false, xlabel = "year", ylabel = "basis * weight")
for y in eachcol(B .* w')
    plot!(d2.year, y)
end
plot!()

# %% 4.78
mu = post.α' .+ B * Array(post[!, w_str])'
mu = meanlowerupper(mu)

scatter(d2.year, d2.doy, alpha = 0.3)
plot!(d2.year, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))

# %% 4.79
@model function spline(D, B = B)
    α ~ Normal(100, 10)
    w ~ filldist(Normal(0, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = [α + sum(Brow .* w) for Brow in eachrow(B)]
    D ~ MvNormal(μ, σ)
end

q4_7alt = quap(spline(d2.doy))


text:
using Ipopt, BenchmarkTools

# options
T = Float64
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    iterations=300)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=40,
    cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-4,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    R_inf=1.0e-8,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    opts=Dict(:max_iter=>10000),
    feasibility_tolerance=max_con_viol)

# iLQR
prob = copy(Problems.quadrotor_static)
U0 = deepcopy(controls(prob))
ilqr = iLQRSolver(prob, opts_ilqr)
initial_controls!(ilqr, U0)
solve!(ilqr)
ilqr.stats.iterations

@btime begin
    initial_controls!($ilqr, $U0)
    solve!($ilqr)
end
# 1.25 ms/iteration @ jarvis
# 0.76 ms/iteration @ bellman

# AL-iLQR
prob = copy(Problems.quadrotor_static)
U0 = deepcopy(controls(prob))
alilqr = AugmentedLagrangianSolver(prob)
initial_controls!(alilqr, U0)
solve!(alilqr)
alilqr.stats.cost
max_violation(alilqr)
cost(alilqr)

@btime begin
    initial_controls!($alilqr, $U0)
    solve!($alilqr)
end

# ALTRO
prob = copy(Problems.quadrotor_static)
altro = ALTROSolver(prob, opts_altro)
initial_controls!(altro, U0)
solve!(altro)
max_violation(altro)

@btime begin
    initial_controls!($altro, $U0)
    solve!($altro)
end

# Ipopt
prob_ipopt = copy(Problems.cartpole_static)
prob_ipopt = TO.change_integration(prob_ipopt, HermiteSimpson)
rollout!(prob_ipopt)
opts_ipopt.verbose = true
ipopt = DIRCOLSolver(prob_ipopt, opts_ipopt)
ipopt.optimizer.options
MOI.optimize!(ipopt.optimizer)
solve!(ipopt)
ipopt.stats
max_violation(ipopt)

@btime solve!($ipopt)


text:
module Coords

include("states.jl")
include("conversions.jl")
include("bodyfixed.jl")
include("time_series.jl")
include("events.jl")
include("trajectories.jl")

end


text:
module INP

export  InpStiffness

using ...TopOptProblems: Metadata, StiffnessTopOptProblem
using Ferrite
using ...Utilities: find_black_and_white, find_varind
import ...TopOptProblems: nnodespercell, getE, getν, getgeomorder, getdensity, getpressuredict, getcloaddict, getfacesets

include(joinpath("Parser", "Parser.jl"))
using .Parser

include("inpstiffness.jl")

end


text:
using JLD2, Polynomials
JLD2.@load "256K.jld2"

f = open("PERMX.txt","w") # do this once

write(f, "PERMX")
write(f, "\n")
ct = 0
for k = 1:256
    println("depth level = $(k)")
    for j = 1:256
        for i = 1:256
            global ct = ct + 1
            write(f, string(K[i,j,k]))
            write(f, " ")
            (ct%5 == 0) && write(f, "\n")
        end
    end
end

write(f, " /")
close(f)

f1 = open("PERMY.txt","w") # do this once

write(f1, "PERMY")
write(f1, "\n")
ct = 0
for k = 1:256
    println("depth level = $(k)")
    for j = 1:256
        for i = 1:256
            global ct = ct + 1
            write(f1, string(K[i,j,k]))
            write(f1, " ")
            (ct%5 == 0) && write(f1, "\n")
        end
    end
end

write(f1, " /")
close(f1)

f2 = open("PERMZ.txt","w") # do this once

write(f2, "PERMZ")
write(f2, "\n")
ct = 0
for k = 1:256
    println("depth level = $(k)")
    for j = 1:256
        for i = 1:256
            global ct = ct + 1
            write(f2, string(K[i,j,k]))
            write(f2, " ")
            (ct%5 == 0) && write(f2, "\n")
        end
    end
end

write(f2, " /")
close(f2)

f3 = open("PORO.txt","w") # do this once

write(f3, "PORO")
write(f3, "\n")
ct = 0
for k = 1:256
    println("depth level = $(k)")
    for j = 1:256
        for i = 1:256
            global ct = ct + 1
            p = Polynomial([-0.0314^2*K[i,j,k],2*0.0314^2*K[i,j,k],-0.0314^2*K[i,j,k],1.527^2])
            phi = minimum(real(roots(p)[findall(real(roots(p)).== roots(p))]))
            write(f3, string(phi))
            write(f3, " ")
            (ct%5 == 0) && write(f3, "\n")
        end
    end
end
write(f3, " /")
close(f3)

max_stars_repo_name: slimgroup/TwoPhaseSpeedTest

max_issues_repo_name: slimgroup/TwoPhaseSpeedTest

max_forks_repo_name: slimgroup/TwoPhaseSpeedTest


text:
#! /usr/bin/julia

# Rosetta Code, Hickerson series of almost integers

function makehickerson{T<:Real}(x::T)
    n = 0
    h = one(T)/2x
    function hickerson()
        n += 1
        h *= n/x
    end
end

function reporthickerson{T<:Real,U<:Integer}(a::T, nmax::U)
    h = makehickerson(a)
    hgm = makehickerson(prevfloat(a))
    hgp = makehickerson(nextfloat(a))

    println()
    print("Performing calculations using ", typeof(a))
    println(", which has ", precision(a), "-bit precision.")
    for i in 1:nmax
        x = h()
        xm = hgm()
        xp = hgp()
        y = ifloor(10x)
        ym = ifloor(10xm)
        yp = ifloor(10xp)
        println()
        println("Hickerson series result for n = ", i)
        println(@sprintf("    ->  %25.4f ", xm))
        println(@sprintf("    0>  %25.4f ", x))
        println(@sprintf("    +>  %25.4f ", xp))
        isprecok =
        isint =
        if ym == y == yp
            print("The precision is adequate, ")
            if  digits(y)[1] in [0, 9]
                println("and the result is an almost integer.")
            else
                println("but the result is not an almost integer.")
            end
        else
            println("The precision is inadequate for a definite result.")
        end
    end
end

a = log(big(2.0))
reporthickerson(a, 17)

a = log(2.0)
reporthickerson(a, 17)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# Newton solver options
@with_kw mutable struct NewtonOptions{T}
    r_tol::T = 1.0e-5            # primal dual residual tolerance
	max_iter::Int = 10           # maximum number of primal dual iter
    max_time::T = 10000.0        # maximum time spent in the Newton solver
    β_init::T = 1.0e-5           # initial dual regularization
    live_plotting::Bool = false  # visualize the trajectory during the solve
    threads::Bool=false
    verbose::Bool = false
    solver::Symbol = :lu_solver  # lu_sparse_solver
end

abstract type NewtonIndices end
abstract type NewtonResidual end
abstract type NewtonJacobian end

mutable struct Newton{T,nq,nu,nw,nc,nb,nz,nθ,nν,NJ,NR,NI,O,LS,NV}
    jac::NJ                         # NewtonJacobian
    res::NR                        # residual
    res_cand::NR                     # candidate residual
    Δ::NR                          # step direction in the Newton solve, it contains: q2-qH+1, u1-uH, γ1-γH, b1-bH, λd1-λdH
    ν::Vector{NV}          # implicit dynamics lagrange multiplier
    ν_cand::Vector{NV}         # candidate implicit dynamics lagrange multiplier
    traj::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}                           # optimized trajectory
    traj_cand::ContactTraj{T,nq,nu,nw,nc,nb,nz,nθ}      # trial trajectory used in line search
    Δq::Vector{Vector{T}}#::Vector{SizedArray{Tuple{nq},T,1,1}}         # difference between the traj and ref_traj
    Δu::Vector{Vector{T}}         # difference between the traj and ref_traj
    Δγ::Vector{Vector{T}}         # difference between the traj and ref_traj
    Δb::Vector{Vector{T}}         # difference between the traj and ref_traj
    ind::NI                                 # indices of a one-time-step block
    obj::O                              # obj function
    solver::LS
    β::T
    opts::NewtonOptions{T}                          # Newton solver options
end

function Newton(s::Simulation{T}, H::Int, h::T,
    traj::ContactTraj{T}, im_traj::ImplicitTrajectory{T};
    obj::Objective = TrackingObjective(s.model, s.env, H),
    opts::NewtonOptions{T} = NewtonOptions(), κ::T=im_traj.ip[1].κ[1]) where T

    model = s.model
    env = s.env

    mode = im_traj.mode

    ind = NewtonIndices(model, env, H, mode = mode)

    nq = model.nq
    nu = model.nu
    nw = model.nw
    nc = model.nc
    nb = nc * friction_dim(env)
    nz = num_var(model, env)
    nθ = num_data(model)
    nd = ind.nd

    jac = NewtonJacobian(model, env, H, mode = mode)

    # precompute Jacobian for pre-factorization
    window = collect(1:(H + 2))
    implicit_dynamics!(im_traj, traj, window=window, threads=opts.threads) #@@@
    jacobian!(jac, im_traj, obj, H, opts.β_init, window)

    res = NewtonResidual(model, env, H, mode = mode)
    res_cand = NewtonResidual(model, env, H, mode = mode)

    Δ = NewtonResidual(model, env, H, mode = mode)

    ν = [zeros(SizedVector{ind.nd,T}) for t = 1:H]
    ν_cand = deepcopy(ν)

    traj = contact_trajectory(model, env, H, h, κ=κ)
    traj_cand = contact_trajectory(model, env, H, h, κ=κ)

    Δq  = [zeros(nq) for t = 1:H]

    Δu  = [zeros(nu) for t = 1:H]
    Δγ  = [zeros(nc) for t = 1:H]
    Δb  = [zeros(nb) for t = 1:H]

    # regularization
    β = copy(opts.β_init)

    # linear solver
    solver = eval(opts.solver)(jac.R)

    return Newton{T,nq,nu,nw,nc,nb,nz,nθ,nd,typeof.([jac,res,ind,obj,solver,ν[1]])...}(
        jac, res, res_cand, Δ, ν, ν_cand, traj, traj_cand,
        Δq, Δu, Δγ, Δb, ind, obj, solver, β, opts)
end

function delta!(Δx::SizedArray{Tuple{nx},T,1,1}, x, x_ref) where {nx,T}
    Δx .= x
    Δx .-= x_ref
    return nothing
end

function delta!(Δx::Vector{T}, x, x_ref) where T
    Δx .= x
    Δx .-= x_ref
    return nothing
end

function copy_traj!(traj::ContactTraj, traj_cand::ContactTraj, H::Int)
    Ht = traj.H
    Hs = traj_cand.H # MAYBE BREAKING TEST

    @assert Hs >= H
    @assert Ht >= H

    traj.κ .= traj_cand.κ

    for t in eachindex(1:H + 2)
        traj.q[t] .= traj_cand.q[t]
    end

    for t in eachindex(1:H)
        traj.u[t] .= traj_cand.u[t]
        traj.w[t] .= traj_cand.w[t]
        traj.γ[t] .= traj_cand.γ[t]
        traj.b[t] .= traj_cand.b[t]
        traj.z[t] .= traj_cand.z[t]
        traj.θ[t] .= traj_cand.θ[t]
    end

    return nothing
end

function reset!(core::Newton, ref_traj::ContactTraj,
    q0::Vector{T}, q1::Vector{T};
    window=collect(1:(core.traj.H + 2)),
    warm_start::Bool = false) where T

    # H = ref_traj.H
    H_mpc = core.traj.H
    opts = core.opts

    # Reset β value
    core.β = opts.β_init

    if !warm_start
		# Reset duals
		for t = 1:H_mpc
			fill!(core.ν[t], 0.0)
			fill!(core.ν_cand[t], 0.0)
		end

        # TODO: not sure this is correct
        # Set up trajectory
        copy_traj!(core.traj, ref_traj, core.traj.H)
	end

    core.traj.q[1] .= q0
    core.traj.q[2] .= q1

    update_θ!(core.traj, 1)
    update_θ!(core.traj, 2)

    # initialized residual Jacobian
    initialize_jacobian!(core.jac, core.obj, core.traj.H)

    # Set up traj cand
    copy_traj!(core.traj_cand, core.traj, core.traj.H)

	return nothing
end

function newton_solve!(
    core::Newton{T},
    s::Simulation{T},
    q0::Vector{T},
    q1::Vector{T},
    window::Vector{Int},
    im_traj::ImplicitTrajectory{T},
    ref_traj::ContactTraj{T};
    warm_start::Bool=false) where T

    elapsed_time = 0.0

    elapsed_time += @elapsed begin
        # reset solver
        reset!(core, ref_traj, q0, q1,
            window=window,
            warm_start=warm_start)
    end
    elapsed_time >= core.opts.max_time && (return nothing)

    elapsed_time += @elapsed begin
        # Compute implicit dynamics about traj
        implicit_dynamics!(im_traj, core.traj, window=window, threads=core.opts.threads)
    end
    elapsed_time >= core.opts.max_time && (return nothing)

    elapsed_time += @elapsed begin
        # Compute residual
        residual!(core.res, core, core.ν, im_traj, core.traj, ref_traj, window)
        r_norm = norm(core.res.r, 1)
    end
    elapsed_time >= core.opts.max_time && (return nothing)

    for l = 1:core.opts.max_iter
		# elapsed_time >= core.opts.max_time && break
		# elapsed_time += @elapsed begin
        # check convergence
        r_norm / length(core.res.r) < core.opts.r_tol && break

        elapsed_time += @elapsed begin
            # Compute NewtonJacobian
            jacobian!(core.jac, im_traj, core.obj, core.traj.H, core.β, window, update_hessian=true)
        end
        elapsed_time >= core.opts.max_time && (return nothing)
        # jacobian!(core.jac, im_traj, core.obj, core.traj.H, core.β, update_hessian=false)
        # @show diag(core.jac.R[1:6,1:6])

        elapsed_time += @elapsed begin
            # Compute Search Direction
            linear_solve!(core.solver, core.Δ.r, core.jac.R, core.res.r)
        end
        elapsed_time >= core.opts.max_time && (return nothing)

        # line search the step direction
        α = 1.0
        iter = 0

        elapsed_time += @elapsed begin
            # candidate step
            update_traj!(core.traj_cand, core.traj, core.ν_cand, core.ν, core.Δ, α)
        end
        elapsed_time >= core.opts.max_time && (return nothing)

        elapsed_time += @elapsed begin
            # Compute implicit dynamics for candidate
            implicit_dynamics!(im_traj, core.traj_cand, window=window, threads=core.opts.threads)
        end
        elapsed_time >= core.opts.max_time && (return nothing)

        elapsed_time += @elapsed begin
            # Compute residual for candidate
            residual!(core.res_cand, core, core.ν_cand, im_traj, core.traj_cand, ref_traj, window)
            r_cand_norm = norm(core.res_cand.r, 1)
        end
        elapsed_time >= core.opts.max_time && (return nothing)

        while r_cand_norm^2.0 >= (1.0 - 0.001 * α) * r_norm^2.0
            α = 0.5 * α

            iter += 1
            if iter > 6
                break
            end

            elapsed_time += @elapsed begin
                update_traj!(core.traj_cand, core.traj, core.ν_cand, core.ν, core.Δ, α)
            end
            elapsed_time >= core.opts.max_time && (return nothing)

            elapsed_time += @elapsed begin
                # Compute implicit dynamics about trial_traj
                implicit_dynamics!(im_traj, core.traj_cand, window=window, threads=core.opts.threads)
            end
            elapsed_time >= core.opts.max_time && (return nothing)

            elapsed_time += @elapsed begin
                residual!(core.res_cand, core, core.ν_cand, im_traj, core.traj_cand, ref_traj, window)
                r_cand_norm = norm(core.res_cand.r, 1)
            end
            elapsed_time >= core.opts.max_time && (return nothing)
        end

        elapsed_time += @elapsed begin
            # update
            update_traj!(core.traj, core.traj, core.ν, core.ν, core.Δ, α)
            core.res.r .= core.res_cand.r
            r_norm = r_cand_norm
        end
        elapsed_time >= core.opts.max_time && (return nothing)

        # regularization update
        iter > 6 ? (core.β = min(core.β * 1.3, 1.0e2)) : (core.β = max(1.0e1, core.β / 1.3))

        # print
        core.opts.verbose && print_status(core, im_traj, elapsed_time, α, l)
		# end
    end

    return nothing
end

function print_status(core::Newton, im_traj, elapsed_time, α, l)
     # print status
     println(
     "     t: ", scn(elapsed_time, digits=0),
     "     r̄: ", scn(norm(core.res_cand.r, 1) / length(core.res_cand.r), digits=0),
     "     r: ", scn(norm(core.res.r, 1) / length(core.res.r), digits=0),
     "     Δ: ", scn(norm(core.Δ.r, 1) / length(core.Δ.r), digits=0),
     "     α: ", -Int(round(log(α))),
     "     l: ", l,
     "     κ: ", scn(im_traj.ip[1].κ[1], digits = 0)
     )
end


text:
export schedule_on

import Base: show

# TODO : Untested and undocumented

schedule_on(scheduler::H) where { H <: AbstractScheduler } = ScheduleOnOperator{H}(scheduler)

struct ScheduleOnOperator{ H <: AbstractScheduler } <: InferableOperator
    scheduler :: H
end

operator_right(operator::ScheduleOnOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::ScheduleOnOperator{H}, source::S) where { L, H, S }
    return scheduled(source, operator.scheduler)
end

Base.show(io::IO, ::ScheduleOnOperator{H}) where H = print(io, "ScheduleOnOperator($H)")


text:
"""
Out Block
"""

export OutBlock

mutable struct OutBlock <: AbstractBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}

    function OutBlock()
        @createblock new(Vector{InPort}(), Vector{OutPort}()) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::OutBlock)
    println(io, "OutBlock()")
end

"""
to expr
"""

function _toexpr(blk::OutBlock)
    _toexpr(blk.inport[1])
end


text:
using Pkg

"""
	rootdir()

Root directory path where objects will be stored.
"""
rootdir() = joinpath(first(DEPOT_PATH), "store")

"""
	itemdir(obj)

Directory path where `obj` will be stored.
"""
itemdir(obj) = joinpath(rootdir(), string(hash(obj), base=16))

"""
	issaved(obj)

Returns true if `obj` is locally stored, false otherwise.
"""
issaved(obj) = isdir(itemdir(obj))

"""
	save(obj)

Stores `obj` locally if it has not been stored already.
Returns path to directory containing `obj`.
"""
function save(obj)
	if !issaved(obj)
		@info "Storing $obj into $(itemdir(obj))"
		mkpath(itemdir(obj))
		save(obj, itemdir(obj))
		#verify(obj)
	end
	return itemdir(obj)
end

"""
	load(obj)

Reads and returns the specified object. Stores it locally if it hasn't been
stored already.
"""
function load(obj)
	issaved(obj) || save(obj)
	load(obj, itemdir(obj))
end

"""
	verify(obj)

Compares the tree_hash of the stored copy of `obj` to the expected tree_hash,
and errors if there is a difference.
"""
function verify(obj)
	issaved(obj) || error("$obj not stored.")
	if Pkg.GitTools.tree_hash(itemdir(obj)) ≠ expected_tree_hash(obj)
		error("Corrupted store for $obj.")
	end
end

"""
	expected_tree_hash(obj)

Expected tree hash for the saved contents of `obj`. Overload this for specific
item types. Otherwise simply returns the hash of the contents of the store
(so no verification happens).
"""
function expected_tree_hash(obj)
	issaved(obj) || error("$obj not stored.")
	Pkg.GitTools.tree_hash(itemdir(obj))
end

"""
	clean(; force=false)

Removes all stored data directories.
"""
function clean(; force::Bool=false)
	if force
		rm(rootdir(), recursive=true, force=true)
	else
		print("Clean store. Type 'yes' to confirm: ")
		if readline() == "yes"
			clean(; force=true)
		end
	end
end

"""
	clean(obj)

Removes the stored data associated with `obj`.
"""
clean(obj) = issaved(obj) && rm(itemdir(obj), recursive=true)


text:
mutable struct NormalDistributionLikelihood <: MyXMLElement
    el::XMLOrNothing
    param::MyXMLElement
    mean::Float64
    sd::Float64
    id::String

    function NormalDistributionLikelihood(param::MyXMLElement, μ::Float64,
                                          σ::Float64, id::String)
        return new(nothing, param, μ, σ, id)
    end
end

function NormalDistributionLikelihood(param::MyXMLElement, id::String)
    return NormalDistributionLikelihood(param, 0.0, 1.0, id)
end

function make_xml(ndl::NormalDistributionLikelihood)
    el = new_element(bn.DISTRIBUTION_LIKELIHOOD)
    set_attribute(el, bn.ID, ndl.id)

    d_el = new_child(el, bn.DATA)
    make_xml(ndl.param)
    add_ref_el(d_el, ndl.param.el)
    dist_el = new_child(el, bn.DISTRIBUTION)

    ndm_el = new_child(dist_el, bn.NORMAL_DISTRIBUTION_MODEL)

    m_el = new_child(ndm_el, bn.MEAN)
    add_parameter(m_el, value = [ndl.mean])
    std_el = new_child(ndm_el, bn.STDEV)
    add_parameter(std_el, value =[ndl.sd], lower=0.0)

    ndl.el = el

    return el
end

function get_normal_prior(ndl::NormalDistributionLikelihood)
    make_xml(ndl)
    return ndl.el
end

function get_priors(ndl::NormalDistributionLikelihood)
    make_xml(ndl)
    return ndl.el
end

function get_loggables(::NormalDistributionLikelihood)
    return MyXMLElement[]
end

function make_loadings_gradient(loadings::MyXMLElement,
                                ndl::NormalDistributionLikelihood)
    el = new_element(bn.GRADIENT)
    make_xml(ndl)
    make_xml(loadings)

    add_ref_el(el, ndl.el)
    add_ref_el(el, loadings.el)

    return el
end

function get_loadings_prior(ndl::NormalDistributionLikelihood)
    return ndl
end


text:
module QuSAT

###################################################################
#
# This file generates SAT's Instance in classical method for tests
# in quantum computing.
# Ref: arxiv:quant-ph/0007071v1
#
###################################################################

import Base: show,call
export AbstractBits, AbstractClause, Bits,
        Clause, ECClause, TruthTable, save,
        generate, Instance, AssignNum, readins

abstract AbstractBits
abstract AbstractClause{N}

immutable Bits <: AbstractBits
    value::UInt32

    Bits(n::Integer) = new(n)
end

function call(b::Bits,n::Integer)
    @assert n>0
    return (b.value>>(n-1))&1
end

# function convert(Bits,Integer)
#     return Bits()

"""
TruthTable Example:

|bits value | truth value|
|-----------|------------|
| 000       | 0          |
| 001       | 1          |
| 010       | 0          |
| 011       | 1          |
| 100       | 0          |
| 101       | 1          |
| 110       | 0          |
| 111       | 0          |

TruthTable(0b00101010)
"""
immutable TruthTable
    value::UInt32
end

function call(truthtable::TruthTable,index::Integer)
    return (truthtable.value>>index)&1
end

immutable Clause{N} <: AbstractClause{N}
    table::TruthTable
    ids::AbstractVector{Integer}

    function Clause(table::TruthTable,ids::AbstractVector{Integer})
        @assert length(ids)==N
        sort!(ids)
        new(table,ids)
    end
end

Clause(table::TruthTable,ids::Integer...) = Clause{length(ids)}(table,[ids...])
Clause(num::Integer,table::TruthTable,ids::Integer...) = Clause{num}(table,[ids...])
Clause(num::Integer,table::Integer,ids::Integer...) = Clause{num}(TruthTable(table),[ids...])

function call{N}(clause::Clause{N},assign::Integer)
    return clause.table(assign)
end

################################################################################
#
# Clauses in exact cover problem
#
################################################################################

type ECClause{N} <: AbstractClause{N}
    ids::Vector{Int}

    function ECClause(ids::Vector{Int})
        @assert length(ids) == N
        sort!(ids)
        new(ids)
    end
end

ECClause(ids::Vector{Int}) = ECClause{length(ids)}(ids)
ECClause(ids::Integer...) = ECClause{length(ids)}([ids...])

function call{N}(c::ECClause{N},assign::Integer)
    res = 0
    for i = 1:N
        res += assign&1
        assign = assign>>1
    end
    return Int(res==1)
end

function save{N}(io::IO,clause::ECClause{N})
    write(io,"$(clause.ids[1])")
    for id in clause.ids[2:end]
        write(io,"\t")
        write(io,"$(id)")
    end
    write(io,"\n")
end

"""
Instance is a collection of M clauses:

~~~TeX
C_1 Λ C_2 Λ... C_M
~~~

construct an instance:

- `num` is the number of bits
- `clauses` is the collection of clause
"""
type Instance{M,N}
    c::AbstractVector{AbstractClause{N}}
end

Instance{N}(num::Integer,clauses::AbstractVector{ECClause{N}}) = Instance{num,N}(clauses)
Instance{N}(num::Integer,clause::ECClause{N},clauses::ECClause{N}...) = Instance(num,[clause,clauses...])
Instance{N}(num::Integer,clauses::AbstractVector{Clause{N}}) = Instance{num,N}(clauses)
Instance{N}(num::Integer,clause::Clause{N},clauses::Clause{N}...) = Instance(num,[clause,clauses...])

function call{M,N}(clauses::Instance{M,N},assign::Bits)
    res = 1
    for clause in clauses.c
        assignment = 0
        digit = 0
        for id in clause.ids
            assignment += assign(id)<<digit
            digit += 1
        end
        res &= clause(assignment)
    end
    return res
end

function show{M,N}(io::IO,ins::Instance{M,N})
    for clause in ins.c[1:end-1]
        print(io,clause.ids)
        print(io,"Λ")
    end
    print(io,ins.c[end].ids)
    print(io,"\n")
end

show{M,N}(ins::Instance{M,N}) = show(STDOUT,ins)

function save{M,N}(io::IO,ins::Instance{M,N})
    for clause in ins.c
        save(io,clause)
    end
end

function readins(io::IO)
    data = convert(Array{Int,2},readdlm(io))
    clauses = ECClause[ECClause(vec(data[i,:])) for i=1:size(data)[1]]
    @show typeof(clauses[1])
    return Instance(size(data)[2],clauses)
end

function AssignNum{M,N}(instance::Instance{M,N})
    count = 0
    ans = Int[]
    for i = 0:2^M-1
        if instance(Bits(i)) == 1
            count+=1
            push!(ans,i)
        end
    end
    return count,ans
end

import Base: push!

function push!{M,N}(instance::Instance{M,N},clause::ECClause{N})
    push!(instance.c,clause)
end

################################################################################
#
# generate instance with only one assignment satisfied
#
################################################################################

function engine(n::Integer,N::Integer,maxtry=100)
    list = collect(1:n)
    ids = (list|>shuffle)[1:N]

    clauses = [ECClause(ids)]
    ins = Instance(n,clauses)

    for i = 1:maxtry
        pre_ids = ids
        ids = (list|>shuffle)[1:N]
        if sort!(ids) != sort!(pre_ids)
            push!(ins,ECClause(ids))
        end

        asign = AssignNum(ins)
        if asign[1]==1
            return ins,asign[2]
        end
    end

    return nothing
end

function generate(n::Integer,maxiter=1000;maxbits=16,N=3)
    for i=1:maxiter
        t = engine(n,N)
        t === nothing || return t
    end
    warn("may not have an assignment\n")
    return nothing,nothing
end

end # module


text:
"""
    Burgers{T,Dim}

Burgers' equation
``
  \\partial_t u + \\partial_x \\frac{u^2}{2} = 0
``
in `Dim` space dimensions using `T` as scalar type.
"""
struct Burgers{T,Dim} <: ScalarBalanceLaw{T,Dim} end

function Burgers(T=Float64, Dim=1)
  Burgers{T,Dim}()
end

function Base.show(io::IO, model::Burgers{T,Dim}) where {T,Dim}
  print(io, "Burgers' equation {T=", T, ", Dim=", Dim, "}",
            " with flux f(u) = u^2 / 2 * (1,...,1)/sqrt(Dim)")
end

"""
    IntegralQuantitiesBurgers{T<:Real}

Some integrated quantities of interest for Burgers' equation. Can be used in
callbacks of DifferentialEquations.jl.
"""
struct IntegralQuantitiesBurgers{T<:Real} <: FieldVector{2,T}
    mass::T
    energy::T
end

function IntegralQuantitiesBurgers(u, model::ScalarBalanceLaw)
    IntegralQuantitiesBurgers(u, u^2)
end

"""
    flux(u, model::Burgers)

Compute the flux of `u` for `model`.
"""
@inline flux(u, model::Burgers) = u^2/2

"""
    speed(u::Real, model::Burgers)

Compute the speed f'(`u`) for `model`.
"""
@inline speed(u, model::Burgers) = u

@inline flux(u, model::Burgers{T,2}, direction) where {T} = u^2 / (2*sqrt(2))

################################################################################

"""
    (::GodunovFlux)(uₗ::T, uᵣ::T, model::Burgers{T,1}) where {T}

Compute Godunov's flux between `uₗ` and `uᵣ` for `model`.
"""
function (::GodunovFlux)(uₗ::T, uᵣ::T, model::Burgers{T,1}) where {T}
  if uₗ < uᵣ
    if uₗ < 0 && 0 < uᵣ
      zero(T)
    else
      min(flux(uₗ, model), flux(uᵣ, model))
    end
  else
    max(flux(uₗ, model), flux(uᵣ, model))
  end
end

function (::GodunovFlux)(uₗ::T, uᵣ::T, model::Burgers{T,2}, direction) where {T}
  if uₗ < uᵣ
    if uₗ < 0 && 0 < uᵣ
      zero(T)
    else
      min(flux(uₗ, model, direction), flux(uᵣ, model, direction))
    end
  else
    max(flux(uₗ, model, direction), flux(uᵣ, model, direction))
  end
end

function (flux::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::Burgers{T,1}) where T<:Real
    (uₗ^2 + uₗ*uᵣ + uᵣ^2) / 6
end

function (flux::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::Burgers{T,2}, direction) where T<:Real
    (uₗ^2 + uₗ*uᵣ + uᵣ^2) / (6*sqrt(6))
end

@inline function (flux::L2L4ConservativeFlux)(uₗ::T, uᵣ::T, model::Burgers{T,1}) where T<:Real
    ( 18 * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{5}()) + 5 * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{3}()) ) /
        ( 60 * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{3}()) + 30 )
end

@generated function (flux::L2L2sConservativeFlux{s})(uₗ, uᵣ, model::Burgers{T,1}) where {s,T<:Real}
    ex = :(
        ( $(6s^2-3s) * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{$(2s+1)}())
        + $(2s+1) * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{3}()) ) /
            ( $(12s^2+6s) * jump_pow_u_r_over_jump_u(uₗ, uᵣ, Val{$(2s-1)}()) + $(12s+6) )
    )
    return :(Base.@_inline_meta; $ex)
end

################################################################################

@inline function add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law::Burgers,
                                                    Nx, basis, boundaries_included::Val{false},
                                                    jacx, parallel)
    Pp1 = length(basis.nodes)
    Rl = basis.interp_left
    Rr = basis.interp_right
    one_3 = 1 / 3
    one_6 = 1 / 6

    # add numerical fluxes
    @inbounds for ix in Base.OneTo(Nx)
        Rul = zero(eltype(u))
        Rur = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            Rul += Rl[nx]*u[nx,ix]
            Rur += Rr[nx]*u[nx,ix]
        end
        Ru2l = zero(eltype(u))
        Ru2r = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            tmp = u[nx,ix]^2
            Ru2l += Rl[nx]*tmp
            Ru2r += Rr[nx]*tmp
        end

        #@. du[:,ix] += ((fluxes[ix] - one_3 * Ru2l - one_6 * Rul^2) * Rl
        #                - (fluxes[ix+1] - one_3 * Ru2r - one_6 * Rur^2) * Rr
        #                ) * jacx / basis.weights
        for nx in Base.OneTo(Pp1)
            du[nx,ix] += ((fluxes[ix] - one_3 * Ru2l - one_6 * Rul^2) * Rl[nx]
                            - (fluxes[ix+1] - one_3 * Ru2r - one_6 * Rur^2) * Rr[nx]
                            ) * jacx / basis.weights[nx]
        end
    end

    nothing
end

################################################################################

"""
    BurgersRiemannSolution{T,T1}

The solution of a Riemann problem `prob` for Burgers' equation.
"""
struct BurgersRiemannSolution{T,T1} <: ScalarRiemannSolution
  prob::RiemannProblem{Burgers{T,1},T,T1}
  σ⁻::T
  σ⁺::T
end

"""
    minmax_speeds(sol::BurgersRiemannSolution)

Return the minimal and maximal speeds `σ⁻, σ⁺` that appear in the solution `sol`.
"""
function minmax_speeds(sol::BurgersRiemannSolution)
  sol.σ⁻, sol.σ⁺
end

"""
    (sol::BurgersRiemannSolution)(ξ::Real)

Evaluate the solution `sol` at the value `ξ` of the self-similarity variable
`ξ = (x - x₀) / (t - t₀)`.
"""
function (sol::BurgersRiemannSolution)(ξ::Real)
  @unpack σ⁻, σ⁺ = sol
  @unpack uₗ, uᵣ = sol.prob

  if ξ < σ⁻
    uₗ
  elseif ξ < σ⁺
    uₗ + (ξ-σ⁻)/(σ⁺-σ⁻) * (uᵣ-uₗ)
  else
    uᵣ
  end
end

"""
    (sol::BurgersRiemannSolution)(t::Real, x::Real)

Evaluate the solution `sol` at the time and space coordinates `t` and `x`.
"""
function (sol::BurgersRiemannSolution)(t::Real, x::Real)
  @unpack x₀, t₀ = sol.prob

  sol((x-x₀)/(t-t₀))
end

"""
    solve(prob::RiemannProblem{Burgers{T,1},T,T1}) where {T,T1}

Compute the solution of the Riemann prolem `prob`.
"""
function solve(prob::RiemannProblem{Burgers{T,1},T,T1}) where {T,T1}
  @unpack uₗ, uᵣ = prob
  if uₗ > uᵣ
    σ⁻ = σ⁺ = (uₗ + uᵣ) / 2
  else
    σ⁻ = uₗ
    σ⁺ = uᵣ
  end
  BurgersRiemannSolution(prob, σ⁻, σ⁺)
end


text:
using KernelMachines
using Test
using FiniteDiff: finite_difference_gradient
using Zygote: gradient

@testset "utils" begin
    s = rand(10, 3)
    slices = KernelMachines.split_matrix(s, (2, 3, 5))
    @test length(slices) == 3
    @test slices[1] == s[1:2, :]
    @test slices[2] == s[3:5, :]
    @test slices[3] == s[6:10, :]
end

@testset "kernelmachine" begin
    dims = (2, 3, 2)
    data = rand(5, 10)
    input = rand(5, 50)
    dm = KernelMachine(data; dims=dims)
    g_auto = gradient(input) do input
        r, n = dm(input)
        return sum(r) + n
    end
    g_num = finite_difference_gradient(input) do input
        r, n = dm(input)
        return sum(r) + n
    end
    @test isapprox(first(g_auto), g_num)
end


text:
"""
Create table of neighbor group tallys (how many neighbors are in each group?)

Inputs:
    adj: sparse adjacency matrix (NxN)
    groups: integer list of each node's group (Nx1)

Outputs:
    counts: MxN array of count of neighbors per each group
"""
function count_neighbor_groups(adj, groups)
    assert(size(adj,1) == length(groups))
    n = size(adj,1)
    counts = zeros(Int64,4,n)
    rows = rowvals(adj)
    for i in 1:n
        for j in nzrange(adj, i)
            row = rows[j]
            counts[groups[i],row] += 1
        end
    end
    return counts
end

"""
Count total degree (stubs) of each group

Inputs:
    ms: dict of edge count between groups

Outputs:
    ks: dict of total degree (stub) count for group
"""
function compile_degree_dict(ms)
    ks = Dict()
    for (a,b) in keys(ms)
        if !(a in keys(ks))
            ks[a] = 0
        end
        ks[a] += ms[a,b]
        if !(b in keys(ks))
            ks[b] = 0
        end
        ks[b] += ms[a,b]
    end
    for i in keys(ks)
        ks[i] = round(Int64, ks[i]/2)
    end
    return ks
end

"""
Compile dicts that count nodes within groups & edges between groups

Inputs:
    counts: MxN array of count of neighbors per each group
        (see count_neighbor_groups)
    groups: Nx1 group assignment of each node
    pairs: list of tuples, each tuple specifying a pair of groups
        Dicts will only be compiled for the group pairs provided. This pairs
        list will ensure that only these provided pairs are used in later
        computations.

Outputs:
    ms: dict of edge count between groups
    ns: dict of node count within groups
    ks: dict of total degree (stub) count for group
"""
function compile_count_dicts(counts, groups, pairs)
    ns = Dict()
    ms = Dict()
    for (a,b) in pairs
        ms[a,b] = sum(counts[b, groups.==a])
        ns[a] = sum(groups.==a)
        ns[b] = sum(groups.==b)
    end
    return ns, ms
end

"""
Calculate log-likelihood after swapping the group of one node
    The updates should be O(1) to prevent from total recalculation

Inputs:
    ns: count dict of current label pair connections
    count_deltas: Mx1 array of swapped node neighbor group counts
    old_group: old group of swapped node
    new_group: new group of swapped node

Outputs:
    new_ns: updated ns
    new_ms: updated ms
"""
function adjust_counts(ns, ms, count_deltas, old_group, new_group)
    new_ns = copy(ns)
    new_ns[old_group] = ns[old_group] - 1
    new_ns[new_group] = ns[new_group] + 1
    new_ms = Dict()
    for (a,b) in keys(ms)
        count_adjust = count_deltas[b]
        if a == new_group
            new_ms[a,b] = ms[a,b] + count_deltas[b]
        elseif a == old_group
            new_ms[a,b] = ms[a,b] - count_deltas[b]
        end
    end
    return new_ns, new_ms
end

"""
Calculate Bernoulli distribution log-likelihood based on count dicts

Inputs:
    ns: dict of node count within groups OR stub count for group (ks)
    ms: dict of edge count between groups

Outputs:
    ll: log-likelihood function value
"""
function calculate_ll(ns::Dict, ms::Dict, method="bernoulli")
    ll = 0
    for (a,b) in keys(ms)
        na = ns[a]
        nb = ns[b]
        m = ms[a,b]
        if na*nb != 0
            if method == "bernoulli"
                ll += m*log10(m/(na*nb)) + (na*nb-m)*log10(1-m/(na*nb))
            elseif method == "poisson"
                ll += m*log10(m/(na*nb))
            end
        end
    end
    return ll
end

"""
Maximize params

Inputs:
    ms: dict of edge count between groups
    ns: dict of node count within groups
    pr: MxM array of connection probabilities between groups

Outputs:
    updated pr
"""
function maximize_parameters!(ns, ms, pr)
    for (a,b) in keys(ms)
        pr[a,b] = ms[a,b] / (ns[a]*ns[b])
    end
end

"""
Get opposite group
"""
function swap_group(old_group)
    if old_group == 1
        return 2
    elseif old_group == 2
        return 1
    elseif old_group == 3
        return 4
    elseif old_group == 4
        return 3
    end
end

"""
Make dummy data to see how well the groups can be recovered
"""
function make_dummy_data(n=40)
    adj = spzeros(n,n)
    # 3/4 pre, 1/4 post; 1/2 excitatory, 1/2 inhibitory
    groups = ones(Int64, n)
    m = round(Int64, 3*n/8)
    groups[1:end-m] = ones(Int64, n-m)*2
    m = round(Int64, 6*n/8)
    groups[1:end-m] = ones(Int64, n-m)*3
    m = round(Int64, 7*n/8)
    groups[1:end-m] = ones(Int64, n-m)*4
    pairs = [(i,j) for i in [1,2] for j in [3,4]]
    # E_pre,I_pre,E_post,I_post x E_pre,I_pre,E_post,I_post
    pr = [0 0 0.15 0.12;
          0 0 0.02 0.02;
          0 0 0 0;
          0 0 0 0]
    for i in 1:n
        for j in 1:n
            ci, cj = groups[i], groups[j]
            p = pr[ci,cj]
            adj[i,j] = rand() < p
        end
    end
    return adj, groups, pr, pairs
end

function run_dummy_test()
    adj, groups, pr, pairs = make_dummy_data(1000)
    method = "poisson"
    groups_guess = copy(groups[r])
    groups_guess[groups_guess .<= 2] = rand([1,2], sum(groups .<= 2))
    groups_guess[groups_guess .>= 3] = rand([3,4], sum(groups .>= 3))
    k = size(adj,1)
    groups_changed = true
    n_iter = 3000
    iter = 0
    counts = count_neighbor_groups(radj, groups[r])
    ns, ms = compile_count_dicts(counts, groups[r], pairs)
    goal_ll = calculate_ll(ns, ms, method)
    results = []
    while groups_changed & (iter <= n_iter)
        groups_changed = false
        counts = count_neighbor_groups(radj, groups_guess)
        ns, ms = compile_count_dicts(counts, groups_guess, pairs)
        ll = calculate_ll(ns, ms, method)
        for i in randperm(k)
            count_deltas = counts[:,i][:]
            old_group = groups_guess[i]
            new_group = swap_group(old_group)
            new_ns, new_ms = adjust_counts(ns, ms, count_deltas, old_group, new_group)
            new_ll = calculate_ll(new_ns, new_ms, method)
            if new_ll - ll > 1e-4
                groups_guess[i] = new_group
                groups_changed = true
                break
            end
        end
        push!(results, [iter, ll, goal_ll, sum(abs(groups[r]-groups_guess))])
        iter += 1
        print("\r $iter")
    end
end


text:
using Test
using AdvancedHMC
using CUDA

@testset "AdvancedHMC GPU" begin
    include("common.jl")

    n_chains = 1000
    n_samples = 1000
    dim = 5

    T = Float32
    m, s, θ₀ = zeros(T, dim), ones(T, dim), rand(T, dim, n_chains)
    m, s, θ₀ = CuArray(m), CuArray(s), CuArray(θ₀)

    target = Gaussian(m, s)
    metric = UnitEuclideanMetric(T, size(θ₀))
    ℓπ, ∇ℓπ = get_ℓπ(target), get_∇ℓπ(target)
    hamiltonian = Hamiltonian(metric, ℓπ, ∇ℓπ)
    integrator = Leapfrog(one(T) / 5)
    proposal = HMCKernel(Trajectory{EndPointTS}(integrator, FixedNSteps(5)))

    samples, stats = sample(hamiltonian, proposal, θ₀, n_samples)
end


text:
function c2ir(c::AbstractVector, len=256)
    h = Array{eltype(c),1}(len)
    m = length(c) - 1

    h[1] = exp(c[1])
    for n=2:len
        d = 0.0
        upperlimit = (n >= m+1) ? m+1 : n
        for k=2:upperlimit
            @inbounds d += (k-1) * c[k] * h[n-k+1]
        end
        @inbounds h[n] = d / (n-1)
    end

    h
end


text:
using Test
using FlightMechanics
using FlightMechanics.Models
using FlightMechanics.Aircrafts

# Integration test for aircraft just making sure everything runs
ac = C310()

@test typeof(ac)<:Aircraft
@test get_name(ac) == "Cessna 310"

@test typeof(ac.propulsion)<:Propulsion

@test isapprox(get_wing_area(ac), 16.258032)

att = Attitude(1/180*pi, 0, 0)
pos = EarthPosition(0, 0, -1000)
state = State(pos, att, [65., 0., 3.], [0., 0., 0.], [0., 0., 0.], [0., 0., 0.])

env = Environment(pos, atmos="ISA1978", wind="NoWind", grav="const")
aerostate = AeroState(state, env)

controls = StickPedalsLeverControls(0.43, 0.562, 0.5, 0.68)

grav = get_gravity(env)
ac = calculate_aircraft(ac, controls, aerostate, state, grav; consume_fuel=false)

tas = 50  # m/s
h = 300  # m
psi = pi/3  # rad
gamma = 5 * DEG2RAD
turn_rate = 0.0

ac, aerostate, state, fcs = steady_state_trim(
    ac, controls, env, tas, pos, psi, gamma, turn_rate, 5.0*DEG2RAD, 0.0*DEG2RAD
    )

@test isapprox(ac.pfm.forces, zeros(3), atol=1e-5)
@test isapprox(ac.pfm.moments, zeros(3), atol=1e-5)
@test isapprox(get_tas(aerostate), tas, atol=1e-5)
@test isapprox(get_beta(aerostate), 0.0, atol=1e-5)


text:
# interface functions

"""
  Numbers all points on the specified surfaces.  All points on the
  surface of the domain not specified are numbered with the number of
  points on the specified surface + 1 (the numbering is 1-based)

  Currently this supports linear coordinate fields only

  **Inputs**

   * mesh
   * bc_nums: array of boundary condition numbers that define the surface(s)
              to number
   * isglobal: if true, do a globally consistent numbering of surface points,
               otherwise, do part-local numbering, default true
   * numbering_name: the name of the Pumi numbering to be created, any existing
                     numbering with this name will be deleted.
                     Default "warp_surf_nums"

  **Outputs**

   * numFaceNodes: number of nodes on the specified surface(s)
   * n_face: a Ptr{Void} (really an apf::Numbering*).  Note that this is
             not an apf::GlobalNumbering*, so be careful not to overflow
             a Cint.
   * face_verts: array of apf::MeshEntity* in the order they appear in the
                 surface numbering, length numFaceNodes
"""
function numberSurfacePoints(mesh::PumiMeshDG, bc_nums::AbstractVector{I}, isglobal::Bool=false, numbering_name::AbstractString="warp_surf_nums") where I<:Integer

  #TODO: should face_verts be only the MeshEntities owned by this part, or
  #      all the ones present on this part?

  @assert mesh.coord_order <= 2
  for i in bc_nums  # check that the BC numbers are valid
    @assert i <= length(mesh.bndry_offsets - 1)
    @assert i > 0
  end
  # number the (unique) coordinate nodes on the faces

  # we can't guarantee an existing numbering with the same name was
  # created using the same bc_nums, so delete any existing numbering
  n_old = apf.findNumbering(mesh.m_ptr, numbering_name)
  if n_old != C_NULL
    apf.destroyNumbering(n_old)
  end
  n_face = apf.createNumberingJ(mesh, numbering_name,
                                1, mesh.coordshape_ptr)
  topo = mesh.topo
  num_i = 1
  verts = Array{Ptr{Void}}(4)
  edges = Array{Ptr{Void}}(12)  # all edges of the element
  face_verts = Array{Ptr{Void}}(0)  # TODO: add sizehint
#  coords = zeros(Float64, 3)
  for i in bc_nums
    start_idx = mesh.bndry_offsets[i]
    end_idx = mesh.bndry_offsets[i+1] - 1
    rng = start_idx:end_idx

    for j in rng
      bndry_j = mesh.bndryfaces[j]
      el_j = bndry_j.element
      el_ptr = mesh.elements[el_j]

      # get the vertices
      apf.getDownward(mesh.m_ptr, el_ptr, 0, verts)

      for k=1:size(mesh.topo.face_verts, 1)
        v_k = verts[mesh.topo.face_verts[k, bndry_j.face]]

        if !apf.isNumbered(n_face, v_k, 0, 0)
          apf.numberJ(n_face, v_k, 0, 0, num_i)
          push!(face_verts, v_k)
          num_i += 1
        end
      end  # end loop k

      # get the edges nodes too
      if apf.hasNodesIn(mesh.coordshape_ptr, 1)
        apf.getDownward(mesh.m_ptr, el_ptr, 1, edges)
        for k=1:size(topo.face_edges, 1)
          edge_k = edges[topo.face_edges[k, bndry_j.face]]

          # up to 2nd order fields, we don't need to worry about edge orientation
          if !apf.isNumbered(n_face, edge_k, 0, 0)

            apf.numberJ(n_face, edge_k, 0, 0, num_i)
            push!(face_verts, edge_k)
            num_i += 1
          end  # end if
        end  # end loop k
      end  # end if

    end  # end loop j
  end  # end loop i

  # number all remaining points with the number of face verts + 1
  #TODO: this should be a loop 1:dim, iterative of meshentities of this dim
  for i=1:mesh.numVert
    v_i = mesh.verts[i]

    if !apf.isNumbered(n_face, v_i, 0, 0)
      apf.numberJ(n_face, v_i, 0, 0, num_i)
    end
  end

  if apf.hasNodesIn(mesh.coordshape_ptr, 1)
    for i=1:mesh.numEdge
      edge_i = mesh.edges[i]
      if !apf.isNumbered(n_face, edge_i, 0, 0)
        apf.numberJ(n_face, edge_i, 0, 0, num_i)
      end
    end
  end

  numFacePts = num_i - 1

  return numFacePts, n_face, face_verts
end

#------------------------------------------------------------------------------
# Reduction operations

"""
  Assignment reduction (always returns the right hand value)
"""
struct AssignReduction{T} <: Reduction{T}
  neutral_element::T

  function AssignReduction{T}() where {T}
    return new(0)
  end
end

function (obj::AssignReduction)(valL::Number, valR::Number)
  return valR
end

"""
  Sum reduction
"""
struct SumReduction{T} <: Reduction{T}
  neutral_element::T

  function SumReduction{T}() where {T}
    return new(0)
  end
end

function (obj::SumReduction)(valL::Number, valR::Number)
  return valL + valR
end

"""
  Minimum reduction
"""
struct MinReduction{T} <: Reduction{T}
  neutral_element::T

  function MinReduction{T}() where {T}
    return new(typemax(T))
  end
end

function (obj::MinReduction)(valL::Number, valR::Number)
  return min(valL, valR)
end

"""
  Max reduction
"""
struct MaxReduction{T}  <: Reduction{T}
  neutral_element::T

  function MaxReduction{T}() where {T}
    return new(typemin(T))
  end
end

function (obj::MaxReduction)(valL::Number, valR::Number)
  return max(valL, valR)
end

"""
  This function applies a user-supplied operation to take an array the same
  shape as `mesh.vert_coords` and transform it into a vector of length
  `mesh.dim*mesh.coord_numNodes` (a flattened array of all the coordinates).
  The user supplied operation is applied pairwise to all the values that are
  mapped to the same location in the vector.  `mesh.coord_nodenums_Nptr` is
  used to define the mapping from the array to the vector.

  The user supplied operation must have the signature

  ```
    val = reduce_op(valL, valR)
  ```

  where `valL` is the value that is already in the vector, and `valR` is
  the new value from `coords_arr`.  The vector is initialized with
  `reduce_op.neutral_element`

  **Inputs**

   * mesh: DG mesh object
   * coords_arr: array, same shape as `mesh.vert_coords`
   * reduce_op: a [`Reduction`](@ref) object, defaults to [`SumReduction`](@ref)

  **Inputs/Outputs**

   * coords_vec: vector to be populated, length `mesh.dim*mesh.coord_numNodes`

  **Keyword Arguments**

   * parallel: if true (default), do the reduction in parallel, sending data
               to the owning process.  The entry in the local vector for
               non-owned entities will be zero (in general,
               `reduce_op.neutral_element`)
"""
function coords3DTo1D(mesh::PumiMeshDG, coords_arr::AbstractArray{T, 3},
                     coords_vec::AbstractVector, reduce_op::Reduction{T}=SumReduction{T}(); parallel=true)  where {T}

  @assert mesh.coord_order <= 2
  @assert size(coords_arr, 3) == mesh.numEl
  @assert size(coords_arr, 2) == mesh.coord_numNodesPerElement
  @assert size(coords_arr, 1) == mesh.dim
  @assert length(coords_vec) == mesh.coord_numNodes*mesh.dim

  _parallel::Bool = parallel

  if _parallel
    sendParallelData(mesh.coordscatter, coords_arr, reduce_op)
  end

  fill!(coords_vec, reduce_op.neutral_element)
  node_entities = apf.ElementNodeEntities(mesh.m_ptr, mesh.coordshape_ptr, mesh.dim)

  #TODO: not sure if this gives enough time for data to arrive, maybe combine
  #      with calcCoordinatesAndMetrics_rev?
  shr = mesh.normalshr_ptr
  for i=1:mesh.numEl
    apf.getNodeEntities(node_entities, mesh.elements[i])
    for j=1:mesh.coord_numNodesPerElement
      entity = node_entities.entities[j]
      if !_parallel || (_parallel && apf.getOwner(shr, entity) == mesh.myrank)
        for k=1:mesh.dim
          idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, 0, k-1)
          orig_val = coords_vec[idx]
          coords_vec[idx] = reduce_op(coords_vec[idx], coords_arr[k, j, i])

        end
      end  # end if
    end  # end j
  end  # end i

  if _parallel
    # lambda function for receiving
    calc_func = (data::PeerData) -> receiveVecFunction(data, mesh, coords_vec,
                                                       reduce_op)

    receiveParallelData(mesh.coordscatter, calc_func)
  end

  return nothing
end

"""
  Like `coords1DTo3D`, but goes from the vector to the 3D array.  Unlike
  `coords1DTo3D`, this is an entirely local operation (does not do parallel
  communications)

  TODO: make this do parallel communication if reduce_op is not assignment

  **Inputs**

   * mesh: DG mesh
   * coords_vec: vector of coordinate-like values, length
                 `mesh.dim*mesh.coord_numNodes`
   * reduce_op: [`Reduction`](@ref) object, default [`AssignmnetReduction`](@ref)

  **Inputs/Outputs**

   * coords_arr: 3D array to be populated (overwritten)

  **Keyword Arguments**

   * parallel: if true, performs a global scatter such that the elements
               in coords_vec are overwritten by the owner's value before
               unpacking into `coords_arr`.
               default false.

  TODO: parallel efficiency of `parallel=true` could be improved
"""
function coords1DTo3D(mesh::PumiMeshDG, coords_vec::AbstractVector,
                      coords_arr::AbstractArray{T, 3},
                      reduce_op::Reduction=AssignReduction{T}();
                      parallel::Bool=false) where {T}

  @assert mesh.coord_order <= 2
  @assert size(coords_arr, 3) == mesh.numEl
  @assert size(coords_arr, 2) == mesh.coord_numNodesPerElement
  @assert size(coords_arr, 1) == mesh.dim
  @assert length(coords_vec) == mesh.coord_numNodes*mesh.dim

  if parallel
    sendParallelData_rev(mesh, mesh.coordscatter, coords_vec)
  end

  fill!(coords_arr, reduce_op.neutral_element)
  node_entities = apf.ElementNodeEntities(mesh.m_ptr, mesh.coordshape_ptr, mesh.dim)

  if parallel
    calc_func = (data::PeerData) -> receiveFromOwner(data, mesh, coords_vec)
    receiveParallelData_rev(mesh.coordscatter, calc_func)
  end

  for i=1:mesh.numEl
    apf.getNodeEntities(node_entities, mesh.elements[i])
    for j=1:mesh.coord_numNodesPerElement
      entity = node_entities.entities[j]
      for k=1:mesh.dim
        idx = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, 0, k-1)

        coords_arr[k, j, i] = reduce_op(coords_arr[k, j, i], coords_vec[idx])
      end
    end
  end

  return nothing
end


text:
module inmet

using HTTP
using CSV
using JSON3
import Dates
using DataFrames

"""
Executa uma chamada GET simples a uma API e retorna uma resposta JSON.

# Examples
```jldoctest
julia> url = "https://apitempo.inmet.gov.br/estacao/2021-01-01/2021-01-02/A756"
julia> resposta = chamarAPI(url)
julia> println(resposta[1])
```
"""
function chamarAPI(url::String)::JSON3.Array
    resposta = HTTP.get(url)
    respostaString = String(resposta.body)
    return JSON3.read(respostaString, allow_inf=true)
end;

"""
Lista as estações meteorológicas no INMET.

# Examples
```jldoctest
julia> estacoes = listarEstacoes(true)
julia> println(estacoes[:, 1])
```
"""
function listarEstacoes(operantes::Bool)::DataFrame
    url = "https://apitempo.inmet.gov.br/estacoes/T"
    listaEstacoesAutomaticasJSON = chamarAPI(url)
    listaEstacoesAutomaticas = DataFrame(listaEstacoesAutomaticasJSON);
    if operantes == true
        listaEstacoesAutomaticasOperantes = listaEstacoesAutomaticas[
            listaEstacoesAutomaticas[:, "CD_SITUACAO"] .== "Operante", :];
    end

    return listaEstacoesAutomaticas
end

"""
Coleta os dados de uma estação INMET.

# Examples
```jldoctest
julia> formato = Dates.DateFormat("y-m-d")
julia> dataInicial = Dates.Date("2021-01-01", formato)
julia> dataFinal = Dates.Date("2021-01-02", formato)
julia> dados = coletarDadosEstacao("A756", dataInicial, dataFinal)
julia> println(resposta[1])
```
"""
function coletarDadosEstacao(
        codigoEstacao::String,
        dataInicial::Dates.Date,
        dataFinal::Dates.Date
)::DataFrame
    urlBaseDadosEstacoes = "https://apitempo.inmet.gov.br/estacao"
    urlQueryDadosEstacoes = "$(dataInicial)/$(dataFinal)/$(codigoEstacao)"
    urlCompleta = "$(urlBaseDadosEstacoes)/$(urlQueryDadosEstacoes)"

    serieEstacaoJSON = chamarAPI(urlCompleta)
    serieEstacao = DataFrame(serieEstacaoJSON)

    # substitui valores 'nothing' por 'missing'
    for coluna in names(serieEstacao)
        serieEstacao[!, coluna] = replace(
            serieEstacao[!, coluna], nothing => missing
        )
    end

    # se o dataframe gerador for vazio, a função irá parar
    if size(dropmissing(serieEstacao))[1] == 0
        error("não há dados disponíveis para a estação '$(codigoEstacao)'.")
    end

    return serieEstacao
end;

"""
Coleta os dados de todas as estações (operantes) e persiste as saídas em formato CSV.

# Examples
```jldoctest
julia> formato = Dates.DateFormat("y-m-d")
julia> dataInicial = Dates.Date("2021-01-01", formato)
julia> dataFinal = Dates.Date("2021-01-02", formato)
julia> coletarDadosEstacoesCSV(
    ["A756", "A223"], dados_estacoes", dataInicial, dataFinal
)
```
"""
function coletarDadosEstacoesCSV(
    dataInicial::Dates.Date,
    dataFinal::Dates.Date,
    diretorioPersistencia::String,
    codigosEstacoes::Union{Vector, Nothing} = nothing,
)

    if !Base.isdir(diretorioPersistencia)
        Base.mkpath(diretorioPersistencia)
    end

    if codigosEstacoes === nothing
        codigosEstacoes = listarEstacoes(true)
    end

    for codigoEstacao in codigosEstacoes
        try
            serieEstacao = coletarDadosEstacao(codigoEstacao, dataInicial, dataFinal)
            CSV.write("$(diretorioPersistencia)/$(codigoEstacao).csv", serieEstacao)
            println("Dados da estação '$(codigoEstacao)' obtidos.")
        catch e
            println("Não foi possível processar os dados de '$(codigoEstacao)': $(e)")
        end
        #=
        Intervalos de tempo são necessários para que a coleta seja feita com sucesso.
        Caso contrário, a API INMET passa a retornar código 200 com mensagens de erro,
        vetando todas as requisições por um período de, normalmente, 1 minuto.
        =#
        Base.sleep(1)
    end
end;

end # module

max_stars_repo_name: ericmiguel/inmet-julia

max_issues_repo_name: ericmiguel/inmet-julia

max_forks_repo_name: ericmiguel/inmet-julia


text:
# Julia wrapper for header: /usr/include/libavutil/pixdesc.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_bits_per_pixel,
    av_get_padded_bits_per_pixel,
    av_pix_fmt_desc_get,
    av_pix_fmt_desc_next,
    av_pix_fmt_desc_get_id,
    av_pix_fmt_get_chroma_sub_sample,
    av_pix_fmt_count_planes,
    av_color_range_name,
    av_color_range_from_name,
    av_color_primaries_name,
    av_color_primaries_from_name,
    av_color_transfer_name,
    av_color_transfer_from_name,
    av_color_space_name,
    av_color_space_from_name,
    av_chroma_location_name,
    av_chroma_location_from_name,
    av_get_pix_fmt,
    av_get_pix_fmt_name,
    av_get_pix_fmt_string,
    av_read_image_line2,
    av_read_image_line,
    av_write_image_line2,
    av_write_image_line,
    av_pix_fmt_swap_endianness,
    av_get_pix_fmt_loss,
    av_find_best_pix_fmt_of_2

function av_get_bits_per_pixel(pixdesc)
    ccall((:av_get_bits_per_pixel, libavutil), Cint, (Ptr{AVPixFmtDescriptor},), pixdesc)
end

function av_get_padded_bits_per_pixel(pixdesc)
    ccall((:av_get_padded_bits_per_pixel, libavutil), Cint, (Ptr{AVPixFmtDescriptor},), pixdesc)
end

function av_pix_fmt_desc_get(pix_fmt::AVPixelFormat)
    ccall((:av_pix_fmt_desc_get, libavutil), Ptr{AVPixFmtDescriptor}, (AVPixelFormat,), pix_fmt)
end

function av_pix_fmt_desc_next(prev)
    ccall((:av_pix_fmt_desc_next, libavutil), Ptr{AVPixFmtDescriptor}, (Ptr{AVPixFmtDescriptor},), prev)
end

function av_pix_fmt_desc_get_id(desc)
    ccall((:av_pix_fmt_desc_get_id, libavutil), AVPixelFormat, (Ptr{AVPixFmtDescriptor},), desc)
end

function av_pix_fmt_get_chroma_sub_sample(pix_fmt::AVPixelFormat, h_shift, v_shift)
    ccall((:av_pix_fmt_get_chroma_sub_sample, libavutil), Cint, (AVPixelFormat, Ptr{Cint}, Ptr{Cint}), pix_fmt, h_shift, v_shift)
end

function av_pix_fmt_count_planes(pix_fmt::AVPixelFormat)
    ccall((:av_pix_fmt_count_planes, libavutil), Cint, (AVPixelFormat,), pix_fmt)
end

function av_color_range_name(range::AVColorRange)
    ccall((:av_color_range_name, libavutil), Cstring, (AVColorRange,), range)
end

function av_color_range_from_name(name)
    ccall((:av_color_range_from_name, libavutil), Cint, (Cstring,), name)
end

function av_color_primaries_name(primaries::AVColorPrimaries)
    ccall((:av_color_primaries_name, libavutil), Cstring, (AVColorPrimaries,), primaries)
end

function av_color_primaries_from_name(name)
    ccall((:av_color_primaries_from_name, libavutil), Cint, (Cstring,), name)
end

function av_color_transfer_name(transfer::AVColorTransferCharacteristic)
    ccall((:av_color_transfer_name, libavutil), Cstring, (AVColorTransferCharacteristic,), transfer)
end

function av_color_transfer_from_name(name)
    ccall((:av_color_transfer_from_name, libavutil), Cint, (Cstring,), name)
end

function av_color_space_name(space::AVColorSpace)
    ccall((:av_color_space_name, libavutil), Cstring, (AVColorSpace,), space)
end

function av_color_space_from_name(name)
    ccall((:av_color_space_from_name, libavutil), Cint, (Cstring,), name)
end

function av_chroma_location_name(location::AVChromaLocation)
    ccall((:av_chroma_location_name, libavutil), Cstring, (AVChromaLocation,), location)
end

function av_chroma_location_from_name(name)
    ccall((:av_chroma_location_from_name, libavutil), Cint, (Cstring,), name)
end

function av_get_pix_fmt(name)
    ccall((:av_get_pix_fmt, libavutil), AVPixelFormat, (Cstring,), name)
end

function av_get_pix_fmt_name(pix_fmt::AVPixelFormat)
    ccall((:av_get_pix_fmt_name, libavutil), Cstring, (AVPixelFormat,), pix_fmt)
end

function av_get_pix_fmt_string(buf, buf_size::Integer, pix_fmt::AVPixelFormat)
    ccall((:av_get_pix_fmt_string, libavutil), Cstring, (Cstring, Cint, AVPixelFormat), buf, buf_size, pix_fmt)
end

function av_read_image_line2(dst, data::NTuple{4, Ptr{UInt8}}, linesize::NTuple{4, Cint}, desc, x::Integer, y::Integer, c::Integer, w::Integer, read_pal_component::Integer, dst_element_size::Integer)
    ccall((:av_read_image_line2, libavutil), Cvoid, (Ptr{Cvoid}, NTuple{4, Ptr{UInt8}}, NTuple{4, Cint}, Ptr{AVPixFmtDescriptor}, Cint, Cint, Cint, Cint, Cint, Cint), dst, data, linesize, desc, x, y, c, w, read_pal_component, dst_element_size)
end

function av_read_image_line(dst, data::NTuple{4, Ptr{UInt8}}, linesize::NTuple{4, Cint}, desc, x::Integer, y::Integer, c::Integer, w::Integer, read_pal_component::Integer)
    ccall((:av_read_image_line, libavutil), Cvoid, (Ptr{UInt16}, NTuple{4, Ptr{UInt8}}, NTuple{4, Cint}, Ptr{AVPixFmtDescriptor}, Cint, Cint, Cint, Cint, Cint), dst, data, linesize, desc, x, y, c, w, read_pal_component)
end

function av_write_image_line2(src, data::NTuple{4, Ptr{UInt8}}, linesize::NTuple{4, Cint}, desc, x::Integer, y::Integer, c::Integer, w::Integer, src_element_size::Integer)
    ccall((:av_write_image_line2, libavutil), Cvoid, (Ptr{Cvoid}, NTuple{4, Ptr{UInt8}}, NTuple{4, Cint}, Ptr{AVPixFmtDescriptor}, Cint, Cint, Cint, Cint, Cint), src, data, linesize, desc, x, y, c, w, src_element_size)
end

function av_write_image_line(src, data::NTuple{4, Ptr{UInt8}}, linesize::NTuple{4, Cint}, desc, x::Integer, y::Integer, c::Integer, w::Integer)
    ccall((:av_write_image_line, libavutil), Cvoid, (Ptr{UInt16}, NTuple{4, Ptr{UInt8}}, NTuple{4, Cint}, Ptr{AVPixFmtDescriptor}, Cint, Cint, Cint, Cint), src, data, linesize, desc, x, y, c, w)
end

function av_pix_fmt_swap_endianness(pix_fmt::AVPixelFormat)
    ccall((:av_pix_fmt_swap_endianness, libavutil), AVPixelFormat, (AVPixelFormat,), pix_fmt)
end

function av_get_pix_fmt_loss(dst_pix_fmt::AVPixelFormat, src_pix_fmt::AVPixelFormat, has_alpha::Integer)
    ccall((:av_get_pix_fmt_loss, libavutil), Cint, (AVPixelFormat, AVPixelFormat, Cint), dst_pix_fmt, src_pix_fmt, has_alpha)
end

function av_find_best_pix_fmt_of_2(dst_pix_fmt1::AVPixelFormat, dst_pix_fmt2::AVPixelFormat, src_pix_fmt::AVPixelFormat, has_alpha::Integer, loss_ptr)
    ccall((:av_find_best_pix_fmt_of_2, libavutil), AVPixelFormat, (AVPixelFormat, AVPixelFormat, AVPixelFormat, Cint, Ptr{Cint}), dst_pix_fmt1, dst_pix_fmt2, src_pix_fmt, has_alpha, loss_ptr)
end


text:
using HPAT

@acc hpat function unique_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Int64}, HDF5, file_name)
    df1 = aggregate(df, :id, :sx=length(unique(:x)))
    return sum(df1[:sx])
end

unique_test(ENV["HOME"]*"/tmp/unique.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
module PlayingCards

using Random: randperm
import Random: shuffle!

import Base

# Suits
export ♣, ♠, ♡, ♢ # aliases

# Card, and Suit
export Card, Suit

# Card properties
export suit, rank, high_value, low_value, color

# Lists of all ranks / suits
export ranks, suits

# Deck & deck-related methods
export Deck, shuffle!, full_deck, ordered_deck

#####
##### Types
#####

"""
    Suit

Encode a suit as a 2-bit value (low bits of a `UInt8`):
- 0 = ♣ (clubs)
- 1 = ♢ (diamonds)
- 2 = ♡ (hearts)
- 3 = ♠ (spades)

Suits have global constant bindings: `♣`, `♢`, `♡`, `♠`.
"""
struct Suit
    i::UInt8
    Suit(s::Integer) = 0 ≤ s ≤ 3 ? new(s) :
        throw(ArgumentError("invalid suit number: $s"))
end

#=
```
    char(::Suit)

julia> Char.(0x2663 .- UInt8.(0:3))
4-element Array{Char,1}:
 '♣': Unicode U+2663 (category So: Symbol, other)
 '♢': Unicode U+2662 (category So: Symbol, other)
 '♡': Unicode U+2661 (category So: Symbol, other)
 '♠': Unicode U+2660 (category So: Symbol, other)
```
=#

"""
    char

Return the unicode characters:
"""
char(s::Suit) = Char(0x2663-s.i)
Base.string(s::Suit) = string(char(s))
Base.show(io::IO, s::Suit) = print(io, char(s))

"""
    Card

Encode a playing card as a 6-bit integer (low bits of a `UInt8`):
- low bits represent rank from 0 to 15
- high bits represent suit (♣, ♢, ♡ or ♠)

Ranks are assigned as follows:
- numbered cards (2 to 10) have rank equal to their number
- jacks, queens and kings have ranks 11, 12 and 13
- there are low and high aces with ranks 1 and 14
- there are low and high jokers with ranks 0 and 15

This allows any of the standard orderings of cards ranks to be
achieved simply by choosing which aces to use.

There are a total of 64 possible card values with this scheme,
represented by `UInt8` values `0x00` through `0x3f`.
"""
struct Card
    value::UInt8
    function Card(r::Integer, s::Integer)
        1 ≤ r ≤ 13 || throw(ArgumentError("invalid card rank: $r"))
        left_bits = UInt8(s << 4)
        right_bits = UInt8(r)
        or_bits = left_bits | right_bits
        return new(or_bits)
    end
end

Card(r::Integer, s::Suit) = Card(r, s.i)

#=
```julia
value(r, s) = UInt8(s << 4) | UInt8(r)
suit(r, s) = (0x30 & value(r, s)) >>> 4
julia> suit(1, 0)+0 # 0
julia> suit(1, 1)+0 # 1
julia> suit(1, 2)+0 # 2
julia> suit(1, 3)+0 # 3
```
=#

"""
    suit(::Card)

The suit of a card
"""
suit(c::Card) = Suit((0x30 & c.value) >>> 4)

#=
```julia
value(r, s) = UInt8(s << 4) | UInt8(r)
rank(r, s) = Int8((0x0f & value(r, s)))
julia> rank(1, 0)+0 # 1
julia> rank(1, 3)+0 # 1
julia> rank(9, 3)+0 # 9
```
=#

"""
    rank(::Card)

The rank of a card
"""
rank(c::Card) = Int8((c.value & 0x0f))

const ♣ = Suit(0)
const ♢ = Suit(1)
const ♡ = Suit(2)
const ♠ = Suit(3)

# Allow constructing cards with, e.g., `3♡`
Base.:*(r::Integer, s::Suit) = Card(r, s)

function Base.show(io::IO, c::Card)
    r = rank(c)
    @assert 1 ≤ r ≤ 14
    if r == 1
        print(io, 'A')
    else
        print(io, "123456789TJQK"[r])
    end
    print(io, suit(c))
end

# And for face cards:
# Not to be confused with other unicode characters:
# ♡, ♡
# ♢, ♢
for s in "♣♢♡♠", (f,typ) in zip((:T,:J,:Q,:K,:A),(10,11,12,13,1))
    ss, sc = Symbol(s), Symbol("$f$s")
    @eval (export $sc; const $sc = Card($typ,$ss))
end

#####
##### Methods
#####

function rank_string(r::Int8)
    @assert 1 ≤ r ≤ 13
    if r ≤ 9
        2 ≤ r && return "$(r)"
        return "A"
    else
        if r < 12
            r == 11 && return "J"
            return "T"
        else
            r == 12 && return "Q"
            return "K"
        end
    end
end

Base.string(card::Card) = rank_string(rank(card))*string(suit(card))

"""
    high_value(::Card)
    high_value(::Rank)

The high rank value. For example:
 - `Rank(1)` -> 14 (use [`low_value`](@ref) for the low Ace value.)
 - `Rank(5)` -> 5
"""
high_value(c::Card) = rank(c) == 1 ? 14 : rank(c)

"""
    low_value(::Card)
    low_value(::Rank)

The low rank value. For example:
 - `Rank(1)` -> 1 (use [`high_value`](@ref) for the high Ace value.)
 - `Rank(5)` -> 5
"""
low_value(c::Card) = rank(c)

"""
    color(::Card)

A `Symbol` (`:red`, or `:black`) indicating
the color of the suit or card.
"""
function color(s::Suit)
    if s == ♣ || s == ♠
        return :black
    else
        return :red
    end
end
color(card::Card) = color(suit(card))

#####
##### Full deck/suit/rank methods
#####

"""
    ranks

A Tuple of ranks `1:13`.
"""
ranks() = 1:13

"""
    suits

A Tuple of all suits
"""
suits() = (♣, ♠, ♡, ♢)

"""
    full_deck

A vector of a cards
containing a full deck
"""
full_deck() = Card[Card(r,s) for s in suits() for r in ranks()]

#### Deck

"""
    Deck

Deck of cards (backed by a `Vector{Card}`)
"""
struct Deck{C <: Vector}
    cards::C
end

Base.length(deck::Deck) = length(deck.cards)

Base.iterate(deck::Deck, state=1) = Base.iterate(deck.cards, state)

function Base.show(io::IO, deck::Deck)
    for (i, card) in enumerate(deck)
        Base.show(io, card)
        if mod(i, 13) == 0
            println(io)
        else
            print(io, " ")
        end
    end
end

"""
    pop!(deck::Deck, n::Int = 1)
    pop!(deck::Deck, card::Card)

Remove `n` cards from the `deck`.
or
Remove `card` from the `deck`.
"""
Base.pop!(deck::Deck, n::Integer = 1) = ntuple(i->pop!(deck.cards), n)
function Base.pop!(deck::Deck, card::Card)
    L0 = length(deck)
    filter!(x -> x ≠ card, deck.cards)
    L0 == length(deck)+1 || error("Could not pop $(card) from deck.")
    return card
end

"""
    ordered_deck

An ordered `Deck` of cards.
"""
ordered_deck() = Deck(full_deck())

"""
    shuffle!

Shuffle the deck! `shuffle!` uses
`Random.randperm` to shuffle the deck.
"""
function shuffle!(deck::Deck)
    deck.cards .= deck.cards[randperm(length(deck.cards))]
    deck
end

end # module


text:
# indicator of the Stiefel manifold

export IndStiefel

@doc raw"""
**Indicator of the Stiefel manifold**

    IndStiefel()

Returns the indicator of the Stiefel manifold
```math
S_{n,p} = \left\{ X \in \mathbb{F}^{n \times p} : X^*X = I \right\}.
```
where ``\mathbb{F}`` is the real or complex field, and parameters ``n`` and ``p``
are inferred from the matrix provided as input.
"""
struct IndStiefel <: ProximableFunction end

is_set(f::IndStiefel) = true

function (f::IndStiefel)(X::AbstractMatrix{T}) where {R <: Real, T <: Union{R, Complex{R}}}
    F = svd(X)
    if all(F.S .≈ R(1))
        return R(0)
    end
    return R(Inf)
end

function prox!(Y::AbstractMatrix{T}, f::IndStiefel, X::AbstractMatrix{T}, gamma::R=R(1)) where {R <: Real, T <: Union{R, Complex{R}}}
    n, p = size(X)
    F = svd(X)
    U_sliced = view(F.U, :, 1:p)
    mul!(Y, U_sliced, F.Vt)
    return R(0)
end

function prox_naive(f::IndStiefel, X::AbstractMatrix{T}, gamma::R=R(1)) where {R, T <: Union{R, Complex{R}}}
    n, p = size(X)
    F = svd(X)
    Y = F.U[:, 1:p] * F.Vt
    return Y, R(0)
end


text:
using FourierTools, FFTW, PGFPlotsX

N_low = 64
x_min = 0.0
x_max = 8*2π

xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
xs_high = range(x_min, x_max, length=5000)[1:end-1]
f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
arr_low = f.(xs_low)
arr_high = f.(xs_high)

N = 1000
xs_interp = range(x_min, x_max, length=N+1)[1:N]
arr_interp = resample(arr_low, N)

plt1 = @pgf PGFPlotsX.Axis(
    {
        xmin=0, xmax=50, ymin=-3, ymax=4.5,
        legend_pos="north east", legend_entries={"low sampling", "sinc interpolated", "high sampling"},
        width = "10cm",
        height= "6cm",
    },

    Plot({mark="x", style="only marks", color="green",mark_options = {scale=1.3},
        },
        Table([xs_low, arr_low])
    ),

    Plot({"sharp plot", color="red"
        },
        Table([xs_interp, arr_interp])
    ),

    Plot({"sharp plot", style="dashed", color="blue"
        },
        Table([xs_high, arr_high])
    ),
)
pgfsave("../../paper/figures/resampling.pdf", plt1)


text:
using Test
using NQCDynamics
using LinearAlgebra
using Random
using Distributions
using NQCDynamics: DynamicsMethods, DynamicsUtils, Calculators
using NQCDynamics.DynamicsMethods: SurfaceHoppingMethods
using ComponentArrays

n_states = 20
temperature = 9.5e-4
model = MiaoSubotnik(M=n_states-1)
atoms = Atoms(2000)
r = zeros(1, 1)
v = randn(1,1) * sqrt(5*temperature / atoms.masses[1])
n_electrons = nstates(model) ÷ 2

sim = Simulation{IESH}(atoms, model; n_electrons=n_electrons)
u = DynamicsVariables(sim, v, r)
sim.method.state .= u.state

@testset "compute_overlap!" begin
    S = zeros(n_electrons, n_electrons)
    state = collect(1:n_electrons)
    ψ = DynamicsUtils.get_quantum_subsystem(u)
    SurfaceHoppingMethods.compute_overlap!(S, ψ, state)
    @test S == I # Check overlap is identity

    # All electrons in state 1, so they all have unity overlap
    # with first electron wavefunction, which is in state 1.
    state = ones(Int, n_electrons)
    SurfaceHoppingMethods.compute_overlap!(S, ψ, state)
    @test all(S[:,1] .== 1) # Check only first column ones
    @test all(S[:,2:end] .== 0)
end

@testset "calculate_Akj" begin
    S = zeros(n_electrons, n_electrons)
    state = collect(1:n_electrons)
    ψ = DynamicsUtils.get_quantum_subsystem(u)
    Akj = SurfaceHoppingMethods.calculate_Akj(S, ψ, 1.0, state)
    @test Akj ≈ 1.0
end

@testset "evaluate_hopping_probability!" begin
    Calculators.update_electronics!(sim.calculator, hcat(20.0))
    z = deepcopy(u)
    rand!(DynamicsUtils.get_quantum_subsystem(z))
    SurfaceHoppingMethods.evaluate_hopping_probability!(sim, z, 1.0)
end

BLAS.set_num_threads(1)

temp = 5*temperature
harm = Harmonic(m=model.m, ω=model.ω)
z = ComponentVector(v=zeros(1,1), r=hcat(10.0))
thermal = Simulation{Langevin}(atoms, harm; temperature=temp, γ=0.1)
Δ = Dict(:X => 1.0)
sample = InitialConditions.MetropolisHastings.run_monte_carlo_sampling(thermal, hcat(0.0), Δ, 1e6)

vel = Normal(0, sqrt(temp / atoms.masses[1]))

dist = DynamicalDistribution(vel, sample.R, (1,1))

using Plots
plotly()

@time res = run_ensemble(sim, (0.0, 5e4), dist;
    output=(:position, :population, :adiabatic_population), trajectories=2, abstol=1e-7, reltol=1e-4, saveat=0.0:50:5e4)

# This should plot the impurity population for each of the trajectories.
# It should look like the inset of fig 5.c in Miao Subotnik.
# Unfortunately it does not :(
impurity_population = zeros(length(0:50:5e4), length(res))
for (i, traj) in enumerate(res)
    # plot!(traj.t.*model.ω, 1 .- [p[1] for p in traj.population])
    impurity_population[:,i] .= 1 .- [p[1] for p in traj.population]
end

avg = mean(impurity_population, dims=2)
plot(res[1].t.*model.ω, avg)


text:
using SafeTestsets

@time begin
@time @safetestset "Fast Broadcast" begin include("fastbc.jl") end
@time @safetestset "Number of Parameters Calculation" begin include("numargs_test.jl") end
@time @safetestset "Data Arrays" begin include("data_array_tests.jl") end
@time @safetestset "Extended Functions" begin include("extended_function_tests.jl") end
@time @safetestset "Callbacks" begin include("callbacks.jl") end
@time @safetestset "Plot Variables" begin include("plot_vars.jl") end
@time @safetestset "Problem Creation Tests" begin include("problem_creation_tests.jl") end
@time @safetestset "Remake tests" begin include("remake_tests.jl") end
@time @safetestset "Affine differential equation operators" begin include("affine_operators_tests.jl") end
@time @safetestset "TableTraits" begin include("tabletraits_tests.jl") end
@time @safetestset "Integrator interface" begin include("integrator_tests.jl") end
@time @safetestset "Export tests" begin include("export_tests.jl") end
@time @safetestset "High Level solve Interface" begin include("high_level_solve.jl") end
end