text:
"""
Sampling from DPP and k-DPP.
References:
-----------
    [1] Kulesza, A., and B. Taskar. Determinantal point processes for machine
        learning. arXiv preprint arXiv:1207.6083, 2012.
    [2] Kang, B. Fast determinantal point process sampling with application to
        clustering. NIPS, 2013.
"""

"""
    _sample_mask!(rng, M, Λ, i)

Sample a mask for an elementary DPP.
"""
function _sample_mask!(
    rng::AbstractRNG, M::AbstractMatrix{Bool}, Λ::AbstractArray{<:Real}, i::Int
)
    for j in 1:length(Λ)
        M[j, i] = rand(rng) < (Λ[j] / (Λ[j] + 1))
    end
    return M
end

"""
    _sample_k_mask!(rng, M, Λ, E, k, i)

Sample a mask for an elementary k-DPP.
"""
function _sample_k_mask!(
    rng::AbstractRNG,
    M::AbstractMatrix{Bool},
    Λ::AbstractArray{<:Real},
    E::AbstractMatrix{<:Real},
    k::Int,
    i::Int,
)
    j = length(Λ)
    remaining = k

    # iteratively sample a k-mask
    while remaining > 0
        # compute marginal of j given that we choose remaining values from 1:j
        if j == remaining
            marg = 1
        else
            marg = Λ[j] * E[remaining, j] / E[remaining + 1, j + 1]
        end

        # sample marginal
        if rand(rng) <= marg
            M[j, i] = true
            remaining -= 1
        end
        j -= 1
    end
    return M
end

"""
    _sample_from_elementary(rng, V, M, i)

Exact sampling from an elementary DPP. The algorithm based on [1].
"""
function _sample_from_elementary(
    rng::AbstractRNG, V::AbstractMatrix{T}, M::AbstractMatrix{Bool}, i::Int
) where {T<:Real}
    # select the elementary DPP
    V_mask = M[:, i]

    # edge case: empty sample
    if isempty(V_mask) # !any(V_mask)
        return Int[]
    end

    # select the kernel of the elementary DPP
    L = V[:, V_mask]

    Y = Int[]
    mask = ones(Bool, size(L, 2))
    prob = zeros(T, size(L, 1))

    for i in 1:size(L, 2)
        # compute probabilities
        fill!(prob, zero(T))
        for c in 1:size(L, 2)
            !mask[c] && continue
            for r in 1:size(L, 1)
                prob[r] += L[r, c] .^ 2
            end
        end
        prob ./= sum(prob)

        # sample a point in the original space
        h = findfirst(rand(rng) .<= cumsum(prob))
        push!(Y, h)

        # select and mask-out an element
        j = get_first_nz_idx(L[h, :], mask)
        mask[j] = false

        if any(mask)
            # Subtract scaled Lj from other columns so that their
            # projections on e_s[i] turns into 0. This operation
            # preserves the rank of L_{-j}.
            for c in 1:size(L, 2)
                !mask[c] && continue
                for r in 1:size(L, 1)
                    L[r, c] -= L[r, j] * L[h, c] / L[h, j]
                end
            end

            # Gram-Schmidt orthogonalization
            L[:, mask] = Matrix(qr(L[:, mask]).Q)
        end
    end

    return sort(Y)
end

Random.rand(pp::PointProcess, n::Int) = rand(GLOBAL_RNG, pp, n)
Random.rand(pp::PointProcess) = rand(GLOBAL_RNG, pp)
Random.rand(rng::AbstractRNG, pp::PointProcess) = first(rand(rng, pp, 1))

"""
    rand([rng::AbstractRNG], dpp::DeterminantalPointProcess)::Vector{Int}
    rand([rng::AbstractRNG], dpp::DeterminantalPointProcess, n::Int)::Vector{Vector{Int}}

Exact sampling from a DPP [1].
Returns a vector of indices with respect to the `L` matrix passed to the `DeterminantalPointProcess`.
The length of each vector can vary from 0 to the `size(L,1)`
"""
function Random.rand(
    rng::AbstractRNG, dpp::DeterminantalPointProcess{T,true}, N::Int) where {T<:Real}
    Λ = SharedArray{T}(dpp.Lfact.values)
    V = SharedMatrix{T}(dpp.Lfact.vectors)
    M = SharedMatrix{Bool}(zeros(Bool, dpp.size, N))

    # step I: sample masks for elementary DPPs
    pmap(
        (i, seed) -> _sample_mask!(MersenneTwister(seed), M, Λ, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )

    # step II: iteratively sample from a mixture of elementary DPPs
    return pmap(
        (i, seed) -> _sample_from_elementary(MersenneTwister(seed), V, M, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )
end

function Random.rand(
    rng::AbstractRNG, dpp::DeterminantalPointProcess{T,false}, N::Int) where {T<:Real}
    Λ = Array{T}(dpp.Lfact.values)
    V = Matrix{T}(dpp.Lfact.vectors)
    M = Matrix{Bool}(zeros(Bool, dpp.size, N))

    # step I: sample masks for elementary DPPs
    map(
        (i, seed) -> _sample_mask!(MersenneTwister(seed), M, Λ, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )

    # step II: iteratively sample from a mixture of elementary DPPs
    return map(
        (i, seed) -> _sample_from_elementary(MersenneTwister(seed), V, M, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )
end

"""
    rand([rng::AbstractRNG], kdpp::kDeterminantalPointProcess)::Vector{Int}
    rand([rng::AbstractRNG], kdpp::kDeterminantalPointProcess, n::Int)::Vector{Vector{Int}}

Exact sampling from a k-DPP [1].
Returns a vector of indices with respect to the `L` matrix passed to the `kDeterminantalPointProcess`
Each vector is of size `k`.
"""
function Random.rand(rng::AbstractRNG, kdpp::kDPP{T,true}, N::Int) where {T<:Real}
    dpp = kdpp.dpp
    Λ = SharedArray{T}(dpp.Lfact.values)
    V = SharedMatrix{T}(dpp.Lfact.vectors)
    M = SharedMatrix{Bool}(zeros(Bool, dpp.size, N))

    # compute elementary symmetric polynomials
    E = SharedMatrix{T}(elem_symm_poly(dpp.Lfact.values, kdpp.k))

    # step I: sample masks for elementary DPPs
    pmap(
        (i, seed) -> _sample_k_mask!(MersenneTwister(seed), M, Λ, E, kdpp.k, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )

    # step II: iteratively sample from a mixture of elementary DPPs
    return pmap(
        (i, seed) -> _sample_from_elementary(MersenneTwister(seed), V, M, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )
end

function Random.rand(rng::AbstractRNG, kdpp::kDPP{T,false}, N::Int) where {T<:Real}
    dpp = kdpp.dpp
    Λ = Array{T}(dpp.Lfact.values)
    V = Matrix{T}(dpp.Lfact.vectors)
    M = Matrix{Bool}(zeros(Bool, dpp.size, N))

    # compute elementary symmetric polynomials
    E = Matrix{T}(elem_symm_poly(dpp.Lfact.values, kdpp.k))

    # step I: sample masks for elementary DPPs
    map(
        (i, seed) -> _sample_k_mask!(MersenneTwister(seed), M, Λ, E, kdpp.k, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )

    # step II: iteratively sample from a mixture of elementary DPPs
    return map(
        (i, seed) -> _sample_from_elementary(MersenneTwister(seed), V, M, i),
        1:N,
        abs.(rand(rng, Int, N)),
    )
end

"""
    Perform one MCMC accept-reject transition for DPP.
"""
function _do_mcmc_step!(rng::AbstractRNG, dpp::DeterminantalPointProcess, state::MCMCState)
    # propose an element to swap
    u = rand(rng, 1:(dpp.size))
    insert = !state[1][u]

    # attempt to make a transition
    if insert
        p = _comp_accept_prob(dpp, state, u, insert)
        rand(rng) < p && _update_mcmc_state!(dpp, state, u, insert)
    else  # delete
        new_state = _update_mcmc_state(dpp, state, u, insert)
        p = _comp_accept_prob(dpp, new_state, u, insert)
        if rand(rng) < p
            state[1] .= new_state[1]
            state[2] .= new_state[2]
        end
    end
end

"""
    Perform one MCMC accept-reject transition for k-DPP.
"""
function _do_mcmc_k_step!(rng::AbstractRNG, kdpp::kDPP, state::MCMCState)
    z, L_z_inv = state

    # propose the elements to swap
    u, v = rand(rng, findall(z)), rand(rng, findall(z .== true))

    # copy the state and delete the u element
    new_state = _update_mcmc_state(kdpp.dpp, state, u, false)

    # attempt to make a transition
    p = _comp_accept_prob(kdpp.dpp, new_state, u, v)
    if rand(rng) < p
        # insert the v element into the new state
        _update_mcmc_state!(kdpp.dpp, new_state, v, true)
        state[1] .= new_state[1]
        state[2] .= new_state[2]
    end
end

"""
    Compute accept probability to insert / delete u from the state.
"""
function _comp_accept_prob(
    dpp::DeterminantalPointProcess, state::MCMCState, u::Int, insert::Bool
)
    z, L_z_inv = state

    d_u = dpp.L[u, u]
    if any(z)
        b_u = dpp.L[z, u]
        d_u -= dot(b_u, L_z_inv[z, z] * b_u)
    end

    return insert ? min(1.0, d_u) : min(1.0, 1.0 / d_u)
end

"""
    Compute accept probability to swap u and v.
"""
function _comp_accept_prob(dpp::DeterminantalPointProcess, state::MCMCState, u::Int, v::Int)
    z, L_z_inv = state

    d_u, d_v = dpp.L[u, u], dpp.L[v, v]
    if any(z)
        b_u, b_v = dpp.L[z, u], dpp.L[z, v]
        d_u -= dot(b_u, L_z_inv[z, z] * b_u)
        d_v -= dot(b_v, L_z_inv[z, z] * b_v)
    end

    return min(1.0, d_v / d_u)
end

"""
    Update the state after u is inserted / deleted.
"""
function _update_mcmc_state!(
    dpp::DeterminantalPointProcess, state::MCMCState, u::Int, insert::Bool
)
    z, L_z_inv = state

    if insert
        d_u = dpp.L[u, u]
        if any(z)
            b_u = dpp.L[z, u]
            x_u = L_z_inv[z, z] * b_u
            d_u -= dot(b_u, x_u)

            L_z_inv[z, z] += (x_u * x_u') / d_u
            L_z_inv[z, u] = L_z_inv[u, z] = -x_u / d_u
        end

        L_z_inv[u, u] = 1.0 / d_u
        z[u] = true
    else  # delete
        z[u] = false

        e = L_z_inv[z, u]
        f = L_z_inv[u, u]

        L_z_inv[z, z] -= (e * e') / f
    end
end

"""
    Update the state after u is inserted / deleted.
"""
function _update_mcmc_state(
    dpp::DeterminantalPointProcess, state::MCMCState, u::Int, insert::Bool
)
    new_state = deepcopy(state)
    _update_mcmc_state!(dpp, new_state, u, insert)
    return new_state
end

"""
    MCMC sampling from a DPP [2].

TODO: Add support for running MCMC in parallel, similar as rand.
          Make sure parallelization produces unbiased and consistent samples.
"""
function randmcmc(
    rng::AbstractRNG,
    dpp::DeterminantalPointProcess{T},
    N::Int;
    init_state=nothing,
    return_final_state::Bool=false,
    mix_eps::Real=1e-1,
    mixing_steps::Int=ceil(Int, dpp.size * log(dpp.size / mix_eps)),
    steps_between_samples::Int=mixing_steps,
) where {T}
    # initialize the Markov chain
    state = init_state
    if state === nothing
        L_z_inv = Array{T}(undef, size(dpp.L))
        z = bitrand(rng, dpp.size)  # TODO: improve initialization (?)
        if any(z)
            L_z_inv[z, z] = pinv(dpp.L[z, z])
        end
        state = (z, L_z_inv)
    end

    # sanity check
    @assert state isa MCMCState

    # mix the Markov chain
    for _ in 1:mixing_steps
        _do_mcmc_step!(rng, dpp, state)
    end

    Y = []
    for _ in 1:N
        push!(Y, findall(state[1]))
        for t in 1:steps_between_samples
            _do_mcmc_step!(rng, dpp, state)
        end
    end

    return return_final_state ? (Y, state) : Y
end

function randmcmc(
    dpp::DeterminantalPointProcess,
    N::Int;
    init_state=nothing,
    return_final_state::Bool=false,
    mix_eps::Real=1e-1,
    mixing_steps::Int=ceil(Int, dpp.size * log(dpp.size / mix_eps)),
    steps_between_samples::Int=mixing_steps,
)
    return randmcmc(
        GLOBAL_RNG,
        dpp,
        N;
        init_state=init_state,
        return_final_state=return_final_state,
        mix_eps=mix_eps,
        mixing_steps=mixing_steps,
        steps_between_samples=steps_between_samples,
    )
end

"""
    randmcmc([rng::AbstractRNG], kdpp::kDPP, N::Int; kwargs...)

MCMC sampling from a k-DPP [2].

## Arguments

- `rng` : Random number generator (by default Random.GLOBAL_RNG is used)
- `kdpp` : k-DeterminantalPointProcess
- `N` : Number of samples

## Keyword Arguments

- `init_state`
- `return_final_state`
- `mix_eps`
- `mixing_steps`
- `steps_between_samples`

TODO:
- Add support for running MCMC in parallel, similar as rand.
- Make sure parallelization produces unbiased and consistent samples.
"""
function randmcmc(
    rng::AbstractRNG,
    kdpp::kDPP{T},
    N::Int;
    init_state=nothing,
    return_final_state::Bool=false,
    mix_eps::Real=1e-1,
    mixing_steps::Int=ceil(Int, kdpp.k * log(kdpp.k / mix_eps)),
    steps_between_samples::Int=mixing_steps,
) where {T}

    # initialize the Markov chain
    state = init_state
    if state === nothing
        L_z_inv = Array{T}(undef, size(kdpp.dpp.L))
        z = falses(kdpp.dpp.size)  # TODO: improve initialization (?)
        z[1:(kdpp.k)] .= true
        if any(z)
            L_z_inv[z, z] = pinv(kdpp.dpp.L[z, z])
        end
        state = (z, L_z_inv)
    end

    # sanity check
    @assert typeof(state) == MCMCState
    @assert sum(state[1]) == kdpp.k

    # mix the Markov chain
    for _ in 1:mixing_steps
        _do_mcmc_k_step!(rng, kdpp, state)
    end

    Y = []
    for _ in 1:N
        push!(Y, findall(state[1]))
        for t in 1:steps_between_samples
            _do_mcmc_k_step!(rng, kdpp, state)
        end
    end

    return return_final_state ? (Y, state) : Y
end

function randmcmc(
    kdpp::kDPP,
    N::Int;
    init_state=nothing,
    return_final_state::Bool=false,
    mix_eps::Real=1e-1,
    mixing_steps::Int=ceil(Int, kdpp.k * log(kdpp.k / mix_eps)),
    steps_between_samples::Int=mixing_steps,
)
    return randmcmc(
        GLOBAL_RNG,
        kdpp,
        N;
        init_state=init_state,
        return_final_state=return_final_state,
        mix_eps=mix_eps,
        mixing_steps=mixing_steps,
        steps_between_samples=steps_between_samples,
    )
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 28558eb8-39be-11ec-16d0-6dbc33614616
md"""
# Deploying Genie apps with Heroku Buildpacks

This tutorial shows how to host a Julia/Genie app using a Heroku Buildpack.

## Prerequisites

This guide assumes you have a Heroku account and are signed into the Heroku CLI. [Information on how to setup the Heroku CLI is available here](https://devcenter.heroku.com/articles/heroku-cli).

## The application

In order to try the deployment, you will need a sample application. Either pick one of yours or clone this sample one, as indicated next.

### All Steps (in easy copy-paste format):

Customize your `HEROKU_APP_NAME` to something unique:
"""

# ╔═╡ 0afc5e62-39b1-40e8-a3aa-e3098a9f4904
md"""
```sh
HEROKU_APP_NAME=my-app-name
```
"""

# ╔═╡ fee6e860-24cd-4e03-8f8c-26b82f5e8e2b
md"""
This must be unique among all Heroku projects, and is part of the url where your project is hosted (e.g. https://my-app-name.herokuapp.com/).

If the name is not unique, you will see this error at the `heroku create` step.

```sh
Creating ⬢ my-app-name... !
 ▸    Name my-app-name is already taken
```

#### Clone an example project

```sh
git clone https://github.com/milesfrain/GenieOnHeroku.git
cd GenieOnHeroku
```

You may also point to your own project, but it must be a git repo.

A `Procfile` in the root contains the launch command to load your app.
The contents of the `Procfile` for this project is this single line:

```sh
web: julia --project src/app.jl $PORT
```

You may edit the `Procfile` to point to your own project's launch script. (for example `src/my_app_launch_file.jl` instead of `src/app.jl`),
but be sure to take into account the dynamically changing `$PORT` environment variable which is set by Heroku.

If you're deploying a standard Genie application built with `Genie.newapp`, the launch script will be `bin/server`. Genie will automatically pick the `$PORT` number from the environment.

#### Create a Heroku project

```sh
heroku create $HEROKU_APP_NAME --buildpack https://github.com/Optomatica/heroku-buildpack-julia.git
```

This creates a project on the Heroku platform, which includes a separate git repository.

This `heroku` repository is added to the list of tracked repositories and can be observed with `git remote -v`.

```sh
heroku  https://git.heroku.com/my-app-name.git (fetch)
heroku  https://git.heroku.com/my-app-name.git (push)
origin  https://github.com/milesfrain/GenieOnHeroku.git (fetch)
origin  https://github.com/milesfrain/GenieOnHeroku.git (push)
```

We are using a buildpack for Julia. This runs many of the common deployment operations required for Julia projects.
It relies on the directory layout found in the example project, with `Project.toml`, `Manifest.toml` in the root,
and all Julia code in the `src` directory.

#### Deploy your app

```sh
git push heroku master
```

This pushes your current branch of your local repo to the `heroku` remote repo's `master` branch.

Heroku will automatically execute the commands described in the Julia buildpack and Procfile of this latest push.

You must push to the heroku `master` branch to trigger an automated deploy.

#### Open your app's webpage

```sh
heroku open -a $HEROKU_APP_NAME
```

This is a convenience command to open your app's webpage in your browser.

The webpage is: `https://$HEROKU_APP_NAME.herokuapp.com/`

For example: <https://my-app-name.herokuapp.com/>

#### View app logs

```sh
heroku logs -tail -a $HEROKU_APP_NAME
```

This is another convenience command to launch a log viewer that remains open to show the latest status of your app.

The `println` statements from Julia will also appear here.

Exit this viewer with `Ctrl-C`.

Logs can also be viewed from the Heroku web dashboard.
For example: <https://dashboard.heroku.com/apps/my-app-name/logs>

### Deploy app updates changes

To deploy any changes made to your app, simply commit those changes locally, and re-push to heroku.

```sh
<make changes>
git commit -am "my commit message"
git push heroku master
```
"""

# ╔═╡ Cell order:
# ╟─28558eb8-39be-11ec-16d0-6dbc33614616
# ╟─0afc5e62-39b1-40e8-a3aa-e3098a9f4904
# ╟─fee6e860-24cd-4e03-8f8c-26b82f5e8e2b

ext: jl


text:
using Random
using Images
using FileIO
using DynamicPolynomials

@polyvar x
@polyvar y
@polyvar z
@polyvar w

global a

function progs(ndims, order)
    m = 1
    dims = ["x", "y", "z", "w"]

    progs = Dict{String, Vector{String}}()
    for i in 1:ndims
        d = dims[i]
        progs[d] = Vector{String}()
        push!(progs[d], "a[$m]")
        m += 1
        for i1 in 1:ndims
            push!(progs[d], "a[$m] * $(dims[i1])")
            m += 1
            for i2 in i1:ndims
                push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2])")
                m += 1
                if order > 2
                    for i3 in i2:ndims
                        push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3])")
                        m += 1
                        if order > 3
                            for i4 in i3:ndims
                                push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3]) * $(dims[i4])")
                                m += 1
                                if order > 4
                                    for i5 in i4:ndims
                                        push!(progs[d], "a[$m] * $(dims[i1]) * $(dims[i2]) * $(dims[i3]) * $(dims[i4]) * $(dims[i5])")
                                        m += 1
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    progs, m-1
end

function polys(es::Dict{String, Expr})
    ps = Dict{String, DynamicPolynomials.Polynomial{true,Float64}}()

    for d in keys(es)
        ps[d] = eval(es[d])
    end
    ps
end

function exprs(ps::Dict{String, Vector{String}})

    es = Dict{String, Expr}()
    for d in keys(ps)
        es[d] = Meta.parse(join(ps[d], " + "))
    end

    es
end

function write_raw(iters, w, h)
    xrange = (minimum(iters["x"]), maximum(iters["x"]))
    xscale = (w-1) / (xrange[2]-xrange[1])
    yrange = (minimum(iters["y"]), maximum(iters["y"]))
    yscale = (h-1) / (yrange[2]-yrange[1])

    pixels = zeros(Int64, w, h)

    for i in 1:length(iters["x"])
        fx = iters["x"][i]
        fy = iters["y"][i]
        if abs(fx) == Inf || isnan(fx) || abs(fy) == Inf || isnan(fy)
            break
        end
        x = floor(Int64, 1+(fx - xrange[1]) * xscale)
        y = floor(Int64, 1+(fy - yrange[1]) * yscale)

        pixels[x, y] += 1
    end

    gmax = 1 / maximum(pixels)

    gs =  (pixels .* gmax) |> colorview(Gray)

    save("t.jpg", gs)
    gs
end

function string2args(s::String)
    [(Int(c)-77)/10.0 for c in s]
end

function iterate(ps, states, ns)

    n = 1 + ns[end]-ns[1]

    iters = Dict{String, Vector{Float64}}()
    for d in keys(ps)
        iters[d] = zeros(Float64, n)
    end

    for i in 1:n
         for d in keys(iters)
            iters[d][i] = ps[d](x=>states["x"], y=>states["y"], z=>states["z"], w=>states["w"])
        end
         for d in keys(iters)
            states[d] = iters[d][i]
        end

    end

    iters
end

function lyup(vs)
    L = 0.
    for i in 2:length(vs)
        L += abs(vs[i-1] - vs[i])
    end
    L /= length(vs)
    0.721347 * log(L)
end

function valid(iters)
   for d in keys(iters)
        ld = lyup(iters[d])
        if isnan(ld) || ld == Inf
            return false
        end

        if 1 < abs(ld) < 2
            return false
        end
    end
    return true
end

while true
      global a
    pgs, as = progs(2,2)
    a = rand(as)
    es = exprs(pgs)
    pl = polys(es)

    states = Dict{String, Float64}()
    states["x"] = 0.05
    states["y"] = 0.05
    states["z"] = 0.05
    states["w"] = 0.05

    iters = iterate(pl, states, 1:100)

    if valid(iters)

        states["x"] = 0.05
        states["y"] = 0.05
        states["z"] = 0.05
        states["w"] = 0.05

        write_raw(iterate(pl, states, 1:50000), 800, 600)
        break
    end
end


text:
"""
    classnames_coarse(; [dir]) -> Vector{String}

Return the 20 names for the CIFAR100 superclasses as a vector of
strings. Note that these strings are read from the actual
resource file.

$(download_docstring("CIFAR100", DEPNAME))
"""
function classnames_coarse(; dir = nothing)
    file_path = datafile(DEPNAME, COARSE_FILENAME, dir)
    readlines(file_path)
end

"""
    classnames_fine(; [dir]) -> Vector{String}

Return the 100 names for the CIFAR100 classes as a vector of
strings. Note that these strings are read from the actual
resource file.

$(download_docstring("CIFAR100", DEPNAME))
"""
function classnames_fine(; dir = nothing)
    file_path = datafile(DEPNAME, FINE_FILENAME, dir)
    readlines(file_path)
end

"""
    traintensor([T = N0f8], [indices]; [dir]) -> Array{T}

Return the CIFAR-100 **training** images corresponding to the
given `indices` as a multi-dimensional array of eltype `T`. If
the corresponding labels are required as well, it is recommended
to use [`CIFAR100.traindata`](@ref) instead.

The image(s) is/are returned in the native horizontal-major
memory layout as a single numeric array. If `T <: Integer`, then
all values will be within `0` and `255`, otherwise the values are
scaled to be between `0` and `1`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 4D array (i.e. a `Array{T,4}`) in
WHCN format (width, height, #channels, #images).
For integer `indices` instead, a 3D array in WHC format is returned.

```julia-repl
julia> CIFAR100.traintensor() # load all training images
32×32×3×50000 Array{N0f8,4}:
[...]

julia> CIFAR100.traintensor(Float32, 1:3) # first three images as Float32
32×32×3×3 Array{Float32,4}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Array{T,3}`.

```julia-repl
julia> CIFAR100.traintensor(1) # load first training image
32×32×3 Array{N0f8,3}:
[...]
```

You can use the utility function
[`convert2image`](@ref) to convert an CIFAR-100 array into a
vertical-major Julia image with the appropriate `RGB` eltype.

```julia-repl
julia> CIFAR100.convert2image(CIFAR100.traintensor(1)) # convert to column-major colorant array
32×32 Array{RGB{N0f8},2}:
[...]
```

$(download_docstring("CIFAR100", DEPNAME))
"""
function traintensor(args...; dir = nothing)
    traindata(args...; dir = dir)[1]
end

# needed for type inference for some reason
function traintensor(::Type{T}, args...; dir = nothing) where T
    traindata(T, args...; dir = dir)[1]
end

"""
    testtensor([T = N0f8], [indices]; [dir]) -> Array{T}

Return the CIFAR-100 **test** images corresponding to the given
`indices` as a multi-dimensional array of eltype `T`. If the
corresponding labels are required as well, it is recommended to
use [`CIFAR100.testdata`](@ref) instead.

The image(s) is/are returned in the native horizontal-major
memory layout as a single numeric array. If `T <: Integer`, then
all values will be within `0` and `255`, otherwise the values are
scaled to be between `0` and `1`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 4D array (i.e. a `Array{T,4}`), in which
the first dimension corresponds to the pixel *columns* (x) of the
image, the second dimension to the pixel *rows* (y) of the
image, the third dimension the RGB color channels, and the fourth
dimension denotes the index of the image.

```julia-repl
julia> CIFAR100.testtensor() # load all training images
32×32×3×10000 Array{N0f8,4}:
[...]

julia> CIFAR100.testtensor(Float32, 1:3) # first three images as Float32
32×32×3×3 Array{Float32,4}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Array{T,3}`.

```julia-repl
julia> CIFAR100.testtensor(1) # load first training image
32×32×3 Array{N0f8,3}:
[...]
```

You can use the utility function
[`convert2image`](@ref) to convert an CIFAR-100 array into a
vertical-major Julia image with the appropriate `RGB` eltype.

```julia-repl
julia> CIFAR100.convert2image(CIFAR100.testtensor(1)) # convert to column-major colorant array
32×32 Array{RGB{N0f8},2}:
[...]
```

$(download_docstring("CIFAR100", DEPNAME))
"""
function testtensor(args...; dir = nothing)
    testdata(args...; dir = dir)[1]
end

function testtensor(::Type{T}, args...; dir = nothing) where T
    testdata(T, args...; dir = dir)[1]
end

"""
    trainlabels([indices]; [dir]) -> Yc, Yf

Return the CIFAR-100 **trainset** labels (coarse and fine)
corresponding to the given `indices` as a tuple of two `Int` or
two `Vector{Int}`. The variables returned are the coarse label(s)
(`Yc`) and the fine label(s) (`Yf`) respectively.

```julia
Yc, Yf = CIFAR100.trainlabels(); # full training set
```

The values of the labels denote the zero-based class-index that
they represent (see [`CIFAR100.classnames_coarse`](@ref) and
[`CIFAR100.classnames_fine`](@ref) for the corresponding names).
If `indices` is omitted, all labels are returned.

```julia-repl
julia> Yc, Yf = CIFAR100.trainlabels(1:3) # first three labels
([11, 15, 4], [19, 29, 0])

julia> yc, yf = CIFAR100.trainlabels(1) # first label
(11, 19)

julia> CIFAR100.classnames_coarse()[yc + 1] # corresponding superclass name
"large_omnivores_and_herbivores"

julia> CIFAR100.classnames_fine()[yf + 1] # corresponding class name
"cattle"
```

$(download_docstring("CIFAR100", DEPNAME))
"""
function trainlabels(args...; dir = nothing)
    _, Yc, Yf = traindata(args...; dir = dir)
    Yc, Yf
end

"""
    testlabels([indices]; [dir]) -> Yc, Yf

Return the CIFAR-100 **testset** labels (coarse and fine)
corresponding to the given `indices` as a tuple of two `Int` or
two `Vector{Int}`. The variables returned are the coarse label(s)
(`Yc`) and the fine label(s) (`Yf`) respectively.

```julia
Yc, Yf = CIFAR100.testlabels(); # full training set
```

The values of the labels denote the zero-based class-index that
they represent (see [`CIFAR100.classnames_coarse`](@ref) and
[`CIFAR100.classnames_fine`](@ref) for the corresponding names).
If `indices` is omitted, all labels are returned.

```julia-repl
julia> Yc, Yf = CIFAR100.testlabels(1:3) # first three labels
([10, 10, 0], [49, 33, 72])

julia> yc, yf = CIFAR100.testlabels(1) # first label
(10, 49)

julia> CIFAR100.classnames_coarse()[yc + 1] # corresponding superclass name
"large_natural_outdoor_scenes"

julia> CIFAR100.classnames_fine()[yf + 1] # corresponding class name
"mountain"
```

$(download_docstring("CIFAR100", DEPNAME))
"""
function testlabels(args...; dir = nothing)
    _, Yc, Yf = testdata(args...; dir = dir)
    Yc, Yf
end

"""
    traindata([T = N0f8], [indices]; [dir]) -> X, Yc, Yf

Returns the CIFAR-100 **trainset** corresponding to the given
`indices` as a three-element tuple. If `indices` is omitted the
full trainingset is returned. The first element of the three
return values (`X`) will be the images as a multi-dimensional
array, the second element (`Yc`) the corresponding coarse labels
as integers, and the third element (`Yf`) the fine labels
respectively.

The image(s) is/are returned in the native horizontal-major
memory layout as a single numeric array of eltype `T`. If `T <:
Integer`, then all values will be within `0` and `255`, otherwise
the values are scaled to be between `0` and `1`. The integer
values of the labels correspond 1-to-1 the digit that they
represent.

```julia
X, Yc, Yf = CIFAR100.traindata() # full datatset
X, Yc, Yf = CIFAR100.traindata(dir="./CIFAR100") # custom folder
x, yc, yf = CIFAR100.traindata(2) # only second observation
```

$(download_docstring("CIFAR100", DEPNAME))

Take a look at [`CIFAR100.traintensor`](@ref) and
[`CIFAR100.trainlabels`](@ref) for more information.
"""
function traindata(args...; dir = nothing)
    traindata(N0f8, args...; dir = dir)
end

function traindata(::Type{T}; dir = nothing) where T
    file_path = datafile(DEPNAME, TRAINSET_FILENAME, dir)
    # simply read the complete content of the testset file
    images, labels_c, labels_f = Reader.readdata(file_path, TRAINSET_SIZE)

    # optionally transform the image array before returning
    bytes_to_type(T, images), labels_c, labels_f
end

function traindata(::Type{T}, index::Integer; dir = nothing) where T
    @assert 1 <= index <= TRAINSET_SIZE "parameter \"index\" ($index) not in 1:$(TRAINSET_SIZE)"
    file_path = datafile(DEPNAME, TRAINSET_FILENAME, dir)
    # read the single image+labelc+labelf corresponding to "index"
    image, label_c, label_f = Reader.readdata(file_path, TRAINSET_SIZE, index)
    # optionally transform the image array before returning
    bytes_to_type(T, image), label_c, label_f
end

function traindata(::Type{T}, indices::AbstractVector; dir = nothing) where T
    mi, ma = extrema(indices)
    @assert mi >= 1 && ma <= TRAINSET_SIZE "not all elements in parameter \"indices\" are in 1:$(TRAINSET_SIZE)"
    # preallocate a buffer we will reuse for reading individual
    # images. "buffer" is written to length(indices) times
    buffer = Array{UInt8,3}(undef, Reader.NROW, Reader.NCOL, Reader.NCHAN)
    # we know the types and dimensions of the return values,
    # so we can preallocate them
    images = Array{UInt8,4}(undef, Reader.NROW, Reader.NCOL, Reader.NCHAN, length(indices))
    labels_c = Array{Int,1}(undef, length(indices))
    labels_f = Array{Int,1}(undef, length(indices))
    file_path = datafile(DEPNAME, TRAINSET_FILENAME, dir)
    open(file_path, "r") do io
        # iterate over the given "indices"
        @inbounds for (i, index) in enumerate(indices)
            # read the corresponding image into "buffer" and
            # the corresponding coarse label into "yc" and
            # the corresponding fine label into "yf"
            _, yc, yf = Reader.readdata!(buffer, io, index)
            # write the image into the appropriate position
            # of our preallocated "images" array.
            copyto!(view(images,:,:,:,i), buffer)
            # same with labels
            labels_c[i] = yc
            labels_f[i] = yf
        end
    end

    # optionally transform the image array before returning
    bytes_to_type(T, images::Array{UInt8,4}), labels_c::Vector{Int}, labels_f::Vector{Int}
end

"""
    testdata([T = N0f8], [indices]; [dir]) -> X, Yc, Yf

Returns the CIFAR-100 **testset** corresponding to the given
`indices` as a three-element tuple. If `indices` is omitted the
full testset is returned. The first element of the three return
values (`X`) will be the images as a multi-dimensional array, the
second element (`Yc`) the corresponding coarse labels as
integers, and the third element (`Yf`) the fine labels
respectively.

The image(s) is/are returned in the native horizontal-major
memory layout as a single numeric array of eltype `T`. If `T <:
Integer`, then all values will be within `0` and `255`, otherwise
the values are scaled to be between `0` and `1`. The integer
values of the labels correspond 1-to-1 the digit that they
represent.

```julia
X, Yc, Yf = CIFAR100.testdata() # full datatset
X, Yc, Yf = CIFAR100.testdata(dir="./CIFAR100") # custom folder
x, yc, yf = CIFAR100.testdata(2) # only second observation
```

$(download_docstring("CIFAR100", DEPNAME))

Take a look at [`CIFAR100.testtensor`](@ref) and
[`CIFAR100.testlabels`](@ref) for more information.
"""
function testdata(args...; dir = nothing)
    testdata(N0f8, args...; dir = dir)
end

function testdata(::Type{T}; dir = nothing) where T
    file_path = datafile(DEPNAME, TESTSET_FILENAME, dir)
    # simply read the complete content of the testset file
    images, labels_c, labels_f = Reader.readdata(file_path, TESTSET_SIZE)
    # optionally transform the image array before returning
    bytes_to_type(T, images), labels_c, labels_f
end

function testdata(::Type{T}, index::Integer; dir = nothing) where T
    @assert 1 <= index <= TESTSET_SIZE "parameter \"index\" ($index) not in 1:$(TESTSET_SIZE)"
    file_path = datafile(DEPNAME, TESTSET_FILENAME, dir)
    # read the single image+labelc+labelf corresponding to "index"
    image, label_c, label_f = Reader.readdata(file_path, TESTSET_SIZE, index)
    # optionally transform the image array before returning
    bytes_to_type(T, image), label_c, label_f
end

function testdata(::Type{T}, indices::AbstractVector; dir = nothing) where T
    mi, ma = extrema(indices)
    @assert mi >= 1 && ma <= TESTSET_SIZE "not all elements in parameter \"indices\" are in 1:$(TESTSET_SIZE)"
    # preallocate a buffer we will reuse for reading individual
    # images. "buffer" is written to length(indices) times
    buffer = Array{UInt8,3}(undef, Reader.NROW, Reader.NCOL, Reader.NCHAN)
    # we know the types and dimensions of the return values,
    # so we can preallocate them
    images = Array{UInt8,4}(undef, Reader.NROW, Reader.NCOL, Reader.NCHAN, length(indices))
    labels_c = Array{Int,1}(undef, length(indices))
    labels_f = Array{Int,1}(undef, length(indices))
    file_path = datafile(DEPNAME, TESTSET_FILENAME, dir)
    open(file_path, "r") do io
        # iterate over the given "indices"
        @inbounds for (i, index) in enumerate(indices)
            # read the corresponding image into "buffer" and
            # the corresponding coarse label into "yc" and
            # the corresponding fine label into "yf"
            _, yc, yf = Reader.readdata!(buffer, io, index)
            # write the image into the appropriate position
            # of our preallocated "images" array.
            copyto!(view(images,:,:,:,i), buffer)
            # same with labels
            labels_c[i] = yc
            labels_f[i] = yf
        end
    end
    # optionally transform the image array before returning
    bytes_to_type(T, images::Array{UInt8,4}), labels_c::Vector{Int}, labels_f::Vector{Int}
end


text:
export  coeffs,
        nw,
        dim,
        deg,
        multi2uni,
        getentry,
        issymmetric,
        integrate

dim(op::OrthoPoly)::Int64 = op.deg + 1
dim(opq::OrthoPolyQ) = dim(opq.op)
dim(mop::MultiOrthoPoly) = mop.dim

deg(op::OrthoPoly) = op.deg
deg(opq::OrthoPolyQ) = deg(opq.op)
deg(mop::MultiOrthoPoly) = mop.deg

"""
```
nw(q::Quad)
nw(opq::OrthoPolyQ)
nw(opq::Vector{OrthoPolyQ})
nw(mOP::MultiOrthoPoly)
```
returns nodes and weights in matrix form
"""
function nw(q::Quad)
    [q.nodes q.weights]
end
nw(opq::OrthoPolyQ) = nw(opq.quad)
function nw(qq::Vector{Quad})
    n = [ q.nodes for q in qq]
    w = [ q.weights for q in qq]
    return n,w
end
function nw(opq::Vector{OrthoPolyQ})
    q = [ p.quad for p in opq ]
    nw(q)
end

nw(mOP::MultiOrthoPoly) = nw(mOP.uni)

"""
```
coeffs(op::OrthoPoly)
coeffs(opq::OrthoPolyQ)
coeffs(op::Vector{OrthoPoly})
coeffs(opq::Vector{OrthoPolyQ})
coeffs(mop::MultiOrthoPoly)
```
returns recurrence coefficients of in matrix form
"""
function coeffs(op::OrthoPoly)
    [op.α op.β]
end
coeffs(opq::OrthoPolyQ) = coeffs(opq.op)

function coeffs(op::Vector{OrthoPoly})
    a = [ p.α for p in op]
    b = [ p.β for p in op]
    return a,b
end
function coeffs(opq::Vector{OrthoPolyQ})
    a = [ p.op.α for p in opq]
    b = [ p.op.β for p in opq]
    return a,b
end
coeffs(mop::MultiOrthoPoly) = coeffs(mop.uni)

"""
```
integrate(f::Function,nodes::Vector{Float64},weights::Vector{Float64})
integrate(f::Function,q::Quad)
integrate(f::Function,opq::OrthogonalPolyQ)
```
integrate function `f` using quadrature rule specified via `nodes`, `weights`.
For example ``\\int_0^1 6x^5 = 1`` can be solved as follows:

```@repl
julia> opq = OrthoPolyQ("uniform01",3)
julia> integrate(x->6x^5,opq)
1.0000000000000002
```

!!! note
- function ``f`` is assumed to return a scalar.
- interval of integration is "hidden" in `nodes`.
"""
function integrate(f::Function,nodes::Vector{Float64},weights::Vector{Float64})::Float64
    dot(weights,f.(nodes))
end
integrate(f::Function,q::Quad)::Float64 = integrate(f,q.nodes,q.weights)
integrate(f::Function,opq::OrthoPolyQ)::Float64 = integrate(f,opq.quad)

"""
```
issymmetric(m::Measure)::Bool
issymmetric(op::OrthoPoly)::Bool
issymmetric(q::Quad)::Bool
issymmetric(opq::OrthoPolyQ)::Bool
```
Is the measure symmetric (around any point in the domain)?
"""
function issymmetric(m::Measure)::Bool
    m.symmetric
end
issymmetric(op::OrthoPoly)::Bool = issymmetric(op.meas)
issymmetric(q::Quad)::Bool = issymmetric(q.meas)
function issymmetric(opq::OrthoPolyQ)::Bool
    @assert issymmetric(opq.op)==issymmetric(opq.quad) "inconsistent symmetries"
    issymmetric(opq.op)
end

function multi2uni(a::Vector{Int64},ind::Matrix{Int64})
    @assert minimum(a)>=0 "no negative degrees allowed"
    l,p = size(ind) # p-variate basis
    m = length(a) # dimension of scalar product
    l=l-1 # (l+1)-dimensional basis
    @assert maximum(a)<=l "not enough elements in multi-index"
    A = zeros(Int64,p,m)
    for (i,a_) in enumerate(a)
        A[:,i] = ind[a_+1,:]
    end
    return A
end

function getentry(a::Vector{Int64},T::SparseVector{Float64,Int64},ind::Matrix{Int64},dim::Int64)
    m::Int64 = length(a)
    l::Int64 = size(ind,1)-1
    @assert minimum(a)>=0 "no negative degrees allowed"
    @assert maximum(a)<=l "not enough elements in multi-index (requested: $(maximum(a)), max: $l)"
    @assert m==dim "length $m of provided index $a is inconsistent with dimension $(dim) of multi-index "
    # a .+= 1
    sort!(a)
    sp_ind::Int64 = 1 + reduce(+,[idx*l^(m-i) for (i,idx) in enumerate(a)])
    return T[sp_ind]::Float64
end


text:
"""
# Get the source of a specific time series.

    get_source(source_id)

## Arguments

    - `source_id`: an ID, e.g. 1 for "Board of Governors of the Federal Reserve System (US)"

## Examples

```jldoctest
julia> get_source(1)
julia> get_source(3)
julia> get_source(4)
```
"""
function get_source(source_id::Int)

    url = "https://api.stlouisfed.org/fred/source"
    parameters = Dict("api_key" => ENV["API_FRED"],
                      "file_type" => "json",
                      "source_id" => source_id)

    response = HTTP.request("GET", url; query = parameters)
    results = JSON.parse(String(response.body))["sources"][1]

    return(results)
end


text:
# ============================================================
# Training & Optimization Methods
# ============================================================

"""
    sci_train(
        ann::T,
        data,
        θ = initial_params(ann),
        opt = ADAM(),
        adtype = GalacticOptim.AutoZygote(),
        args...;
        maxiters = 200::Integer,
        kwargs...,
    ) where {T<:DiffEqFlux.FastLayer}

    sci_train(loss::Function, θ, opt, adtype = GalacticOptim.AutoZygote(), args...; kwargs...)

Scientific machine learning training function

- @arg ann: neural network model
- @arg data: tuple (X, Y) of dataset
- @arg θ: parameters of neural network
- @arg opt: optimizer
- @arg adtype: automatical differentiation type
- @arg args: rest arguments
- @arg device: cpu / gpu
- @arg maxiters: maximal iteration number
- @arg kwargs: keyword arguments

"""
function sci_train(
    ann::T,
    data,
    θ = initial_params(ann),
    opt = ADAM(),
    adtype = GalacticOptim.AutoZygote(),
    args...;
    device = cpu,
    maxiters = 200::Integer,
    kwargs...,
) where {T<:DiffEqFlux.FastLayer}
    data = data |> device
    θ = θ |> device
    L = size(data[1], 2)
    loss(p) = sum(abs2, ann(data[1], p) - data[2]) / L

    cb = function (p, l)
        println("loss: $(loss(p))")
        return false
    end

    f = GalacticOptim.OptimizationFunction((x, p) -> loss(x), adtype)
    fi = GalacticOptim.instantiate_function(f, θ, adtype, nothing)
    prob = GalacticOptim.OptimizationProblem(fi, θ; kwargs...)

    return GalacticOptim.solve(prob, opt, args...; cb = Flux.throttle(cb, 1), maxiters = maxiters, kwargs...)
end

function sci_train(loss::Function, θ, opt = ADAM(), adtype = GalacticOptim.AutoZygote(), args...; maxiters = 200::Integer, kwargs...)
    f = GalacticOptim.OptimizationFunction((x, p) -> loss(x), adtype)
    fi = GalacticOptim.instantiate_function(f, θ, adtype, nothing)
    prob = GalacticOptim.OptimizationProblem(fi, θ; kwargs...)
    return GalacticOptim.solve(prob, opt, args...; maxiters = maxiters, kwargs...)
end

"""
    sci_train!(ann, data::Tuple, opt = ADAM(); device = cpu, epoch = 1, batch = 1)
    sci_train!(ann, dl::Flux.Data.DataLoader, opt = ADAM(); device = cpu, epoch = 1)

Scientific machine learning training function

- @arg ann: neural network model
- @arg data: tuple (X, Y) of dataset
- @arg opt: optimizer
- @arg epoch: epoch number
- @arg batch: batch size
- @arg device: cpu / gpu

"""
function sci_train!(ann, data::Tuple, opt = ADAM(); device = cpu, epoch = 1, batch = 1)
    X, Y = data |> device
    L = size(X, 2)
    data = Flux.Data.DataLoader((X, Y), batchsize = batch, shuffle = true) |> device

    ann = device(ann)
    ps = Flux.params(ann)
    loss(x, y) = sum(abs2, ann(x) - y) / L
    cb = () -> println("loss: $(loss(X, Y))")

    Flux.@epochs epoch Flux.train!(loss, ps, data, opt, cb = Flux.throttle(cb, 1))

    return nothing
end

function sci_train!(ann, dl::Flux.Data.DataLoader, opt = ADAM(); device = cpu, epoch = 1)
    X, Y = dl.data |> device
    L = size(X, 2)
    dl = dl |> device

    ann = device(ann)
    ps = Flux.params(ann)
    loss(x, y) = sum(abs2, ann(x) - y) / L
    cb = () -> println("loss: $(loss(X, Y))")

    Flux.@epochs epoch Flux.train!(loss, ps, dl, opt, cb = Flux.throttle(cb, 1))

    return nothing
end

# ------------------------------------------------------------
# TensorFlow
# ------------------------------------------------------------
function sci_train!(ann::PyObject, data::Tuple; device = cpu, split = 0.0, epoch = 1, batch = 64, verbose = 1)
    X, Y = data
    ann.fit(X, Y, validation_split=split, epochs=epoch, batch_size=batch, verbose=verbose)

    return nothing
end


text:
using DataFrames, Query, Turing, StatsFuns

using RDatasets

cars = dataset("datasets", "cars")

@model model7_20(dist, speed) = begin
    σ ~ Exponential(1)
    α ~ Normal(0, 100)
    β ~ Normal(0, 10)
    for i ∈ eachindex(dist)
        dist[i] ~ Normal(α + β * speed[i], σ)
    end
end

dist = cars.Dist
speed = cars.Speed

model = model7_20(dist, speed)

m7_20 = sample(model, NUTS(), 2000, save_state = true)

describe(m7_20)

loglik =
[logprob"dist = dist[i] | speed = speed[i], chain = m7_20" for i ∈ eachindex(dist)]

ll = DataFrame(loglik)

function WAIC(ll::DataFrame)
    n_samples, n_obs = size(ll)
    logsums = [logsumexp(col) for col ∈ eachcol(ll)]
    lppd = logsums .- log(n_samples)
    pWAIC = [var(col) for col ∈ eachcol(ll)]
    sum_lppd = sum(lppd)
    sum_pWAIC =sum(pWAIC)
    waic = -2 * (sum_lppd - sum_pWAIC)
    waic_vec = -2 .* (lppd .- pWAIC)
    se = sqrt(n_obs * var(waic_vec))

    return (waic = waic, se = se, lppd = sum_lppd, pwaic = sum_pWAIC)
end

waic, se, lppd, pWAIC = WAIC(ll)


text:
using CUDA.CURAND

@test CURAND.version() isa VersionNumber

rng = CURAND.default_rng()
Random.seed!(rng)
Random.seed!(rng, nothing)
Random.seed!(rng, 1)
Random.seed!(rng, 1, 0)


text:
using Transformers
using Transformers.Basic
using Transformers.Pretrain
using Transformers.Datasets
using Transformers.Datasets: GLUE
using Transformers.BidirectionalEncoder

using Flux
using Flux: onehotbatch, gradient
import Flux.Optimise: update!
using WordTokenizers

include("./args.jl")

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

const Epoch = 20
const Batch = 4

if args["task"] == "cola"
    const task = GLUE.CoLA()

    markline(sent) = ["[CLS]"; sent; "[SEP]"]
    function preprocess(batch)
        sentence = markline.(wordpiece.(tokenizer.(batch[1])))
        mask = getmask(sentence)
        tok = vocab(sentence)
        segment = fill!(similar(tok), 1)

        label = onehotbatch(batch[2], labels)
        return (tok=tok, segment=segment), label, mask
    end
else
    if args["task"] == "mnli"
        const task = GLUE.MNLI(false)
    elseif args["task"] == "mrpc"
        const task = GLUE.MRPC()
    end
    markline(s1, s2) = ["[CLS]"; s1; "[SEP]"; s2; "[SEP]"]
    function preprocess(batch)
        s1 = wordpiece.(tokenizer.(batch[1]))
        s2 = wordpiece.(tokenizer.(batch[2]))
        sentence = markline.(s1, s2)
        mask = getmask(sentence)
        tok = vocab(sentence)

        segment = fill!(similar(tok), 1)
        for (i, sent) ∈ enumerate(sentence)
            j = findfirst(isequal("[SEP]"), sent)
            if j !== nothing
                @view(segment[j+1:end, i]) .= 2
            end
        end

        label = onehotbatch(batch[3], labels)
        return (tok=tok, segment=segment), label, mask
    end
end

const labels = get_labels(task)

const _bert_model, wordpiece, tokenizer = pretrain"Bert-uncased_L-12_H-768_A-12"
const vocab = Vocabulary(wordpiece)

const hidden_size = size(_bert_model.classifier.pooler.W ,1)
const clf = gpu(Chain(
    Dropout(0.1),
    Dense(hidden_size, length(labels)),
    logsoftmax
))

const bert_model = gpu(
    set_classifier(_bert_model,
                   (
                       pooler = _bert_model.classifier.pooler,
                       clf = clf
                   )
                  )
)

const ps = params(bert_model)
const opt = ADAM(1e-6)

function acc(p, label)
    pred = Flux.onecold(collect(p))
    label = Flux.onecold(collect(label))
    sum(pred .== label) / length(label)
end

function loss(data, label, mask=nothing)
    e = bert_model.embed(data)
    t = bert_model.transformers(e, mask)

    p = bert_model.classifier.clf(
        bert_model.classifier.pooler(
            t[:,1,:]
        )
    )

    l = Basic.logcrossentropy(label, p)
    return l, p
end

function train!()
    global Batch
    global Epoch
    @info "start training: $(args["task"])"
    for e = 1:Epoch
        @info "epoch: $e"
        datas = dataset(Train, task)

        i = 1
        al::Float64 = 0.
        while (batch = get_batch(datas, Batch)) !== nothing
            data, label, mask = todevice(preprocess(batch))
            l, p = loss(data, label, mask)
            # @show l
            a = acc(p, label)
            al += a
            grad = gradient(()->l, ps)
            i+=1
            update!(opt, ps, grad)
            i%16==0 && @show al/i
        end

        test()
    end
end

function test()
    Flux.testmode!(bert_model)
    i = 1
    al::Float64 = 0.
    datas = dataset(Dev, task)
    while (batch = get_batch(datas, Batch)) !== nothing
      data, label, mask = todevice(preprocess(batch))
      _, p = loss(data, label, mask)
      # @show l
      a = acc(p, label)
      al += a
      i+=1
    end
    al /= i
    Flux.testmode!(bert_model, false)
    @show al
end

train!()


text:
using Test
using RxJulia

include("./support.jl")

@testset CustomTestSet "Basics" begin

  @test 1 == 1

  @test [1, 2, 3] == begin
    evts = @rx() do
      [1, 2, 3]
    end
    [evt for evt in evts]
  end

  function plus3(observers, x)
    evt = ValueEvent(x + 3)
    notify!(observers, evt)
  end

  @test [4, 5, 6] == begin
    evts = @rx() do
      [1, 2, 3]
      react(plus3)
    end
    [evt for evt in evts]
  end

  @test [4, 5, 6] == begin
    evts1 = @rx() do
      [1, 2, 3]
    end
    evts2 = @rx(evts1) do
      react(plus3)
    end
    [evt for evt in evts2]
  end

  @test [4, 5, 6] == begin
    evts = @rx() do
      [1, 2, 3]
      react(plus3)
    end
    [evt for evt in evts]
  end

  @test [1] == begin
    evts = @rx() do
      1
    end
    [evt for evt in evts]
  end

  @test [1, 3, 5] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6]
      detect(isodd)
    end
    [evt for evt in evts]
  end

  @test [2, 4, 6] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6]
      ignore(isodd)
    end
    [evt for evt in evts]
  end

  @test [2, 3, 4] == begin
    evts = @rx() do
      [1, 2, 3]
      fmap() do value
        value + 1
      end
    end
    [evt for evt in evts]
  end

  @test [1, 2, 3] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6, 7, 8, 9]
      take(3)
    end
    [evt for evt in evts]
  end

  @test [7, 8, 9] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6, 7, 8, 9]
      keep(3)
    end
    [evt for evt in evts]
  end

  @test [1, 3, 5] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6, 7, 8, 9]
      detect(isodd)
      take(3)
    end
    [evt for evt in evts]
  end

  @test [4, 5, 6] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6]
      drop(3)
    end
    [evt for evt in evts]
  end

  @test [1, 2, 3, 4, 5, 6] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6, 7, 8, 9]
      cut(3)
    end
    [evt for evt in evts]
  end

  @test [] == begin
    evts = @rx() do
      [1, 2, 3, 4, 5, 6]
      cut(9)
    end
    [evt for evt in evts]
  end

  @test [1, 2, 3] == begin
    evts = @rx() do
      [1, 1, 2, 1, 3, 3, 2, 3]
      distinct()
    end
    [evt for evt in evts]
  end

  @test [2, 3, 4, 5] == begin
    evts = @rx() do
      span(2, 5)
    end
    [evt for evt in evts]
  end

  @test Set([1 2 3 4 5 6 7]) == begin
    evts3 = @rx() do
      let evts1 = @rx(() -> [1 2 3])
        evts2 = @rx(() -> [4 5 6 7])
        RxJulia.merge(evts1, evts2)
      end
    end
    Set([evt for evt in evts3])
  end

  @test [(1, 4), (2, 5), (3, 6)] == begin
    evts3 = @rx() do
      let evts1 = @rx(() -> [1 2 3])
        evts2 = @rx(() -> [4 5 6])
        RxJulia.zip(evts1, evts2)
      end
    end
    [evt for evt in evts3]
  end

  @test [6] == begin
    evts = @rx() do
      [1 2 3]
      RxJulia.sum()
    end
    [evt for evt in evts]
  end

  @test [2] == begin
    evts = @rx() do
      [1 2 3]
      RxJulia.average()
    end
    [evt for evt in evts]
  end

  @test [3] == begin
    evts = @rx() do
      [1 2 3]
      RxJulia.max()
    end
    [evt for evt in evts]
  end

  @test [1] == begin
    evts = @rx() do
      [1 2 3]
      RxJulia.min()
    end
    [evt for evt in evts]
  end

end


text:
tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_gwf_model(config)
@unpack network = model

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    row = csv_first_row(model.writer.csv_path)

    @test row.time == DateTime("2000-06-01T00:00:00")
    @test row.Q_av ≈ 0.01620324716944374f0
    @test row.head ≈ 1.8416896245327632f0
end

@testset "first timestep" begin
    sbm = model.vertical

    @test model.clock.iteration == 2
    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
    @test sbm.runoff[1] == 0.0
    @test sbm.soilevap[1] == 0.0
    @test sbm.transpiration[1] ≈ 0.6117526566330049f0
end

# run the second timestep
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "second timestep" begin
    sbm = model.vertical
    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
    @test sbm.runoff[1] == 0.0
    @test sbm.soilevap[1] == 0.0
    @test sbm.transpiration[1] ≈ 1.0122634204681036f0
end

@testset "overland flow" begin
    q = model.lateral.land.q_av
    @test sum(q) ≈ 2.2298616f-7
end

@testset "river domain" begin
    q = model.lateral.river.q_av
    river = model.lateral.river
    @test sum(q) ≈ 0.034846861707851576f0
    @test q[6] ≈ 0.007866587223009643f0
    @test river.volume[6] ≈ 4.474556355925399f0
    @test river.inwater[6] ≈ 0.00036392604840276924f0
    @test q[13] ≈ 0.0005952292000597283f0
    @test q[network.river.order[end]] ≈ 0.008387430132453135f0
end

@testset "groundwater" begin
    gw = model.lateral.subsurface
    @test gw.river.stage[1] ≈ 1.2123636929067039f0
    @test gw.flow.aquifer.head[19] ≈ 1.7999999523162842f0
    @test gw.river.flux[1] ≈ -50.46515313302901f0
    @test gw.drain.flux[1] ≈ 0.0
    @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
end

Wflow.close_files(model)

@testset "no drains" begin
    config.model.drains = false
    delete!(Dict(config.output.lateral.subsurface), "drain")
    model = Wflow.initialize_sbm_gwf_model(config)
    @test collect(keys(model.lateral.subsurface)) == [:flow, :recharge, :river]
    Wflow.close_files(model)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

coords2point(coords) = Point(coords)

coords2multipoint(coords) = Multi(coords2point.(coords))

coords2chain(coords) = Chain(Point.(coords))

coords2multichain(coords) = Multi(coords2chain.(coords))

function coords2poly(coords)
  chains = [Point.(coord) for coord in coords]
  PolyArea(chains[begin], chains[begin+1:end])
end

coords2multipoly(coords) = Multi(coords2poly.(coords))

function geom2meshes(geom)
  gtype  = GI.geotype(geom)
  coords = GI.coordinates(geom)
  if gtype == :Point
    coords2point(coords)
  elseif gtype == :LineString
    coords2chain(coords)
  elseif gtype == :MultiLineString
    coords2multichain(coords)
  elseif gtype == :Polygon
    coords2poly(coords)
  elseif gtype == :MultiPolygon
    coords2multipoly(coords)
  end
end


text:
module Duffing

using OrdinaryDiffEq
using NLsolve
using StaticArrays

function duffing!(dxdt, x, p, t)
	dxdt[1] = x[2]
	dxdt[2] = sin(p.ω*t) - p.c*x[2] - p.k*x[1] - p.β*x[1]^3
	nothing
end

function zeroproblem!(res, x, p)
	c, k, β, ω = p
	prob = ODEProblem(duffing!, x, (0.0, 2π/ω), (c=c, k=k, β=β, ω=ω))
	sol = solve(prob, Tsit5(), dtmax=one(eltype(x))/10)
    res[1:2] .= x[1:2] .- sol[1:2, end]
    nothing
end

function runtest()
	# Generate a brute-force bifurcation diagram
	p = [0.05, 1.0, 0.1, 0.0]
	Ω = range(0.1, stop=4, length=51)
	x = Vector{Vector{Float64}}()
	lastx = [0.0, 0.0]
	for ω in Ω
		p[4] = ω
		sol = nlsolve((res, u) -> zeroproblem!(res, u, p), lastx)
		if converged(sol)
			lastx = sol.zero
			push!(x, lastx)
		else
			@warn "Failed to converge" ω
			break
		end
	end
	(Ω, x)
end

end


text:
"""
    is_radiologic(x) -> Bool

Test to see if `x` is in radiological orientation.
"""
is_radiologic(x) = is_radiologic(spatialorder(x))
function is_radiologic(x::NTuple{3,Symbol})
    return is_left(first(x)) & is_anterior(@inbounds(x[2])) & is_superior(last(x))
end

"""
    is_neurologic(x) -> Bool

Test to see if `x` is in neurological orientation.
"""
is_neurologic(x) = is_neurologic(spatial_order(x))
function is_neurologic(x::NTuple{3,Symbol})
    return is_right(first(x)) & is_anterior(@inbounds(x[2])) & is_superior(last(x))
end

"""
    is_anatomical(x) -> Bool

Returns `true` if `T` represents anatomical data.
"""
@inline function is_anatomical(x)
    if has_sagittaldim(x)
        return true
    elseif has_axialdim(x)
        return true
    elseif has_coronaldim(x)
        return true
    else
        return false
    end
end

"""
    is_electrophysiology(::T) -> Bool

Returns `true` if `T` represents electrophysiology data.
"""
is_electrophysiology(::T) where {T} = is_electrophysiology(T)
is_electrophysiology(::Type{T}) where {T} = has_channel_axis(T) & has_time_axis(T)


text:
export MISSolverType

"""
# Description
    MISSolverType(;
        splitting_type = SlowFastSplitting(),
        fast_model = AtmosAcousticGravityLinearModel,
        mis_method = MIS2,
        fast_method = LSRK54CarpenterKennedy,
        nsubsteps = 50,
    )

This solver type constructs an ODE solver using a generalization
of the split-explicit Runge-Kutta method. Known as the Multirate
Infinitesimal Step (MIS) method, this solver solves ODEs with
the partitioned form:

```math
    \\dot{Q} = f_{fast}(Q, t) + f_{slow}(Q, t)
```

where the right-hand-side functions `f_fast` and `f_slow` denote
fast and slow dynamics respectively, depending on the state `Q`.

# Arguments
- `splitting_type` (DiscreteSplittingType): The type of discrete
    splitting to apply to the right-hand side.
    Default: `SlowFastSplitting()`
- `fast_model` (Type): The model describing fast dynamics.
    Default: `AtmosAcousticGravityLinearModel`
- `mis_method` (Function): Function defining the particular MIS
    method to be used.
    Default: `MIS2`
- `fast_method` (Function): Function defining the fast solver.
    Default: `LSRK54CarpenterKennedy`
- `nsubsteps` (Int): Integer denoting the total number of times
    to substep the fast process.
    Default: `50`
- `discrete_splitting` (Boolean): Boolean denoting whether a PDE level or
    discretized level splitting should be used. If `true` then the PDE is
    discretized in such a way that `f_fast + f_slow` is equivalent to
    discretizing the original PDE directly.
    Default: `false`

### References
    @article{KnothWensch2014,
        title={Generalized split-explicit Runge--Kutta methods for
            the compressible Euler equations},
        author={Knoth, Oswald and Wensch, Joerg},
        journal={Monthly Weather Review},
        volume={142},
        number={5},
        pages={2067--2081},
        year={2014}
    }
"""
struct MISSolverType{DS} <: AbstractSolverType
    # The type of discrete splitting to apply to the right-hand side
    splitting_type::DS
    # The model describing fast dynamics
    fast_model::Type
    # Main MIS function
    mis_method::Function
    # Fast RK solver
    fast_method::Function
    # Substepping parameter for the fast processes
    nsubsteps::Int
    # Whether to use a PDE level or discrete splitting
    discrete_splitting::Bool

    function MISSolverType(;
        splitting_type = SlowFastSplitting(),
        fast_model = AtmosAcousticGravityLinearModel,
        mis_method = MIS2,
        fast_method = LSRK54CarpenterKennedy,
        nsubsteps = 50,
        discrete_splitting = false,
    )

        DS = typeof(splitting_type)

        return new{DS}(
            splitting_type,
            fast_model,
            mis_method,
            fast_method,
            nsubsteps,
            discrete_splitting,
        )
    end
end

"""
    getdtmodel(ode_solver::MISSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(ode_solver::MISSolverType, bl)
    # Most restrictive dynamics are part of the fast model
    return ode_solver.fast_model(bl)
end

"""
# Description
    function solversetup(
        ode_solver::MISSolverType{SlowFastSplitting},
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an ODE solver for the partition slow-fast ODE
using an MIS method with explicit time-integration.
The splitting of the fast (acoustic and gravity waves)
dynamics is done in _all_ spatial directions.
"""
function solversetup(
    ode_solver::MISSolverType{SlowFastSplitting},
    dg,
    Q,
    dt,
    t0,
    diffusion_direction,
)

    # Extract fast model and define a DG model
    # for the fast processes (acoustic/gravity waves
    # in all spatial directions)
    fast_model = ode_solver.fast_model(dg.balance_law)
    fast_dg = DGModel(
        fast_model,
        dg.grid,
        dg.numerical_flux_first_order,
        dg.numerical_flux_second_order,
        dg.numerical_flux_gradient,
        state_auxiliary = dg.state_auxiliary,
        state_gradient_flux = dg.state_gradient_flux,
        states_higher_order = dg.states_higher_order,
        direction = EveryDirection(),
    )

    # Using the RemainderModel, we subtract away the
    # fast processes and define a DG model for the
    # slower processes (advection and diffusion)
    if ode_solver.discrete_splitting
        numerical_flux_first_order =
            (dg.numerical_flux_first_order, (dg.numerical_flux_first_order,))
    else
        numerical_flux_first_order = dg.numerical_flux_first_order
    end
    slow_dg = remainder_DGModel(
        dg,
        (fast_dg,);
        numerical_flux_first_order = numerical_flux_first_order,
    )

    solver = ode_solver.mis_method(
        slow_dg,
        fast_dg,
        ode_solver.fast_method,
        ode_solver.nsubsteps,
        Q;
        dt = dt,
        t0 = t0,
    )

    return solver
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module HydroPowerSimulations

#################################################################################
# Exports
export HydroEnergyCascade
export HydroDispatchCascade
export HydroDispatchReservoirCascade
export HydroDispatchRunOfRiverCascade
export HydroDispatchReservoirBudgetUpperBound

#################################################################################
# Imports
using PowerSystems
import InfrastructureSystems
import Dates
import TimeZones
import PowerSimulations
import PowerModels
import JuMP
import CSV
import DataFrames
import ParameterJuMP

const PSY = PowerSystems
const IS = InfrastructureSystems
const PM = PowerModels
const PSI = PowerSimulations
const PJ = ParameterJuMP

#################################################################################
# Includes
include("models/HydroCascade.jl")
include("models/generated/includes.jl")
include("devices_models/devices/HydroReservoirCascade.jl")
include("devices_models/devices/HydroDispatchCascade.jl")
include("devices_models/devices/HydroReservoirBudgetUpperBound.jl")
include("devices_models/device_constructors/HydroReservoirCascade_constructor.jl")
include("devices_models/device_constructors/HydroDispatchCascade_constructor.jl")
include("devices_models/device_constructors/HydroReservoirBudgetUpperBound_constructor.jl")
include("parsers/set_upstream.jl")

end # module


text:
"""
    abstract type BasisSet

An abstract type that represents a set of orthonormal states that can be used to generate
a vector space.

The basis sets are assumed to be indexed with indices in `1:length(b)`.

# Interface

Child types must implement `Base.length`, which should return the number of basis elements
(or, equivalently, the dimension of the generated vector space).
"""
abstract type BasisSet end

"""
    abstract type LinearOperator{B <: BasisSet}

# Interface

Where `L <: LinearOperator` is the user-defined linear operator type:

* `apply(i::Integer, op::L, j::Integer)`: return the `i`-th compoment of the vector that results
  from applying operator `op` on the basis state `j`.
* `basis(op::L)`: return the basis set
"""
abstract type LinearOperator{B <: BasisSet} end

"""
    abstract type BasisTransformation{B1 <: BasisSet, B2 <: BasisSet}

Abstract type for representations of basis transformations between two basis sets, from `B1` to `B2`.

# Interface

* `basisfrom(::BasisTransformation)`: return the basis `B1`
* `basisto(::BasisTransformation)`: return the basis `B2`
* `transform(i::Integer, ::BasisTransformation, j::Integer)`: return the `i`-th component of the vector
  in `B2` corresponding to the transformed `j`-th vector of `B1`.
"""
abstract type BasisTransformation{B1 <: BasisSet, B2 <: BasisSet} end

function apply end
function basis end
function basisfrom end
function basisto end
function transform end

struct VSVector{B <: BasisSet}
    cs :: Vector{ComplexF64}
end
Base.length(v::VSVector) = length(v.cs)
Base.:(+)(v1::VSVector{B}, v2::VSVector{B}) where {B <: BasisSet} = VSVector{B}(v1.cs .+ v2.cs)

function apply(op::LinearOperator{B}, x::VSVector{B}) where B <: BasisSet
    N = length(basis(op))
    VSVector{B}(ComplexF64[
        sum(apply(i, op, j) * x.cs[j] for j = 1:N)
        for i = 1:N
    ])
end

function apply(op::LinearOperator{B}, cs::Vector) where B <: BasisSet
    @assert length(cs) == length(basis(op))
    apply(op, VSVector{B}(cs))
end

"""
    matrix(op::LinearOperator)

Returns a matrix representation of the operator `op`.
"""
function matrix(op::LinearOperator)
    N = length(basis(op))
    [apply(i, op, j) for i = 1:N, j = 1:N]
end

"""
    matrix(bt::BasisTransformation)

Return the matrix representation of the basis transformation `bt`.
"""
function matrix(bt::BasisTransformation)
    N1, N2 = length(basisfrom(bt)), length(basisto(bt))
    [transform(i, bt, j) for i = 1:N2, j = 1:N1]
end

struct IdentityOperator{B <: BasisSet} <: LinearOperator{B}
    b :: B
end
basis(op::IdentityOperator) = op.b
apply(i::Integer, op::IdentityOperator, j::Integer) = ComplexF64(i == j ? 1.0 : 0.0)


text:
export LibGitError

module GitErrorConst
    import ..cint
    import ..cuint

    const GIT_OK          = cint(0)
    const ERROR           = cint(-01)
    const ENOTFOUND       = cint(-03)
    const EEXISTS         = cint(-04)
    const EAMBIGUOUS      = cint(-05)
    const EBUFS           = cint(-06)
    const EUSER           = cint(-07)
    const EBAREREPO       = cint(-08)
    const EUNBORNBRANCH   = cint(-09)
    const EUNMERGED       = cint(-10)
    const ENONFASTFORWARD = cint(-11)
    const EINVALIDSPEC    = cint(-12)
    const EMERGECONFLICT  = cint(-13)
    const ELOCKED         = cint(-14)
    const PASSTHROUGH     = cint(-30)
    const ITEROVER        = cint(-31)
end

const git_error_code = (Int => Symbol)[
     00 => :OK,             # no error
    -01 => :Error,          # generic error
    -03 => :NotFound,       # requested object could not be found
    -04 => :Exists,         # object exits preventing op
    -05 => :Ambiguous,      # more than one object matches
    -06 => :Bufs,           # output buffer too small to hold data
    -07 => :User,           # user callback generated error
    -08 => :BareRepo,       # operation not allowed on bare repo
    -09 => :UnbornBranch,   # HEAD refers to branch with 0 commits
    -10 => :Unmerged,       # merge in progress prevented op
    -11 => :NonFastForward, # ref not fast-forwardable
    -12 => :InvalidSpec,    # name / ref not in valid format
    -13 => :MergeConflict,  # merge conflict prevented op
    -14 => :Locked,         # lock file prevented op
    -15 => :Modified,       # ref value does not match expected
    -31 => :Iterover        # signals end of iteration
]

const git_error_class = (Int => Symbol)[
     0 => :None,
     1 => :NoMemory,
     2 => :OS,
     3 => :Invalid,
     4 => :Ref,
     5 => :Zlib,
     6 => :Repo,
     7 => :Config,
     8 => :Regex,
     9 => :Odb,
    10 => :Index,
    11 => :Object,
    12 => :Net,
    13 => :Tag,
    14 => :Tree,
    15 => :Indexer,
    16 => :SSL,
    17 => :Submodule,
    18 => :Thread,
    19 => :Stash,
    20 => :Checkout,
    21 => :FetchHead,
    22 => :Merge,
    23 => :SSH,
    24 => :Filter,
    25 => :Revert,
    26 => :Callback,
    27 => :CherryPick
]

immutable ErrorStruct
    message::Ptr{Uint8}
    class::Cint
end

immutable LibGitError{Class, Code}
    msg::UTF8String
end

function last_error()
    err = ccall((:giterr_last, libgit2), Ptr{ErrorStruct}, ())
    err_obj   = unsafe_load(err)
    err_class = git_error_class[int(err_obj.class)]
    err_msg   = bytestring(err_obj.message)
    return (err_class, err_msg)
end

LibGitError(code::Integer) = begin
    err_code = git_error_code[int(code)]
    err_class, err_msg = last_error()
    return LibGitError{err_class, err_code}(err_msg)
end

macro check(git_func)
    quote
        local err::Cint
        err = $(esc(git_func::Expr))
        if err < 0
            throw(LibGitError(err))
        end
        err
    end
end


text:
module FiniteElementInterface

abstract type AbstFE{T} end

# TODO
# * integrate

end

max_stars_repo_name: OptimalDesignLab/CSEInterfaces

max_issues_repo_name: OptimalDesignLab/CSEInterfaces

max_forks_repo_name: OptimalDesignLab/CSEInterfaces


text:
import Pkg;
# use a separate Project for building docs; see
# https://discourse.julialang.org/t/psa-use-a-project-for-building-your-docs/14974
Pkg.activate("docs")

# to find parent package when executing this file using
# `include("docs/make.jl")`
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end

# to find parent package from examples, doctests, etc. (executed in
# working directory docs/build/)
if !("../.." in LOAD_PATH)
    push!(LOAD_PATH, "../..")
end

using Documenter, ABCDBeamTrace

makedocs(
    sitename = "ABCDBeamTrace",
    format = Documenter.HTML(
        prettyurls = false # put everything in one file
    ),
    modules = [ABCDBeamTrace]
)

# re-activate parent Project
Pkg.activate(".")


text:
# Device type and auxiliary functions

export
    CuDevice, name, totalmem, attribute

"""
    CuDevice(i::Integer)

Get a handle to a compute device.
"""
struct CuDevice
    handle::CUdevice

    # CuDevice is just an integer, but we need (?) to call cuDeviceGet to make sure this
    # integer is valid. to avoid ambiguity, add a bogus argument (cfr. `checkbounds`)
    CuDevice(::Type{Bool}, handle::CUdevice) = new(handle)
end

const DEVICE_CPU = CuDevice(Bool, CUdevice(-1))
const DEVICE_INVALID = CuDevice(Bool, CUdevice(-2))

function CuDevice(ordinal::Integer)
    device_ref = Ref{CUdevice}()
    cuDeviceGet(device_ref, ordinal)
    CuDevice(Bool, device_ref[])
end

Base.convert(::Type{CUdevice}, dev::CuDevice) = dev.handle

Base.:(==)(a::CuDevice, b::CuDevice) = a.handle == b.handle
Base.hash(dev::CuDevice, h::UInt) = hash(dev.handle, h)

function Base.show(io::IO, dev::CuDevice)
  print(io, "CuDevice($(dev.handle)): ")
  if dev == DEVICE_CPU
      print(io, "CPU")
  elseif dev == DEVICE_INVALID
      print(io, "INVALID")
  else
      print(io, "$(name(dev))")
  end
end

"""
    name(dev::CuDevice)

Returns an identifier string for the device.
"""
function name(dev::CuDevice)
    buflen = 256
    buf = Vector{Cchar}(undef, buflen)
    cuDeviceGetName(pointer(buf), buflen, dev)
    buf[end] = 0
    return unsafe_string(pointer(buf))
end

"""
    totalmem(dev::CuDevice)

Returns the total amount of memory (in bytes) on the device.
"""
function totalmem(dev::CuDevice)
    mem_ref = Ref{Csize_t}()
    cuDeviceTotalMem(mem_ref, dev)
    return mem_ref[]
end

"""
    attribute(dev::CuDevice, code)

Returns information about the device.
"""
function attribute(dev::CuDevice, code::CUdevice_attribute)
    value_ref = Ref{Cint}()
    cuDeviceGetAttribute(value_ref, code, dev)
    return value_ref[]
end

@enum_without_prefix CUdevice_attribute CU_

## device iteration

export devices

struct DeviceSet end

"""
    devices()

Get an iterator for the compute devices.
"""
devices() = DeviceSet()

Base.eltype(::DeviceSet) = CuDevice

function Base.iterate(iter::DeviceSet, i=1)
    i >= length(iter) + 1 ? nothing : (CuDevice(i-1), i+1)
end

function Base.length(::DeviceSet)
    count_ref = Ref{Cint}()
    cuDeviceGetCount(count_ref)
    return count_ref[]
end

Base.IteratorSize(::DeviceSet) = Base.HasLength()

## convenience attribute getters

export warpsize, capability

"""
    warpsize(dev::CuDevice)

Returns the warp size (in threads) of the device.
"""
warpsize(dev::CuDevice) = attribute(dev, DEVICE_ATTRIBUTE_WARP_SIZE)

"""
    capability(dev::CuDevice)

Returns the compute capability of the device.
"""
function capability(dev::CuDevice)
    return VersionNumber(attribute(dev, DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MAJOR),
                         attribute(dev, DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MINOR))
end


text:
# ==============================================================================
# TODO Maxplus eigenvalues eigenvectors (Howard algorithm)


text:
################################################################################
#
#   NfMaxOrd/FracIdl.jl : Fractional ideals of maximal orders in
#                         absolute number fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2015, 2016 Tommy Hofmann
# (C) 2016, 2016 Claus Fieker
#
################################################################################

export basis_mat, norm, inv, ==, *

################################################################################
#
#  Deepcopy
#
################################################################################

function deepcopy(x::NfMaxOrdFracIdl)
  z = NfMaxOrdFracIdl(deepcopy(x.num), deepcopy(x.den))
  if isdefined(x, :basis_mat)
    z.basis_mat = deepcopy(x.basis_mat)
  end
  return z
end

################################################################################
#
#  Basis matrix
#
################################################################################

doc"""
***
    basis_mat(I::NfOrdFracIdl) -> FakeFmpqMat

> Returns the basis matrix of $I$ with respect to the basis of the order.
"""
function basis_mat(x::NfMaxOrdFracIdl)
  return FakeFmpqMat(basis_mat(num(x)), den(x))
end

################################################################################
#
#  Numerator and denominator
#
################################################################################

num(x::NfMaxOrdFracIdl) = x.num

den(x::NfMaxOrdFracIdl) = x.den

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, s::NfMaxOrdFracIdlSet)
   print(io, "Set of fractional ideals of ")
   print(io, s.order)
end

function show(io::IO, id::NfMaxOrdFracIdl)
  print(io, "1//", id.den, " * ")
  show(io, id.num)
end

################################################################################
#
#  Norm
#
################################################################################

doc"""
***
    norm(I::NfOrdFracIdl) -> fmpq

> Returns the norm of $I$
"""
function norm(A::NfMaxOrdFracIdl)
  return norm(A.num)//A.den^degree(A.num.parent.order)
end

################################################################################
#
#  Minimum
#
################################################################################

function minimum(A::NfMaxOrdFracIdl)
  return minimum(A.num)//A.den
end

################################################################################
#
#  Inverse
#
################################################################################

doc"""
***
    inv(A::NfMaxOrdFracIdl) -> NfMaxOrdFracIdl

> Returns the fractional ideal $B$ such that $AB = \mathcal O$.
"""
function inv(A::NfMaxOrdFracIdl)
  B = inv(A.num)
  g = gcd(B.den, A.den)
  B.den = divexact(B.den, g)
  a = divexact(A.den, g)
  return prod_by_int(B, a)
end

################################################################################
#
#  Simplification
#
################################################################################

function simplify_exact(A::NfMaxOrdFracIdl)
  g = A.den

  A.den = divexact(A.den, g)

  @assert A.den == 1

  if has_2_elem(A.num)
    A.num.gen_one = divexact(A.num.gen_one, g)
    A.num.gen_two = divexact(A.num.gen_two, g)
  end

  if has_minimum(A.num)
    A.num.minimum = divexact(A.num.minimum, g)
  end
  if has_norm(A.num)
    A.num.norm = divexact(A.num.norm, g^degree(A.num.parent.order))
  end
  if has_basis_mat(A.num)
    A.num.basis_mat = divexact(A.num.basis_mat, g)
  end
end

function simplify(A::NfMaxOrdFracIdl)
  simplify(A.num)
  b = basis_mat(A.num)
  g = den(A)
  for i in 1:rows(b)
    for j in 1:cols(b)
      g = gcd(g, b[i, j])
    end
  end

  if g != 1
    if has_2_elem(A.num)
      A.num.gen_one = divexact(A.num.gen_one, g)
      A.num.gen_two = divexact(A.num.gen_two, g)
    end
    A.den = divexact(A.den, g)
    if has_minimum(A.num)
      A.num.minimum = divexact(A.num.minimum, g)
    end
    if has_norm(A.num)
      A.num.norm = divexact(A.num.norm, g^degree(A.num.parent.order))
    end
    if has_basis_mat(A.num)
      A.num.basis_mat = divexact(A.num.basis_mat, g)
    end
  end

#  m = minimum(A.num)
#  println("minimum is $m")
#  g = gcd(m, A.den)
#  println("gcd is $g")
#  println("now do something with $(A.num)")
#  if g != 1
#    if has_2_elem(A.num)
#      A.num.gen_one = divexact(A.num.gen_one, g)
#      A.num.gen_two = divexact(A.num.gen_two, g)
#    end
#    A.den = divexact(A.den, g)
#    if has_minimum(A.num)
#      A.num.minimum = divexact(A.num.minimum, g)
#    end
#    if has_norm(A.num)
#      A.num.norm = divexact(A.num.norm, g^degree(A.num.parent.order))
#    end
#    if has_basis_mat(A.num)
#      A.num.basis_mat = divexact(A.num.basis_mat, g)
#    end
#    simplify(A.num)
#  end
  return A
end

################################################################################
#
#  Primeness
#
################################################################################

is_prime_known(A::NfMaxOrdFracIdl) = is_prime_known(A.num)

################################################################################
#
#  Equality
#
################################################################################

doc"""
***
    ==(x::NfOrdFracIdl, y::NfOrdFracIdl) -> Bool

> Returns whether $x$ and $y$ are equal.
"""
function ==(A::NfMaxOrdFracIdl, B::NfMaxOrdFracIdl)
  #return B.den * basis_mat(A.num) == A.den * basis_mat(B.num)
  D = inv(B)
  E = prod(A, D)
  C = simplify(E)
  return norm(C) == 1 && C.den == 1
end

################################################################################
#
#  Binary operations
#
################################################################################

function prod(a::NfMaxOrdFracIdl, b::NfMaxOrdFracIdl)
  A = a.num* b.num
  return NfMaxOrdFracIdl(A, a.den*b.den)
end

doc"""
***
    *(I::NfMaxOrdFracIdl, J::NfMaxOrdFracIdl) -> NfMaxOrdFracIdl

> Returns $IJ$.
"""
*(A::NfMaxOrdFracIdl, B::NfMaxOrdFracIdl) = prod(A, B)

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function prod_by_int(A::NfMaxOrdFracIdl, a::fmpz)
  return NfMaxOrdFracIdl(A.num * a, A.den)
end

function *(A::NfMaxOrdIdl, B::NfMaxOrdFracIdl)
  z = NfMaxOrdFracIdl(A*B.num, B.den)
  return z
end

*(A::NfMaxOrdFracIdl, B::NfMaxOrdIdl) = NfMaxOrdFracIdl(A.num*B, A.den)

function *(A::NfMaxOrdFracIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfMaxOrdFracIdl) = A*a

function *(A::NfMaxOrdIdl, a::nf_elem)
  C = *(A, a*order(A))
  return C
end

*(a::nf_elem, A::NfMaxOrdIdl) = A*a

function /(A::NfMaxOrdFracIdl, B::NfMaxOrdIdl)
  C = prod(A, inv(B))
  return C
end

function /(A::NfMaxOrdFracIdl, a::nf_elem)
  C = prod(A, Idl((order(A), inv(a))))
  return C
end

function *(x::nf_elem, y::NfMaxOrd)
  b, z = _check_elem_in_order(den(x, y)*x, y)
  return NfMaxOrdFracIdl(ideal(y, y(z)), den(x, y))
end

################################################################################
#
#  Conversion
#
################################################################################

function call(ord::NfMaxOrdIdlSet, b::NfMaxOrdFracIdl)
   b.den > 1 && error("not integral")
   return b.num
end


text:
using ArgParse
using Random
using JLD
using NPZ
include("./DEAC.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--temperature", "-T"
            help = "Temperature of system."
            arg_type = Float64
            default = 0.0
        "--number_of_generations", "-N"
            help = "Number of generations before genetic algorithm quits."
            arg_type = Int64
            default = 100000
        "--population_size", "-P"
            help = "Size of initial population"
            arg_type = Int64
            default = 512
        "--genome_size", "-M"
            help = "Size of genome."
            arg_type = Int64
            default = 600
        "--omega_max"
            help = "Maximum frequency to explore."
            arg_type = Float64
            default = 60.0
        "--smooth"
            help = "Use smoothing factor in fitness."
            action = :store_true
        "--normalize"
            help = "Normalize spectrum to the zeroeth moment."
            action = :store_true
        "--use_inverse_first_moment"
            help = "Calculate inverse first moment from ISF data and use it in fitness."
            action = :store_true
        "--first_moment"
            help = "FIXME First moment."
            arg_type = Float64
            default = -1.0
        "--third_moment"
            help = "FIXME Third moment."
            arg_type = Float64
            default = -1.0
        "--third_moment_error"
            help = "FIXME Third moment."
            arg_type = Float64
            default = 0.01
        "--crossover_probability", "-r"
            help = "Initial probability for parent gene to become mutant vector gene."
            arg_type = Float64
            default = 0.90
        "--self_adapting_crossover_probability", "-u"
            help = "Probability for `crossover_probability` to mutate."
            arg_type = Float64
            default = 0.10
        "--differential_weight", "-F"
            help = "Initial weight factor when creating mutant vector."
            arg_type = Float64
            default = 0.90
        "--self_adapting_differential_weight_probability", "-v"
            help = "Probability for `differential_weight` to mutate."
            arg_type = Float64
            default = 0.10
        "--self_adapting_differential_weight_shift", "-l"
            help = "If `self_adapting_differential_weight_probability` mutate, new value is `l + m*rand()`."
            arg_type = Float64
            default = 0.10
        "--self_adapting_differential_weight", "-m"
            help = "If `self_adapting_differential_weight_probability` mutate, new value is `l + m*rand()`."
            arg_type = Float64
            default = 0.90
        "--reject"
            help = "Name of rejection function."
            arg_type = String
            default = "smallerFit"
        "--frequency_bins_method"
            help = "Name of function to generate frequency bins."
            arg_type = String
            default = "fixed_frequency_bins"
        "--load_frequency_bins_filename"
            help = "A JLD file containing frequency_bins to load."
            arg_type = String
            default = "./frequency.jld"
        "--stop_minimum_fitness"
            help = "Stopping criteria, if minimum fitness is below `stop_minimum_fitness` stop evolving."
            arg_type = Float64
            default = 1.0
        "--seed"
            help = "Seed to pass to random number generator."
            arg_type = Int64
            default = 1
        "--save_state"
            help = "Save state of DEAC algorithm. Saves the random number generator, population, and population fitness."
            action = :store_true
        "--save_file_dir"
            help = "Directory to save results in."
            arg_type = String
            default = "./deacresults"
        "--track_stats"
            help = "Track minimum fitness and other stats."
            action = :store_true
        "--number_of_blas_threads"
            help = "Number of BLAS threads to utilize if > 0 otherwise uses Julia default value."
            arg_type = Int64
            default = 0
        "qmc_data"
            help = "File containing quantum Monte Carlo data with columns: IMAGINARY_TIME, INTERMEDIATE_SCATTERING_FUNCTION, ERROR"
            arg_type = String
            required = true
    end

    return parse_args(s)
end

function checkMoves(argname::String,s::Array{String,1},l::Array{String,1})
    for ss in s
        checkMoves(argname,ss,l)
    end
end

function checkMoves(argname::String,s::String,l::Array{String,1})
    if !(s in l)
        print("$s is not a valid parameter for $argname. Valid parameters are: ")
        println(l)
        error("Failed to validate argument: $argname")
    end
nothing
end

function main()
    start = time();
    parsed_args = parse_commandline()

    # create data directory
    try
        mkpath(parsed_args["save_file_dir"]);
    catch
        nothing
    end
    save_dir = parsed_args["save_file_dir"];

    checkMoves("reject",parsed_args["reject"],DEAC.DEACreject.functionNames)
    checkMoves("frequency_bins_method",parsed_args["frequency_bins_method"],DEAC.DEACfrequency.functionNames)
    get_frequency_bins = getfield(DEAC.DEACfrequency, Symbol(parsed_args["frequency_bins_method"]));

    if parsed_args["frequency_bins_method"] == "load_frequency_bins"
        frequency_bins_filename = parsed_args["load_frequency_bins_filename"];
        println("loading frequency_bins from $(frequency_bins_filename) (will override genome_size and omega_max");
        frequency_bins = get_frequency_bins(frequency_bins_file);
    else
        frequency_bins = get_frequency_bins(0.0,parsed_args["omega_max"],parsed_args["genome_size"]);
    end
    #FIXME ensure frequency odd
    if size(frequency_bins,1)%2 == 0
        throw(AssertionError("frequency space must be odd FIXME"))
    end

    _ext = splitext(parsed_args["qmc_data"])[2];
    if _ext == ".npz"
        qmcdata=NPZ.npzread(parsed_args["qmc_data"]);
    elseif _ext == ".jld"
        qmcdata = load(parsed_args["qmc_data"]);
    else
        throw(AssertionError("qmc_data must be *.jld or *.npz"));
    end
    tau = qmcdata["tau"];
    isf = qmcdata["isf"];
    err = qmcdata["error"];

    u4,omega,minS,minP,minFit,total_generations,
    avgFitness,minFitness,
    stdFitness,P,
    fitnessP,rng,
    crossover_probs,differential_weights = DEAC.deac( tau,isf,err,frequency_bins,
                           use_inverse_first_moment = parsed_args["use_inverse_first_moment"],
                           smooth = parsed_args["smooth"],
                           normalize = parsed_args["normalize"],
                           first_moment = parsed_args["first_moment"],
                           third_moment = parsed_args["third_moment"],
                           third_moment_error = parsed_args["third_moment_error"],
                           temperature = parsed_args["temperature"],
                           number_of_generations = parsed_args["number_of_generations"],
                           population_size = parsed_args["population_size"],
                           crossoverProb = parsed_args["crossover_probability"],
                           SAcrossoverProb = parsed_args["self_adapting_crossover_probability"],
                           differentialWeight = parsed_args["differential_weight"],
                           SAdifferentialWeightProb = parsed_args["self_adapting_differential_weight_probability"],
                           SAdifferentialWeightShift = parsed_args["self_adapting_differential_weight_shift"],
                           SAdifferentialWeight = parsed_args["self_adapting_differential_weight"],
                           rejectFunc = parsed_args["reject"],
                           stop_minimum_fitness = parsed_args["stop_minimum_fitness"],
                           seed = parsed_args["seed"],
                           track_stats = parsed_args["track_stats"],
                           number_of_blas_threads = parsed_args["number_of_blas_threads"])
    seed = lpad(parsed_args["seed"],4,'0');
    elapsed = time() - start;
    filename = "$(save_dir)/deac_results_$(seed)_$u4.jld";
    println("Saving results to $filename");
    save(filename,
         "u4",u4,
         "frequency",omega,
         "dsf",minS,
         "minP",minP,
         "minFit",minFit,
         "total_generations",total_generations,
         "elapsed_time",elapsed);

    if parsed_args["track_stats"]
        filename = "$(save_dir)/deac_stats_$(seed)_$u4.jld";
        println("Saving stats to $filename");
        save(filename,
             "u4",u4,
             "avgFitness",avgFitness,
             "minFitness",minFitness,
             "stdFitness",stdFitness)
    end

    filename = "$(save_dir)/deac_params_$(seed)_$u4.jld";
    println("Saving parameters to $filename");
    save(filename,
         "u4",u4,
         "tau",tau,
         "isf",isf,
         "isf_error",err,
         "frequency",frequency_bins,
         "first_moment",parsed_args["first_moment"],
         "third_moment",parsed_args["third_moment"],
         "third_moment_error",parsed_args["third_moment_error"],
         "temperature",parsed_args["temperature"],
         "number_of_generations",parsed_args["number_of_generations"],
         "population_size",parsed_args["population_size"],
         "genome_size",parsed_args["genome_size"],
         "omega_max",parsed_args["omega_max"],
         "crossoverProb",parsed_args["crossover_probability"],
         "SAcrossoverProb",parsed_args["self_adapting_crossover_probability"],
         "differentialWeight",parsed_args["differential_weight"],
         "SAdifferentialWeightProb",parsed_args["self_adapting_differential_weight_probability"],
         "SAdifferentialWeightShift",parsed_args["self_adapting_differential_weight_shift"],
         "SAdifferentialWeight",parsed_args["self_adapting_differential_weight"],
         "rejectFunc",parsed_args["reject"],
         "frequency_bins_method",parsed_args["frequency_bins_method"],
         "load_frequency_bins_filename",parsed_args["load_frequency_bins_filename"],
         "stop_minimum_fitness",parsed_args["stop_minimum_fitness"],
         "seed",parsed_args["seed"])
    if parsed_args["save_state"]
        filename = "$(save_dir)/deac_state_$(seed)_$(seed_resample)_$u4.jld";
        println("Saving state to $filename");
        save(filename,
             "u4",u4,
             "P",P,
             "fitness",fitnessP,
             "rng",rng,
             "crossover_probs",crossover_probs,
             "differential_weights",differential_weights)
    end
    nothing
end

main()


text:
using Printf

function print_header(ws::COSMO.Workspace)

	m, n = ws.p.model_size
	settings = ws.settings
	settings.scaling > 0 ? scaling_status = "on" : scaling_status = "off"
	nnz_in_P = count(!iszero,ws.p.P) - count(!iszero,diag(ws.p.P)) + n
	nnz_in_M = 2*count(!iszero,ws.p.A) + nnz_in_P + m
	println("-"^66 * "\n" * " "^10 * "COSMO v0.4.4 - A Quadratic Objective Conic Solver\n" * " "^25 * "Michael Garstka\n"  * " "^16 * "University of Oxford, 2017 - 2019\n" * "-"^66 * "\n")

	println("Problem:  x ∈ R^{$(n)},\n          constraints: A ∈ R^{$(m)x$(n)} ($(count(!iszero, ws.p.A)) nnz),\n          matrix size to factor: $(n + m)x$(n + m) ($(nnz_in_M) nnz)")
	for (iii, set) in enumerate(sort(ws.p.C.sets, by = x -> -x.dim))
		set_name = split(string(typeof(set).name), ".")[end]
		iii == 1 ? println("Sets:"*" "^5*"$(set_name) of dim: $(set.dim)") : println(" "^10*"$(set_name) of dim: $(set.dim)")
		if iii > 5
			println(" "^10*"... and $(length(ws.p.C.sets)-5) more")
			break
		end
	end
	if ws.settings.decompose && ws.ci.num_decomposable > 0
		println("Decomp:   Num of original PSD cones: $(ws.ci.num_psd_cones)\n" * " "^10*"Num decomposable PSD cones: $(ws.ci.num_decomposable)\n" * " "^10*"Num PSD cones after decomposition: $(ws.ci.num_decom_psd_cones)")
	end
	println("Settings: ϵ_abs = $(@sprintf("%.1e", settings.eps_abs)), ϵ_rel = $(@sprintf("%.1e", settings.eps_rel)),\n" * " "^10 * "ϵ_prim_inf = $(@sprintf("%.1e", settings.eps_prim_inf)), ϵ_dual_inf = $(@sprintf("%.1e", settings.eps_dual_inf)),\n" * " "^10 * "ρ = $(settings.rho), σ = $(settings.sigma), α = $(settings.alpha),\n" * " "^10 * "max_iter = $(settings.max_iter),\n" * " "^10 * "scaling iter = $(settings.scaling) ($(scaling_status)),\n" * " "^10 * "check termination every $(settings.check_termination) iter,\n" * " "^10 * "check infeasibility every $(settings.check_infeasibility) iter,\n" *	 " "^10 * "KKT system solver: $(print_lin_sys(settings.kkt_solver.ObjectType))")

	println("Setup Time: $(round.(ws.times.setup_time*1000; digits=2))ms\n")
	println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:")
	nothing
end

function print_iteration(ws::COSMO.Workspace, iter::Int64, cost::Float64, r_prim::Float64, r_dual::Float64)
	settings = ws.settings
	if mod(iter, 1) == 0 || iter == 1 || iter == 2 || iter == settings.max_iter
		if mod(iter, settings.check_termination) == 0
			@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\n", iter, cost, r_prim, r_dual, ws.ρ)
		else
			@printf("%d\t%.4e\t ---\t\t\t---\n", iter, cost)
		end
	end
	nothing
end

function print_result(status::Symbol, iter::Int64, cost::Float64, rt::Float64)
	println("\n" * "-"^66 * "\n>>> Results\nStatus: $(status)\nIterations: $(iter)\nOptimal objective: $(round.(cost; digits = 4))\nRuntime: $(round.(rt; digits = 3))s ($(round.(rt * 1000; digits = 2))ms)\n")
	nothing
end

function print_lin_sys(s::Type{<:AbstractKKTSolver})
	lin_sys_dict = Dict("QdldlKKTSolver" => "QDLDL", "CholmodKKTSolver" => "CHOLMOD", "PardisoDirectKKTSolver" => "Pardiso (direct)", "PardisoIndirectKKTSolver" => "Pardiso (indirect)",  "MKLPardisoKKTSolver" => "MKL Pardiso")
	return get(lin_sys_dict, string(s), string(s))
end


text:
using SubmodularMaximization
using POMDPs
using MCTS
using D3Trees

grid_size = 10
iterations = 10
horizon = SubmodularMaximization.default_horizon

grid = Grid(grid_size, grid_size)
sensor = RangingSensor()

@show robot_state = random_state(grid)

histogram_filter = Filter(grid)

problem = SingleRobotTargetTrackingProblem(grid, sensor, horizon,
                                           [histogram_filter])

solver = generate_solver(horizon, n_iterations = iterations,
                         enable_tree_vis=true)
policy = solve(solver, problem)

a, info = action_info(policy, MDPState(robot_state))

D3Tree(info[:tree])

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using Test, Pkg

@testset "Common" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false
    for submodule in ["Thermodynamics", "SurfaceFluxes"]
        if all_tests || "$submodule" in ARGS || "Common" in ARGS
            include_test(submodule)
        end
    end

end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
using Gadfly

n = 10
plot(x=fill(1, n), y=rand(n), Geom.bar)


text:
using Documenter;
push!(LOAD_PATH, "../src/");
using SharedArrays, SpecialFunctions, Distributed, Printf
using DelimitedFiles, LinearAlgebra
include("../src/AARMBEMC_ReadCMDARGS.jl");
include("../src/E1func.jl");
include("../src/GenericMatrixVectorUtils.jl");
include("../src/MolData.jl");
include("../src/TransData.jl");
include("../src/PeriodicData.jl");
include("../src/MolTransUtils.jl");
include("../src/MolAlpha.jl");
include("../src/DeltaKIdata.jl");
include("../src/GFVACGG.jl");
include("../src/GFPECGG.jl");
include("../src/MolAssembleGFGG.jl");
include("../src/heattransfer.jl");
include("../src/vdWenergy.jl");

makedocs(
    sitename = "AARMBEM-C.jl",
    authors = string("Prashanth S. Venkataram, ",
                     "Jan Hermann, ",
                     "Alexandre Tkatchenko, ",
                     "and Alejandro W. Rodriguez"),
    pages = [
        "Home" => "index.md",
        "Background" => [
            "History" => "theory/historicalbackground.md",
            "Theory" => "theory/theoreticalbackground.md",
            "Common acronyms" => "theory/acronyms.md"
        ],
        "Implementation Notes" => [
            "Units" => "implementation/units.md",
            "Morse potential for intermolecular internuclear couplings" => "implementation/morse.md",
            "Matrix storage conventions" => "implementation/matrixstorage.md",
            "Data structures and rigid transformations" => "implementation/datastructuresandtransformations.md",
            "Performance notes and required packages" => "implementation/performancenotes.md",
            "Output quantities" => "implementation/outputquantities.md",
            "Warning: Ewald summation" => "implementation/ewaldwarning.md",
            "Warning: vdW interactions" => "implementation/vdWwarning.md"
        ],
        "Usage" => [
            "Configuration files" => "usage/configfiles.md",
            "Transformation files" => "usage/transfiles.md",
            "Bloch periodicity files" => "usage/periodicfiles.md",
            "Frequency and Bloch wavevector files" => "usage/freqkfiles.md",
            "Command-line arguments" => "usage/CMDARGS.md",
            "Output file structure" => "usage/outfiles.md"
        ],
        "Examples" => [
            "Two identical mirrored collinear carbyne wires (vacuum)" => "examples/2C250carbynecollinearVAC.md"
        ],
        "API" => [
            "Command-line arguments" => "api/CMDARGS.md",
            "Exponential integral function" => "api/expintegral.md",
            "Generic matrix and vector utilities" => "api/genericmatrixvectorutils.md",
            "Maxwell Green's functions" => "api/greensfunctions.md",
            "Data structures and functions for Bloch periodicity" => "api/periodicdata.md",
            "Molecular data structures and methods" => "api/moldata.md",
            "Intermolecular internuclear coupling data structures and methods" => "api/DeltaKIdata.md",
            "Transformation data structures and methods" => "api/transdata.md",
            "vdW interactions and thermal radiation functions" => "api/vdWRHT.md",
        ]
    ],
    format = Documenter.HTML(prettyurls = false) # keep this line only for local builds
);