text:
"""
    IdentityAxis(start, stop) -> axis
    IdentityAxis(keys::AbstractUnitRange) -> axis
    IdentityAxis(keys::AbstractUnitRange, indices::AbstractUnitRange) -> axis

`AbstractAxis` subtype that preserves indices after indexing.

# Examples
```julia
julia> using AxisIndices

julia> axis = AxisIndices.IdentityAxis(3:5)
idaxis(3:5)(SimpleAxis(1:3))

julia> axis[4:5]
idaxis(4:5)(SimpleAxis(2:3))

```
"""
struct IdentityAxis{I,Inds<:AbstractAxis,F} <: AbstractOffsetAxis{I,Inds,F}
    offset::F
    parent::Inds

    @inline function IdentityAxis{I,Inds,F}(f::Integer, inds::AbstractUnitRange) where {I,Inds,F}
        if inds isa Inds && f isa F
            return new{I,Inds,F}(f, inds)
        else
            return IdentityAxis(f, Inds(inds))
        end
    end

    function IdentityAxis{I,Inds}(f::Integer, inds::AbstractAxis) where {I,Inds}
        return IdentityAxis{I,Inds,typeof(f)}(f, inds)
    end
    function IdentityAxis{I,Inds}(f::Integer, inds::AbstractArray) where {I,Inds}
        return IdentityAxis{I,Inds}(f, compose_axis(inds))
    end
    function IdentityAxis{I,Inds}(f::AbstractUnitRange, inds::AbstractAxis) where {I,Inds}
        return IdentityAxis{I,Inds}(static_first(ks) - static_first(inds), inds)
    end

    function IdentityAxis{I,Inds}(f::AbstractUnitRange, inds::AbstractArray) where {I,Inds}
        return IdentityAxis{I,Inds}(f, compose_axis(inds))
    end
    function IdentityAxis{I,Inds}(ks::AbstractUnitRange) where {I,Inds}
        if can_change_size(ks)
            inds = OneToMRange(length(ks))
        else
            inds = OneTo(static_length(ks))
        end
        f = static_first(ks) - one(F)
        return IdentityAxis{eltype(inds),typeof(f),typeof(inds)}(f, inds)
    end

    # IdentityAxis{I}
    function IdentityAxis{I}(f::Integer, inds::AbstractOffsetAxis) where {I}
        return IdentityAxis(f + offsets(inds, 1), parent(inds))
    end
    function IdentityAxis{I}(f::Integer, inds::AbstractArray) where {I}
        return IdentityAxis{I}(f, compose_axis(inds))
    end
    function IdentityAxis{I}(f::Integer, inds::AbstractAxis) where {I}
        return IdentityAxis{I,typeof(inds)}(f, inds)
    end
    function IdentityAxis{I}(ks::AbstractRange{<:Integer}) where {I}
        return IdentityAxis{I}(ks, OneTo{I}(length(ks)))
    end
    function IdentityAxis{I}(start::Integer, stop::Integer) where {I}
        return IdentityAxis{I}(start:stop)
    end
    function IdentityAxis{I}(ks::AbstractRange, inds::AbstractAxis) where {I}
        return IdentityAxis{I}(static_first(ks) - static_first(inds), inds)
    end
    function IdentityAxis{I}(ks::AbstractRange, inds::AbstractOffsetAxis) where {I}
        return IdentityAxis{I}(ks, parent(inds))
    end
    function IdentityAxis{I}(ks::AbstractRange, inds::AbstractArray) where {I}
        return IdentityAxis{I}(ks, compose_axis(inds))
    end

    # IdentityAxis
    function IdentityAxis(f::Integer, inds::AbstractAxis)
        return IdentityAxis{eltype(inds)}(f, inds)
    end
    function IdentityAxis(ks::AbstractRange, inds::AbstractAxis)
        return IdentityAxis{eltype(inds)}(ks, compose_axis(inds))
    end
    function IdentityAxis(ks::AbstractRange, inds::AbstractArray)
        return IdentityAxis(ks, compose_axis(inds))
    end
    IdentityAxis(start::Integer, stop::Integer) = IdentityAxis(start:stop)
    function IdentityAxis(ks::Ks) where {Ks}
        if can_change_size(ks)
            inds = SimpleAxis(OneToMRange(length(ks)))
        else
            inds = SimpleAxis(StaticInt(1):static_length(ks))
        end
        f = static_first(ks)
        f = f - one(f)
        return new{eltype(inds),typeof(inds),typeof(f)}(f, inds)
    end
end

@inline Base.getproperty(axis::IdentityAxis, k::Symbol) = getproperty(parent(axis), k)

ArrayInterface.known_first(::Type{T}) where {T<:IdentityAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_first(::Type{T}) where {Inds,F,T<:IdentityAxis{<:Any,Inds,StaticInt{F}}}
    if known_first(Inds) === nothing
        return nothing
    else
        return known_first(Inds) + F
    end
end
Base.first(axis::IdentityAxis) = first(parent(axis)) + getfield(axis, :offset)

ArrayInterface.known_last(::Type{T}) where {T<:IdentityAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_last(::Type{T}) where {Inds,F,T<:IdentityAxis{<:Any,Inds,StaticInt{F}}}
    if known_last(Inds) === nothing
        return nothing
    else
        return known_last(Inds) + F
    end
end
Base.last(axis::IdentityAxis) = last(parent(axis)) + getfield(axis, :offset)

function ArrayInterface.unsafe_reconstruct(axis::IdentityAxis, inds; keys=nothing, kwargs...)
    if keys === nothing
        if inds isa AbstractOffsetAxis
            f_axis = offsets(axis, 1)
            f_inds = offsets(inds, 1)
            if f_axis === f_inds
                return OffsetAxis(offsets(axis, 1), unsafe_reconstruct(parent(axis), parent(inds); kwargs...))
            else
                return OffsetAxis(f_axis + f_inds, unsafe_reconstruct(parent(axis), parent(inds); kwargs...))
            end
        else
            return OffsetAxis(offsets(axis, 1), unsafe_reconstruct(parent(axis), inds; kwargs...))
        end
    else
        return IdentityAxis(keys, inds)
    end
end

"""
    idaxis(inds::AbstractUnitRange{<:Integer})

Shortcut for creating [`IdentityAxis`](@ref).

# Examples

```jldoctest
julia> using AxisIndices

julia> AxisArray(ones(3), idaxis)[2:3]
2-element AxisArray(::Vector{Float64}
  • axes:
     1 = 2:3
)
     1
  2  1.0
  3  1.0

```
"""
struct IdAxis <: AxisInitializer end
axis_method(::IdAxis, x, inds) = IdentityAxis(x, inds)
const idaxis = IdAxis()
function idaxis(collection)
    if known_step(collection) === 1
        return IdentityAxis(collection)
    else
        return idaxis(collection, axes(collection))
    end
end

"""
    IdentityArray(A::AbstractArray)

Provides [`IdentityAxis`](@ref)s for indexing `A`.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.IdentityArray(ones(3,3))[2:3, 2:3]
2×2 AxisArray(::Matrix{Float64}
  • axes:
     1 = 2:3
     2 = 2:3
)
     2    3
  2  1.0  1.0
  3  1.0  1.0

```
"""
const IdentityArray{T,N,P,A<:Tuple{Vararg{<:IdentityAxis}}} = AxisArray{T,N,P,A}

IdentityArray(A::AbstractArray{T,N}) where {T,N} = IdentityArray{T,N}(A)

IdentityArray(A::AbstractArray{T,0}) where {T} = AxisArray(A)

IdentityArray{T}(A::AbstractArray) where {T} = IdentityArray{T,ndims(A)}(A)

IdentityArray{T,N}(A::AbstractArray) where {T,N} = IdentityArray{T,N,typeof(A)}(A)

function IdentityArray{T,N,P}(x::P) where {T,N,P<:AbstractArray{T,N}}
    return initialize_axis_array(x, map(IdentityAxis, axes(x)))
end
function AxisArray{T,N,P}(x::P, axs::Tuple) where {T,N,P}
    axs = map(IdentityAxis, axs)
    return AxisArray{T,N,P}(convert(P, x), axs)
end

IdentityArray{T,N,P}(A::IdentityArray{T,N,P}) where {T,N,P} = A

"""
    IdentityVector(v::AbstractVector)

Provides an [`IdentityAxis`](@ref) for indexing `v`.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.IdentityVector(ones(4))[3:4]
2-element AxisArray(::Vector{Float64}
  • axes:
     1 = 3:4
)
     1
  3  1.0
  4  1.0

```
"""
const IdentityVector{T,P<:AbstractVector{T},Ax1<:IdentityAxis} = IdentityArray{T,1,P,Tuple{Ax1}}

IdentityVector(A::AbstractVector) = IdentityArray{eltype(A)}(A)
IdentityVector{T}(A::AbstractVector) where {T} = IdentityArray{T,1}(A)

"""
    IdentityVector{T}(init::ArrayInitializer, sz::Integer)

Creates a vector with elements of type `T` of size `sz` an [`IdentityAxis`](@ref).

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.IdentityVector{Union{Missing, Int}}(missing, 3)[2:3]
2-element AxisArray(::Vector{Union{Missing, Int64}}
  • axes:
     1 = 2:3
)
     1
  2   missing
  3   missing

```
"""
function IdentityVector{T}(init::ArrayInitializer, arg) where {T}
    return IdentityVector{T}(Vector{T}(init, arg))
end

@inline function ArrayInterface.to_axis(::IndexStyle, axis::IdentityAxis, inds)
    return unsafe_reconstruct(axis, StaticInt(1):static_length(inds); keys=inds)
end

function print_axis(io::IO, axis::IdentityAxis)
    print(io, "idaxis(")
    print(io, Int(first(axis)))
    print(io, ")(")
    print(io, Int(last(axis)))
    print(io, " parent=$(parent(axis)))")
end


text:
@testset "Step" begin
    function _getbucket(t, buckets)
        i = findlast(x -> t > x, cumsum(buckets))
        i = isnothing(i) ? 1 : i + 1

        return i
    end

    λ = rand()
    γ = rand()
    step_sizes = [rand(1:10), rand(1:10)]
    s = Step(λ = λ, γ = γ, step_sizes = step_sizes)
    @test s == Step(λ, γ, step_sizes)
    @test Step(λ, γ, 1).step_sizes == [1]
    @test all(λ == s(t) for t in 1:step_sizes[1])
    @test all(λ * γ == s(t) for t in (step_sizes[1] + 1):(step_sizes[1] + step_sizes[2] + 1))
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "Exp" begin
    λ = rand()
    γ = rand()
    s = Exp(λ = λ, γ = γ)
    @test s == Exp(λ, γ)
    @test all(λ * γ^(t - 1) == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

@testset "Poly" begin
    λ = rand()
    p = rand(1:20)
    max_iter = rand(1:100)
    s = Poly(λ = λ, p = p, max_iter = max_iter)
    @test s == Poly(λ, p, max_iter)
    @test all(λ * (1 - (t - 1) / max_iter)^p == s(t) for t in 1:max_iter)
    @test all(p == s(t) for (t, p) in zip(1:max_iter, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ)
    @test Base.IteratorSize(typeof(s)) == Base.HasLength()
    @test length(s) == max_iter
    @test_throws BoundsError s(max_iter + 1)
    @test axes(s) == 1:length(s)
end

@testset "Inv" begin
    λ = rand()
    γ = rand()
    p = rand(1:20)
    s = Inv(λ = λ, p = p, γ = γ)
    @test s == Inv(λ, γ, p)
    @test all(λ / (1 + (t - 1) * γ)^p == s(t) for t in 1:100)
    @test all(p == s(t) for (t, p) in zip(1:100, s))
    @test Base.IteratorEltype(typeof(s)) == Base.HasEltype()
    @test eltype(s) == eltype(λ)
    @test Base.IteratorSize(typeof(s)) == Base.IsInfinite()
    @test axes(s) == (OneToInf(),)
end

ext: jl


text:
using Flux, DiffEqFlux, DifferentialEquations, Plots

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.0f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
mp = Float32[0.2,0.2]
function true_noise_func(du,u,p,t)
    du .= mp.*u
end
prob = SDEProblem(trueODEfunc,true_noise_func,u0,tspan)
sde_data = Array(solve(prob,SOSRI(),saveat=t))

dudt = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
ps = Flux.params(dudt)
n_sde = x->neural_dmsde(dudt,x,mp,tspan,SOSRI(),saveat=t,reltol=1e-2,abstol=1e-2)

pred = n_sde(u0) # Get the prediction using the correct initial condition
scatter(t,sde_data[1,:],label="data")
scatter!(t,Flux.data(pred[1,:]),label="prediction")

function predict_n_sde()
  n_sde(u0)
end
loss_n_sde() = sum(abs2,sde_data .- predict_n_sde())
Flux.back!(loss_n_sde())

data = Iterators.repeated((), 100)
opt = ADAM(0.025)
cb = function () #callback function to observe training
  sample = predict_n_sde()
  # loss against current data
  display(sum(abs2,sde_data .- sample))
  # plot current prediction against data
  cur_pred = Flux.data(sample)
  pl = scatter(t,sde_data[1,:],label="data")
  scatter!(pl,t,cur_pred[1,:],label="prediction")
  display(plot(pl))
end

# Display the SDE with the initial parameter values.
cb()

Flux.train!(loss_n_sde, ps, data, opt, cb = cb)

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using DataFrames
using Query
using PyPlot
using LaTeXStrings

include("schema.jl")

const STDERR_ALPHA = 0.3
const PLOT_DIR="plots"
const DATA_DIR = "data"
const WIDTH = 3
const HEIGHT = 3
const SHARED_YLABEL = "AIDE estimate (nats)"
const SHARED_YLIM = [-0.25, 9]
const TITLE_FONTSIZE=12

function add_baseline()
    ax = plt[:gca]()
    ax[:set_xscale]("log")
    xlim = ax[:get_xlim]()
    plt[:plot](xlim, [0, 0], color="black", alpha=0.5)
    ax[:set_xlim](xlim)
end

function plot_smc_data(data::DataFrame)
    num_rejuvenation = 2
    for k in unique(data[COL_NUM_METAINFERENCE])
        subdata  = @from i in data begin
            @where i.num_metainference == k
            @select {i.num_particles,
                     i.aide_estimate,
                     i.aide_estimate_standard_error}
            @collect DataFrame
        end
        exponent = Int(log10(k))
        num_particles = subdata[COL_NUM_PARTICLES]
        est = subdata[COL_AIDE_ESTIMATE]
        err = subdata[COL_AIDE_STDERR]
        label=latexstring("\$M_t = 10^$exponent\$")
        plt[:plot](num_particles, est, label=label)
        plt[:fill_between](num_particles, est-err, est+err, alpha=STDERR_ALPHA)
    end
    add_baseline()
    plt[:gca]()[:set_ylim](SHARED_YLIM)
    plt[:xlabel]("Number of particles")
    plt[:ylabel](SHARED_YLABEL)
    plt[:legend]()
    plt[:title]("Sequential Monte Carlo", fontsize=TITLE_FONTSIZE)
end

function plot_bbvi_data(data::DataFrame)
    est = data[COL_AIDE_ESTIMATE]
    err = data[COL_AIDE_STDERR]
    num_vi_iters = data[COL_NUM_SGD_ITERS]
    plt[:plot](num_vi_iters, est, label="Variational")
    kl_asymptote = est[end] # assumes sorted.
    plt[:fill_between](num_vi_iters, est - err, est + err, alpha=STDERR_ALPHA)
    add_baseline()
    xlim = [minimum(num_vi_iters), maximum(num_vi_iters)]
    plt[:plot](xlim, [kl_asymptote, kl_asymptote], color="black", alpha=0.5)
    plt[:gca]()[:set_xlim](xlim)
    plt[:gca]()[:set_ylim](SHARED_YLIM)
    plt[:xlabel]("Number of gradient steps")
    plt[:ylabel](SHARED_YLABEL)
    plt[:title]("Variational Inference", fontsize=TITLE_FONTSIZE)
end

function plot_mcmc_data(data::DataFrame)
    for k in unique(data[COL_NUM_METAINFERENCE])
        subdata  = @from i in data begin
            @where i.num_metainference == k
            @select {i.num_mcmc_iters,
                     i.aide_estimate,
                     i.aide_estimate_standard_error}
            @collect DataFrame
        end
        exponent = Int(log10(k))
        num_mcmc_iters = subdata[COL_NUM_MCMC_ITERS]
        est = subdata[COL_AIDE_ESTIMATE]
        err = subdata[COL_AIDE_STDERR]
        label=latexstring("\$M_t = 10^$exponent\$")
        plt[:plot](1 + num_mcmc_iters, est, label=label)
        plt[:fill_between](1 + num_mcmc_iters, est-err, est+err, alpha=STDERR_ALPHA)
    end
    add_baseline()
    plt[:gca]()[:set_ylim](SHARED_YLIM)
    plt[:xlabel]("1 + Number of transitions")
    plt[:ylabel](SHARED_YLABEL)
    plt[:title]("Metropolis-Hastings", fontsize=TITLE_FONTSIZE)
end

println("plotting divergence bounds..")

plt[:figure](figsize=(WIDTH*3,HEIGHT))

println("plotting SMC data...")
smc_data = readtable("$DATA_DIR/smc_data.csv")
plt[:subplot](1, 3, 1)
plot_smc_data(smc_data)

println("plotting MCMC data...")
mcmc_data = readtable("$DATA_DIR/mcmc_data.csv")
plt[:subplot](1, 3, 2)
plot_mcmc_data(mcmc_data)

println("plotting BBVI data...")
bbvi_data = readtable("$DATA_DIR/bbvi_data.csv")
plt[:subplot](1, 3, 3)
plot_bbvi_data(bbvi_data)

println("Saving output plot.")
plt[:tight_layout]()
plt[:savefig]("$PLOT_DIR/linear_regression_combined.pdf")

println("done!")


text:
using PkgServer, Pkg, HTTP, JSON
using Test

# You can either perform the following setup:
#  - Already-running PkgServer, located at $JULIA_PKG_SERVER
#  - PkgServer cache directory located at $PKG_SERVER_CACHE_DIR
#
# Or you can leave those blank, and we'll start up a PkgServer for you, running on a background task.

if VERSION < v"1.4"
    error("These tests require a PkgServer-compatible Julia version!")
end

# If these are not set, we will attempt to auto-initiate them.
server_process = nothing
if isempty(get(ENV, "JULIA_PKG_SERVER", "")) || isempty(get(ENV, "JULIA_PKG_SERVER_CACHE_DIR", ""))
    @info("Automatically starting local PkgServer for testing at http://127.0.0.1:8000")
    # Start a background PkgServer as a separate process
    code_dir = dirname(@__DIR__)
    temp_dir = mktempdir(; cleanup=false)
    cd(temp_dir) do
        global server_process = run(`$(Base.julia_cmd()) --project=$(code_dir) $(code_dir)/bin/run_server.jl`; wait=false)

        ENV["JULIA_PKG_SERVER"] = "http://127.0.0.1:8000"
        ENV["JULIA_PKG_SERVER_CACHE_DIR"] = joinpath(temp_dir, "cache")

        # Tell Julia to cleanup this directory at process exit
        Base.Filesystem.temp_cleanup_later(temp_dir)
    end
end

server_url = ENV["JULIA_PKG_SERVER"]
cache_dir = ENV["JULIA_PKG_SERVER_CACHE_DIR"]

# Server may be sleepy and still waking up; wait until we can ask it for `/registries`
t_start = time()
@info("Waiting for PkgServer liveness")
while true
    # Try to get an HTTP 200 OK on /registries
    response_code = try
        HTTP.get("$(server_url)/registries"; retry = false, readtimeout=1).status
    catch e
        # If it just can't connect, deal with it silently, otherwise rethrow
        if isa(e, HTTP.IOExtras.IOError) && e.e.code in (-Base.Libc.ECONNREFUSED, -Base.Libc.EPIPE, -Base.Libc.ECONNRESET)
            # Return fake status code, this doesn't really matter, just so long as it's not 200
            e.e.code
        else
            rethrow(e)
        end
    end

    if response_code == 200
        break
    end

    # If we've been trying this for more than 20s, error out
    if (time() - t_start) >= 20
        error("Unable to hit testing server at $(server_url)/registries, get HTTP $(response_code)")
    end

    # Sleep a bit between attempts
    sleep(0.2)
end

@testset "Direct HTTP requests" begin
    # Test that we get a sensical answer for /registries
    response = HTTP.get("$(server_url)/registries")
    @test response.status == 200
    registry_url = chomp(String(response.body))
    @test match(r"^/registry/([a-z0-9]+-){4}[a-z0-9]+/[a-z0-9]+$", registry_url) != nothing
    registry_uuid, registry_treehash = split(registry_url, "/")[3:4]

    # Test asking for that registry directly, unpacking it and verifying the treehash
    mktemp() do tarball_path, tarball_io
        response = HTTP.get("$(server_url)/registry/$(registry_uuid)/$(registry_treehash)"; response_stream=tarball_io)
        close(tarball_io)
        @test response.status == 200
        @test registry_treehash == PkgServer.tarball_git_hash(tarball_path)
    end

    # Verify that these files exist within the cache
    @test isfile(joinpath(cache_dir, "registries"))
    @test isfile(joinpath(cache_dir, "registry", registry_uuid, registry_treehash))

    # Next, hit the `/meta` endpoint, ensure that the version it reports matches with what we expect:
    response = HTTP.get("$(server_url)/meta")
    @test response.status == 200
    meta = JSON.parse(String(response.body))
    @test haskey(meta, "pkgserver_version")
    @test meta["pkgserver_version"] == PkgServer.get_pkgserver_version()
end

function with_depot_path(f::Function, dp::Vector{String})
    OLD_DEPOT_PATH = copy(Base.DEPOT_PATH)
    empty!(Base.DEPOT_PATH)
    append!(Base.DEPOT_PATH, dp)
    try
        f()
    finally
        empty!(Base.DEPOT_PATH)
        append!(Base.DEPOT_PATH, OLD_DEPOT_PATH)
    end
end

@testset "Package Installation" begin
    mktempdir() do temp_dir
        with_depot_path([temp_dir]) do
            # Prevent `Pkg.add()` from changing PkgServer.jl's Project.toml
            Pkg.activate(temp_dir)
            # Install something with a huge number of Package and Artifact dependencies
            Pkg.add(Pkg.PackageSpec(;name="Gtk", version=v"1.1.2"))

            # Install something with platform-independent artifacts, so that we can check the hashes
            Pkg.add(Pkg.PackageSpec(;name="DotNET", version=v"0.1.0"))
            Pkg.add(Pkg.PackageSpec(;name="FIGlet", version=v"0.2.1"))
        end
    end

    # Test that the resources we expect to exist are in fact cached on the server.
    # We just test that a few UUIDs and tree hashes exist:
    pkg_uuid_treehashes = [
        # XML2_jll
        "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" => "ed5603a695aefe3e9e404fc7b052e02cc72cfab6",
        # Gtk
        "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" => "ab2676a65345f8813be07675aa657464f98e6704",
        # FixedPointNumbers
        "53c48c17-4a7d-5ca2-90c5-79b7896eea93" => "4aaea64dd0c30ad79037084f8ca2b94348e65eaa",
        # Zlib_jll
        "83775a58-1f1d-513f-b197-d71354ab007a" => "5618a43055eb09377edca21d19d0e99bce24a9c3",
    ]
    for (uuid, treehash) in pkg_uuid_treehashes
        @test isfile(joinpath(cache_dir, "package", uuid, treehash))
    end

    artifact_treehashes = [
        # DotNET clrbridge
        "5c005e142ebba033996dcc97b249e9e5ebcbf138",
        # FIGlet fonts
        "125ac0315d68bbb612f8c2189ea83401f73238f0",
    ]
    for treehash in artifact_treehashes
        @test isfile(joinpath(cache_dir, "artifact", treehash))
    end

    # Ensure that, when we hit `/meta` now, the server knows that it has a bunch of packages and artifacts:
    response = HTTP.get("$(server_url)/meta")
    @test response.status == 200
    meta = JSON.parse(String(response.body))
    @test haskey(meta, "packages_cached")
    @test meta["packages_cached"] >= 80
    @test haskey(meta, "artifacts_cached")
    @test meta["artifacts_cached"] >= 30
end

if server_process != nothing
    @info("Reaping automatically-started local PkgServer...")
    kill(server_process)
    wait(server_process)
end


text:
@testset "Equality of systems" begin
    A = [1. 1; 1 -1]
    A2 = [1.0 1+1e-8; 1 -1]
    c = [1.; 1.]
    X = Line([1., -1], 0.) # line x = y
    Y = Singleton([0., 0.])
    s1 = ConstrainedAffineDiscreteSystem(A, c, X)
    s2 = ConstrainedAffineContinuousSystem(A, c, X)
    s3 = ConstrainedAffineDiscreteSystem(A2, c, X)
    s4 = ConstrainedAffineDiscreteSystem(A, c, Y)

    @test s1 == s1
    @test s1 != s2
    @test s1 != s3
    @test s1 != s4
    @test s1 ≈ s1
    @test !(s1 ≈ s2)
    @test s1 ≈ s3
    @test !(s1 ≈ s4)

    x0 = [1.0, 3.0]
    ivp1 = InitialValueProblem(s1, x0)
    ivp2 = InitialValueProblem(s2, x0)
    ivp3 = InitialValueProblem(s3, x0)
    ivp4 = InitialValueProblem(s4, x0)

    @test ivp1 == ivp1
    @test ivp1 != ivp2
    @test ivp1 != ivp3
    @test ivp1 != ivp4
    @test ivp1 ≈ ivp1
    @test !(ivp1 ≈ ivp2)
    @test ivp1 ≈ ivp3
    @test !(ivp1 ≈ ivp4)
end


text:
using Test
using QuantLib
using Dates

euribor3m = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(3)))

todaysDate = Date(2006, 5, 23)
set_eval_date!(settings, todaysDate)
calendar = euribor3m.fixingCalendar
fixingDays = euribor3m.fixingDays
settlementDate = QuantLib.Time.advance(Dates.Day(fixingDays), calendar, todaysDate)

threeMonthFraQuote = zeros(10)
threeMonthFraQuote[2] = 0.030
threeMonthFraQuote[3] = 0.031
threeMonthFraQuote[4] = 0.032
threeMonthFraQuote[7] = 0.033
threeMonthFraQuote[10] = 0.034

## QUOTES ##
h1x4 = Quote(threeMonthFraQuote[2])
h2x5 = Quote(threeMonthFraQuote[3])
h3x6 = Quote(threeMonthFraQuote[4])
h6x9 = Quote(threeMonthFraQuote[7])
h9x12 = Quote(threeMonthFraQuote[10])

## RATE HELPERS ##
fraDayCount = euribor3m.dc
convention = euribor3m.convention
endOfMonth = euribor3m.endOfMonth

fra1x4 = FraRateHelper(h1x4, 1, 4, fixingDays, calendar, convention, endOfMonth, fraDayCount)
fra2x5 = FraRateHelper(h2x5, 2, 5, fixingDays, calendar, convention, endOfMonth, fraDayCount)
fra3x6 = FraRateHelper(h3x6, 3, 6, fixingDays, calendar, convention, endOfMonth, fraDayCount)
fra6x9 = FraRateHelper(h6x9, 6, 9, fixingDays, calendar, convention, endOfMonth, fraDayCount)
fra9x12 = FraRateHelper(h9x12, 9, 12, fixingDays, calendar, convention, endOfMonth, fraDayCount)

fraInstruments = FraRateHelper[fra1x4, fra2x5, fra3x6, fra6x9, fra9x12]
termStructureDC = QuantLib.Time.ISDAActualActual()
tol = 1.0e-15

fraTermStructure = PiecewiseYieldCurve(settlementDate, fraInstruments, termStructureDC, QuantLib.Math.LogLinear(), Discount(), tol, IterativeBootstrap())

euribor3m = update_termstructure(euribor3m, fraTermStructure)

fraCalendar = euribor3m.fixingCalendar
fraBusinessDayConvention = euribor3m.convention
fraFwdType = LongPosition()
fraNotional = 100.0
monthsToStart = [1, 2, 3, 6, 9]
fraTermMonths = 3
fraValueDate = QuantLib.Time.advance(Dates.Month(monthsToStart[1]), fraCalendar, settlementDate, fraBusinessDayConvention)
fraMaturityDate = QuantLib.Time.advance(Dates.Month(fraTermMonths), fraCalendar, fraValueDate, fraBusinessDayConvention)
fraStrikeRate = threeMonthFraQuote[monthsToStart[1] + 1]
myFRA = ForwardRateAgreement(fraValueDate, fraMaturityDate, fraFwdType, fraStrikeRate, fraNotional, euribor3m, fraTermStructure)

# test calculations
@test forward_rate(myFRA).rate == 0.030000000000000172
@test spot_value(myFRA) == 99.73470289713767
@test forward_value(myFRA) == 100.76666666666667
@test implied_yield(myFRA, spot_value(myFRA), forward_value(myFRA), settlementDate, SimpleCompounding(), fraDayCount).rate == 0.030039933543589872

# other tests
@test maturity_date(fra1x4) == Date(2006, 9, 26)
@test QuantLib.implied_quote(fraTermStructure.instruments[1]) == 0.030000000000000172


text:
function samplehull_discretizer(zs::AbstractVector{<:AbstractFloat},
                                algo::Discretizers.DiscretizationAlgorithm = DiscretizeUniformWidth(:auto);
                                closed = :right)
    grid = binedges(algo, zs)
    if closed === :right
        grid[1] = prevfloat(grid[1])
        return BoundedIntervalDiscretizer{:open,:closed}(grid)
    elseif closed === :left
        grid[end] = nextfloat(grid[end])
        return BoundedIntervalDiscretizer{:closed,:open}(grid)
    else
        throw(ArgumentError("closed can be :right or :left"))
    end
end

#function support_discretizer(Zs::AbstractVector, )
#    get_nbins
#end


text:
@testset "FlexTable" begin
    t = @inferred(FlexTable(a = [1,2,3], b = [2.0, 4.0, 6.0]))::FlexTable

    @test axes(t) === (Base.OneTo(3),)
    @test size(t) === (3,)
    @test length(t) === 3
    @test t[2] === (a = 2, b = 4.0)
    @test @inferred(t[1:2])::FlexTable == FlexTable(a = [1,2], b = [2.0,4.0])
    @test @inferred(t[:])::FlexTable == t
    @test @inferred(view(t, 1:2))::FlexTable == FlexTable(a = [1,2], b = [2.0,4.0])
    @test @inferred(view(t, :))::FlexTable == t
    @test_throws BoundsError t[4]

    @test similar(t) isa typeof(t)
    @test axes(similar(t)) === axes(t)

    @test @inferred(vcat(t, t))::FlexTable{1} == FlexTable(a = [1,2,3,1,2,3], b = [2.0, 4.0, 6.0, 2.0, 4.0, 6.0])
    @test @inferred(hcat(t, t))::FlexTable{2} == FlexTable(a = [1 1;2 2;3 3], b = [2.0 2.0; 4.0 4.0; 6.0 6.0])
    # TODO hvcat

    t2 = empty(t)
    @test t2 isa typeof(t)
    @test isempty(t2)

    t3 = copy(t)
    @test t3 isa typeof(t)
    @test t3 == t
    @test @inferred(t3[1:2])::FlexTable == FlexTable(a = [1,2], b = [2.0,4.0])
    t3[3] = (a = 4, b = 8.0)
    @test t3 == FlexTable(a = [1,2,4], b = [2.0, 4.0, 8.0])
    t3[1:2] .= Ref((a = 0, b = 0.0))
    @test t3 == FlexTable(a = [0,0,4], b = [0.0, 0.0, 8.0])
    t3[1:2] = FlexTable(a = [2,3], b = [4.0, 6.0])
    @test t3 == FlexTable(a = [2,3,4], b = [4.0, 6.0, 8.0])
    v3 = view(t3, :)
    v3[3] = (a = 5, b = 10.0)
    @test t3 == FlexTable(a = [2,3,5], b = [4.0, 6.0, 10.0])

    empty!(t3)
    @test isempty(t3)

    push!(t3, (a=10, b=0.0))
    @test t3 == FlexTable(a = [10], b = [0.0])

    pushfirst!(t3, (a=-10, b=-10.0))
    @test t3 == FlexTable(a = [-10, 10], b = [-10.0, 0.0])

    insert!(t3, 2, (a=5,b=5.0))
    @test t3 == FlexTable(a = [-10, 5, 10], b = [-10.0, 5.0, 0.0])

    @test splice!(t3, 2:1) == empty(t)
    @test t3 == FlexTable(a = [-10, 5, 10], b = [-10.0, 5.0, 0.0])
    @test splice!(t3, 2:1, (a=1, b=1.0)) == empty(t)
    @test t3 == FlexTable(a = [-10, 1, 5, 10], b = [-10.0, 1.0, 5.0, 0.0])
    @test splice!(t3, 3:2, (a=[2], b=[2.0])) == empty(t)
    @test t3 == FlexTable(a = [-10, 1, 2, 5, 10], b = [-10.0, 1.0, 2.0, 5.0, 0.0])
    @test splice!(t3, 2) === (a = 1, b = 1.0)
    @test t3 == FlexTable(a = [-10, 2, 5, 10], b = [-10.0, 2.0, 5.0, 0.0])
    @test splice!(t3, 2, (a = 3, b = 3.0)) === (a = 2, b = 2.0)
    @test t3 == FlexTable(a = [-10, 3, 5, 10], b = [-10.0, 3.0, 5.0, 0.0])
    @test splice!(t3, 2, (a = [4], b = [4.0])) === (a = 3, b = 3.0)
    @test t3 == FlexTable(a = [-10, 4, 5, 10], b = [-10.0, 4.0, 5.0, 0.0])

    @test pop!(t3) === (a = 10, b = 0.0)
    @test t3 == FlexTable(a = [-10, 4, 5], b = [-10.0, 4.0, 5.0])

    @test popfirst!(t3) === (a = -10, b = -10.0)
    @test t3 == FlexTable(a = [4, 5], b = [4.0, 5.0])

    append!(t3, FlexTable(a = [6], b = [6.0]))
    @test t3 == FlexTable(a = [4,5,6], b = [4.0, 5.0, 6.0])

    prepend!(t3, FlexTable(a = [3], b = [3.0]))
    @test t3 == FlexTable(a = [3,4,5,6], b = [3.0, 4.0, 5.0, 6.0])

    @testset "Merging FlexTables" begin
        t1 = FlexTable(a = [1,2,3],)
        t2 = FlexTable(b = [2.0, 4.0, 6.0],)
        t3 = FlexTable(a = [1,2,3], b = [2.0, 4.0, 6.0])

        @test @inferred(map(merge, t1, t2))::FlexTable == t3
        @test @inferred(mapview(merge, t1, t2))::FlexTable == t3
        @test @inferred(broadcast(merge, t1, t2))::FlexTable == t3
    end

    @testset "GetProperty on FlexTables" begin
        t = FlexTable(a = [1,2,3], b = [2.0, 4.0, 6.0])

        @test map(getproperty(:a), t)::Vector == [1,2,3]
        @test mapview(getproperty(:a), t)::Vector == [1,2,3]
        @test broadcast(getproperty(:a), t)::Vector == [1,2,3]
    end

    # setproperty!
    t4 = FlexTable(a = [1,2,3])
    t4.b = [2.0, 4.0, 6.0]
    @test t4 == t
end


text:
# The intakes bounds of the network are determined by the
# medium concentration in the Chemostat model (see Cossios paper)
# This is a base medium for modeling
base_intake_info = Dict(
    "EX_glc__D_e" => Dict("c"=>  15.0, "lb"=> -100.0),
    "EX_h2o_e" => Dict("c"=>  99999.0, "lb"=> -100.0),
    "EX_nh4_e" => Dict("c"=>  99999.0, "lb"=> -100.0),
    "EX_o2_e" => Dict("c"=>  99999.0, "lb"=> -100.0),
    "EX_pi_e" => Dict("c"=>  99999.0, "lb"=> -100.0),
    # "EX_acald_e" => Dict("c"=>  99999.0, "lb"=> -100.0),
)


text:
using BinaryBuilder

version = v"0.1.1-dev"

# Collection of sources required to build libreadstat
sources = [
    "https://github.com/WizardMac/ReadStat.git" =>
    "93a6349ad402faf59547f9d80f21a7955585ce01",
]

script = raw"""
cd $WORKSPACE/srcdir
cd ReadStat/
./autogen.sh
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then ./configure --prefix=${prefix} --host=${target} CFLAGS="-I$prefix/include" LDFLAGS="-L$prefix/lib"; else ./configure --prefix=${prefix} --host=${target}; fi
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms()

dependencies = [
    "https://github.com/davidanthoff/IConvBuilder/releases/download/v1.15%2Bbuild.3/build.jl"
]

products = prefix -> [
    ExecutableProduct(prefix,"readstat", :readstat),
    LibraryProduct(prefix,"libreadstat", :libreadstat)
]

build_tarballs(ARGS, "ReadStat", version, sources, script, platforms, products, dependencies)


text:
# # Basic Tutorial

# Since DFControl is aimed at improving the day to day quality of life of a material's scientist/anyone running DFT codes,
# we will look at a simple demonstration of how by creating and submitting some
# [Quantum-Espresso](https://www.quantum-espresso.org/) calculations on Si
# starting from a cif file specifying the structure.

using DFControl

# First we download the cif file, extract the `Structure` and assign the right pseudos to it.
# In this case Si (F d -3 m :1) from http://www.crystallography.net/cod/9011998.cif
using Downloads
cif_file = Downloads.download("http://www.crystallography.net/cod/9011998.cif", "Si.cif")

structure = Structure(cif_file; name = "Si")
set_pseudos!(structure, :pbesol)

# This assumes that the `:pbesol` pseudopotential set was installed during the
#md # [configuration step](@ref Configuration).

# Next we specify the executables with which to run the QE calculations.
# We assume here that `mpirun` is installed and in the user's PATH, and that
# QE is installed, and to be found in the `/opt/qe/bin`, change this according to your own
# setup.

pw_execs = [Exec("mpirun", "", :np => 4), Exec("pw.x", "/opt/qe/bin/", :nk => 4)]

# Additional executable flags can be passed as varargs to the constructor of `Exec`,
# e.g. `Exec("pw.x", "/opt/qe/bin/", :nk => 4, :ndiag => 2)`.

# Then we create the first calculation for our job, we name it scf, which will be used to reference it later.
# We also pass the executables to be used and additional flags to be set to the constructor.
# Afterwards we set the kpoints to be used in the scf calculation.
scf_calculation = Calculation{QE}("scf", :calculation => "scf"; execs = pw_execs)
set_kpoints!(scf_calculation, (6, 6, 6, 1, 1, 1))

# The code recognizes internally that this 6-Tuple corresponds to a
# `K_POINTS (automatic)` block in QE. Alternatively (leading to an identical final result):

scf_calculation = Calculation{QE}("scf", :calculation => "scf"; execs = pw_execs,
                                  data = [InputData(:k_points, :automatic,
                                                    (6, 6, 6, 1, 1, 1))])

# We can now define our job:
job = DFJob("Si", structure, [scf_calculation], :ecutwfc => 20, :conv_thr => 1e-6;
            local_dir = "job")

# Additional calculations would be be added to the list `[scf_calculation]`.
# The flag => value pairs will set the specified flags to that value for all calculations in the job
# that allow recognize that flag, so it's ideal for things like cutoffs and smearing etc.

# We are now ready to submit the job, which will run in the current working directory
if false #hide
    submit(job)
else #hide
    global job = DFJob(joinpath(@__DIR__, "../../src/assets/job/"))#hide
    pop!(job) #hide
end #hide

# This will generate and save all calculation files, and the corresponding job script (`job.tt`),
# and subsequently run the job.
# First submission through `sbatch job.tt` will be tried, if that fails then the script will run
# through `bash job.tt`.

# After the job finishes the outputs can be parsed through
outputdata(job)
# or for a specific calculation
outputdata(job["scf"])

# This also demonstrates how a calculation can be referenced using its name
# (remember that we named the calculation "scf" when creating it).

# Now that the scf calculation finished succesfully, the next step is usually to
# have a look at the bandstructure. For this we generate a bands calculation,
# using the scf calculation as a template and a generated high symmetry k-point path
# with 20 kpoints per segment.

bands_calc = gencalc_bands(job["scf"], high_symmetry_kpath(job.structure, 20))

# Observe the :calculation => "bands", and automatic setting of the :verbosity => "high" flags.
# We now push! this calculation to the job queue
push!(job, bands_calc)

# However, since we know the scf succeeded there is no need to rerun it.
# To un-schedule it we do
job["scf"].run = false

# Printing the job will now highlight the scheduled calculations differently from the non-scheduled ones
job

# Seeing that all is right we submit the job again
set_localdir!(job, "job"); #hide
if false #hide
    submit(job)
else #hide
    global job = DFJob(joinpath(pathof(DFControl), "..", "..", "docs", "src", "assets",
                                "job"))#hide
end #hide

# We can access the bands through
bands = readbands(job);
# or
bands = outputdata(job["bands"])[:bands];

# They can be plotted too
fermi = readfermi(job) # = outputdata(job["scf"])[:fermi]
using Plots
plot(bands; fermi = fermi)

# Since more info (such as the structure) is available in the job,
# plotting the job leads to a richer plot
plot(job, -10, 1)

# As can be seen in the [Advanced Usage](@ref), additional information
# generated by additional calculations will be picked up by DFControl
# in order to create richer plots.


text:
using AutoGrad
using Knet

include(Knet.dir("data", "housing.jl"))

predict(w,x) = w[1]*x .+ w[2]

loss(w,x,y) = mean(abs2, y - predict(w,x))

lossgradient = grad(loss)

function train(w, data; lr=0.1)
    for (x,y) in data
        dw = lossgradient(w,x,y)
        for i = 1:length(w) # 0 -> weights, 1 -> bias
            w[i] = w[i] - lr*dw[i]
        end
    end
    return w
end

x,y = housing()
w = Any[0.1 * randn(1,13), 0.0]
for i = 1:10
    train(w, [(x,y)])
    println(loss(w,x,y))
end

max_stars_repo_name: jeannekamikaze/ann

max_issues_repo_name: jeannekamikaze/ann

max_forks_repo_name: jeannekamikaze/ann


text:
module Checkers

using Reexport
@reexport using Test
import Test: Result, Pass, Fail, Error, record, get_testset

export @test_formany, @test_cases, @test_forall, @test_exists

for (dir, filename) in [
		(".", "test-formany.jl"),
		(".", "test-cases.jl"),
		(".", "test-exists.jl"),
		(".", "test-forall.jl"),
		(".", "parse-argument-data.jl"),
		(".", "custom-generator.jl")
	]

include(joinpath(dir, filename))
end

end #of module


text:
include("../test.jl")
using GerryChain
using Test

@testset "Loading a partition" begin
    graph = load_graph_from_edge_list("data/2x2_grid.dat", "data/2x2_grid_weight.dat")
    @testset "two districts" begin
        partition = load_partition_from_line(graph,"0 1 1 1")
        @test partition.num_dists == 2
        @test partition.num_cut_edges == 2
        @test partition.assignments == [1, 2, 2, 2]
        @test partition.dist_populations == [1, 14]
    end
    @testset "three districts" begin
        partition = load_partition_from_line(graph,"0 2 1 1")
        @test partition.num_dists == 3
        @test partition.num_cut_edges == 3
        @test partition.assignments == [1, 3, 2, 2]
        @test partition.dist_populations == [1, 12, 2]
    end
    @testset "reassign district id" begin
        partition = load_partition_from_line(graph,"5 6 5 5")
        @test partition.num_dists == 2
        @test partition.num_cut_edges == 2
        @test partition.assignments == [1, 2, 1, 1]
    end
end

@testset "Calculate benchmark" begin
    graph = load_graph_from_edge_list("data/2x2_grid.dat", "data/2x2_grid_weight.dat")
    result = calculate_benchmark(graph, "data/2x2_grid_enum.txt", "")
    @testset "cut edges" begin
        @test result == Dict{Int64,Int64}(2 => 4)
    end
end

max_stars_repo_name: mggg/benchmarks

max_issues_repo_name: mggg/benchmarks

max_forks_repo_name: mggg/benchmarks


text:
using Distributed
addprocs(exeflags="--project=@.")
using Dates
using VITA2D

println("Loaded modules ...")

TASK_ID = parse(Int, ARGS[1])

L = 12
sgd_params = VITA2D.SGDParams(lr=0.2, N_mc=10^3, N_sgd=2000)

repeats = [i for i in 1:2]
# Ps = [2]
Ps = [5, 10]
hs = [2.0]#3.0, 2.0, 2.5, 4.0, 3.5]

list = Iterators.product(Ps, hs, repeats) |> collect
println("length of array jobs is:\t", length(list))

# println(list)
println(list[TASK_ID])

P, hX, r = list[TASK_ID]
path = "../data/" * string(Dates.today())
path *= "/L-$L-p-$P-h-$hX/"
mkpath(path)
path *= "iteration-$r.txt"
touch(path)

θ = P * rand(3P)
# if hX == 2.0
    # if repeats == 1
        # θ = [0.613, 1.632, 0.0, 0.3469, 0.20271, 0.0]
    # else
        # θ = [1.22, 1.632, 0.0, 0.3469, 0.20271, 0.0]
    # end
# end

# Precompile
sgd(L=L,
    P=P,
    hX = hX,
    θ = θ,
    SGD=SGDParams(),
    path=false,
    print=true)

sgd(L=L,
    P=P,
    hX = hX,
    θ = θ,
    SGD=sgd_params,
    path=path,
    print=true)


text:
"""
Module to define an interactive player. (Intended for testing purposes.)
"""
module InteractivePlay

import PlayingCards: Card, createdeck

import Cribbage: Player
import Cribbage.RandomPlay: discardtocrib, pickcardtoplay
import Cribbage.SmartPlay: averagescore, DECK
import Cribbage.Rules: cardvalue, validhands, validplays
import Cribbage.Scoring: scorehand, scoreplay

struct InteractivePlayer <: Player
    name :: String
end

function InteractivePlayer()
    return InteractivePlayer(UUIDs.uuid4())
end

function discardtocrib(player::InteractivePlayer, cards::Tuple{Vararg{Card}})
    println("You were dealt: $(string(cards))")
    subsets = validhands(cards)
    others = setdiff(DECK, cards)
    sort!(subsets, by=x -> averagescore(x, others), rev=true)
    for (i, h) in enumerate(subsets)
        if i ≥ 10
            break
        end
        println("[$i] $(string(h))")
    end
    println("Which hand would you like to keep?"); userinput = readline()
    handint = parse(Int, userinput)
    hand = subsets[handint]
    crib = setdiff(cards, hand)
    @debug "keeping $(string(hand)), discarding $(string(crib))"
    return hand, crib
end

function pickcardtoplay(player::InteractivePlayer, prior::Tuple{Vararg{Card}}, cards::Tuple{Vararg{Card}})::Union{Missing,Card}
    if length(prior) == 0
        println("No prior cards")
    else
        println("Prior cards $(sum(cardvalue(c) for c in prior)): $(string(prior))")
    end
    println("You have: $(string(cards))")
    possible_plays = validplays(prior, cards)
    if length(possible_plays) == 0
        println("You have no possible plays.")
        return missing
    end
    if length(prior) == 0
        priorarray = Card[]
    else
        priorarray = collect(prior)
    end
    sort!(possible_plays, by=x -> scoreplay(priorarray, x), rev=true)
    for (i, c) in enumerate(possible_plays)
        if i ≥ 10
            break
        end
        println("[$i] $(string(c))")
    end
    println("Which card would you like to play?"); userinput = readline()
    cardint = parse(Int, userinput)
    return possible_plays[cardint]
end

end  # module InteractivePlay


text:
# ---
# title: 1505. Minimum Possible Integer After at Most K Adjacent Swaps On Digits
# id: problem1505
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-possible-integer-after-at-most-k-adjacent-swaps-on-digits/description/>
# hidden: true
# ---
#
# Given a string `num` representing **the digits** of a very large integer and
# an integer `k`.
#
# You are allowed to swap any two adjacent digits of the integer **at most** `k`
# times.
#
# Return _the minimum integer_ you can obtain also as a string.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/06/17/q4_1.jpg)
#
#
#
#     Input: num = "4321", k = 4
#     Output: "1342"
#     Explanation: The steps to obtain the minimum integer from 4321 with 4 adjacent swaps are shown.
#
#
# **Example 2:**
#
#
#
#     Input: num = "100", k = 1
#     Output: "010"
#     Explanation: It's ok for the output to have leading zeros, but the input is guaranteed not to have any leading zeros.
#
#
# **Example 3:**
#
#
#
#     Input: num = "36789", k = 1000
#     Output: "36789"
#     Explanation: We can keep the number without any swaps.
#
#
# **Example 4:**
#
#
#
#     Input: num = "22", k = 22
#     Output: "22"
#
#
# **Example 5:**
#
#
#
#     Input: num = "9438957234785635408", k = 23
#     Output: "0345989723478563548"
#
#
#
#
# **Constraints:**
#
#   * `1 <= num.length <= 30000`
#   * `num` contains **digits** only and doesn't have **leading zeros**.
#   * `1 <= k <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# push!(LOAD_PATH, "./")
# push!(LOAD_PATH, "../")

using Distributions

### For building the tensorflow graph

const Q_SCOPE = "active_q"
const TARGET_Q_SCOPE = "target_q"

"""
    TrainGraph
type to store all the tensor of the DQN graph including both the Q network
and the training operations
"""
mutable struct TrainGraph
    sess::Session
    s::Tensor
    a::Tensor
    sp::Tensor
    r::Tensor
    done_mask::Tensor
    importance_weights::Tensor
    q_lo::Tensor
    q_lo_p::Tensor
    q::Tensor
    qp::Tensor
    target_q::Tensor
    loss::Tensor
    td_errors::Tensor
    train_op::Tensor
    grad_norm::Tensor
    update_op::Tensor
end

"""
Create placeholders for DQN training: s, a, sp, r, done
The shape is inferred from the environment
"""
function build_placeholders(env::AbstractEnvironment)
    obs_dim = obs_dimensions(env)
    n_outs = n_actions(env)
    s = placeholder(Float32, shape=[-1, obs_dim...])
    a = placeholder(Int32, shape=[-1])
    sp = placeholder(Float32, shape=[-1, obs_dim...])
    r = placeholder(Float32, shape=[-1])
    done_mask = placeholder(Bool, shape=[-1])
    w = placeholder(Float32, shape=[-1])
    q_lo = placeholder(Float32, shape=[-1, n_actions(env)])
    q_lo_p = placeholder(Float32, shape=[-1, n_actions(env)])
    return s, a, sp, r, done_mask, w, q_lo, q_lo_p
end

function additive_correction(solver::DeepQLearningSolver, q_lo::Q, q_corr::Q, correction_weight::Float64) where Q <:Union{Array{Float64}, Tensor}
    # q_lo is already normalized ? @ outdated comment @
    if solver.weighted_correction_sum == true
        return correction_weight * q_corr + (1 - correction_weight) * q_lo
    else
        return q_corr + q_lo
    end
end

"""
Build the loss operation
relies on the Bellman equation
"""
function build_loss(solver::DeepQLearningSolver,
                    env::AbstractEnvironment,
                    q_lo::Tensor,
                    q_lo_p::Tensor,
                    q::Tensor,
                    target_q::Tensor,
                    a::Tensor,
                    r::Tensor,
                    done_mask::Tensor,
                    importance_weights::Tensor)
    if env.nn_policy == true
        env.correction_weight = 1.0
        loss, td_errors = nothing, nothing
        variable_scope("loss") do
            term = cast(done_mask, Float32)
            A = one_hot(a, n_actions(env))
            q_corr_sa = sum(A .* q, 2)
            q_sa = q_corr_sa
            q_future = target_q
            q_samp = r + (1 - term) .* discount(env.problem) .* maximum(q_future, 2)
            td_errors = q_sa - q_samp
            errors = huber_loss(td_errors)
            loss = mean(importance_weights .* errors)
        end
    else
        loss, td_errors = nothing, nothing
        variable_scope("loss") do
            term = cast(done_mask, Float32)
            A = one_hot(a, n_actions(env))
            q_corr_sa = sum(A .* q, 2)
            q_lo_sa = sum(A .* q_lo, 2)
            q_sa = additive_correction(solver, q_lo_sa, q_corr_sa, env.correction_weight)
            q_future = additive_correction(solver, q_lo_p, target_q, env.correction_weight)
            q_samp = r + (1 - term) .* discount(env.problem) .* maximum(q_future, 2)
            td_errors = q_sa - q_samp
            errors = huber_loss(td_errors)
            loss = mean(importance_weights .* errors)
        end
    end
    return loss, td_errors
end

"""
Build the loss operation with double_q
relies on the Bellman equation
"""
function build_doubleq_loss(solver::DeepQLearningSolver,
                            env::AbstractEnvironment,
                            q_lo::Tensor,
                            q_lo_p::Tensor,
                            q::Tensor,
                            target_q::Tensor,
                            qp::Tensor,
                            a::Tensor,
                            r::Tensor,
                            done_mask::Tensor,
                            importance_weights::Tensor)
    if env.nn_policy == true
        env.correction_weight = 1.0
        loss, td_errors = nothing, nothing
        variable_scope("loss") do
            term = cast(done_mask, Float32)
            A = one_hot(a, n_actions(env))
            q_corr_sa = sum(A .* q, 2)
            q_sa = q_corr_sa
            best_a = indmax(qp, 2)
            best_A = one_hot(best_a, n_actions(env))
            q_future = target_q
            target_q_best = sum(best_A.*q_future, 2)
            q_samp = r + (1 - term) .* discount(env.problem) .* target_q_best
            td_errors = q_sa - q_samp
            errors = huber_loss(td_errors)
            loss = mean(importance_weights.*errors)
        end
    else
        loss, td_errors = nothing, nothing
        variable_scope("loss") do
            term = cast(done_mask, Float32)
            A = one_hot(a, n_actions(env))
            q_corr_sa = sum(A .* q, 2)
            q_lo_sa = sum(A .* q_lo, 2)
            q_sa = additive_correction(solver, q_lo_sa, q_corr_sa, env.correction_weight)
            best_a = indmax(qp, 2)
            best_A = one_hot(best_a, n_actions(env))
            q_future = additive_correction(solver, q_lo_p, target_q, env.correction_weight)
            target_q_best = sum(best_A.*q_future, 2)
            q_samp = r + (1 - term) .* discount(env.problem) .* target_q_best
            td_errors = q_sa - q_samp
            errors = huber_loss(td_errors)
            loss = mean(importance_weights.*errors)
        end
    end
    return loss, td_errors
end

"""
Build train operation
Support gradient clipping
"""
function build_train_op(loss::Tensor;
                        lr::Union{Float64, Tensor} = 0.1,
                        grad_clip::Bool = true,
                        clip_val::Float64 = 10.,
                        optimizer_type=train.AdamOptimizer)

    optimizer = optimizer_type(lr)

    train_var = get_train_vars_by_name(Q_SCOPE)

    grad_vars = train.compute_gradients(optimizer, loss, train_var)
    clip_grads = grad_vars
    if grad_clip
        clip_grads = [(clip_by_norm(gradvar[1], clip_val), gradvar[2]) for gradvar in grad_vars]
    end
    train_op = train.apply_gradients(optimizer, clip_grads)
    grad_norm = global_norm([g[1] for g in clip_grads])
    return train_op, grad_norm
end

"""
returns a tensorflow operation to update a target network
if you run the operation, it will copy the value of the weights and biases in q_scope to
the weights and biases in target_q_scope
"""
function build_update_target_op(q_scope=Q_SCOPE, target_q_scope=TARGET_Q_SCOPE)
    q_weights = get_train_vars_by_name(q_scope)
    target_q_weights = get_train_vars_by_name(target_q_scope)

    all_ops = [tf.assign(target_q_weights[i], q_weights[i]) for i in 1:length(q_weights)]
    return update_target_op = tf.group(all_ops..., name="update_target_op")
end

function build_graph(solver::DeepQLearningSolver, env::AbstractEnvironment, graph=Graph())
    sess = init_session(graph)
    s, a, sp, r, done_mask, importance_weights, q_lo, q_lo_p = build_placeholders(env)

    q = build_q(s, solver.arch, env, scope=Q_SCOPE, dueling=solver.dueling) # not-normalized
    # q = tf.nn.l2_normalize(q, 2)
    qp = build_q(sp, solver.arch, env, scope=Q_SCOPE, reuse=true, dueling=solver.dueling)
    target_q = build_q(sp, solver.arch, env, scope=TARGET_Q_SCOPE, dueling=solver.dueling) # not-normalized

    if solver.double_q
        loss, td_errors = build_doubleq_loss(solver, env, q_lo, q_lo_p, q, target_q, qp, a, r, done_mask,  importance_weights)
    else
        loss, td_errors = build_loss(solver, env, q_lo, q_lo_p, q, target_q, a, r, done_mask, importance_weights)
    end

    train_op, grad_norm = build_train_op(loss, lr=solver.lr, grad_clip=solver.grad_clip, clip_val=solver.clip_val)

    update_op = build_update_target_op(Q_SCOPE, TARGET_Q_SCOPE)

    return TrainGraph(sess, s, a, sp, r, done_mask, importance_weights, q_lo, q_lo_p,
                      q, qp, target_q, loss, td_errors, train_op, grad_norm, update_op)
end

function get_q_lo_batch(graph::TrainGraph, env::AbstractEnvironment, s_batch, sp_batch)
    bs = size(s_batch, 1)
    q_lo_batch = zeros(bs, n_actions(env))
    q_lo_p_batch = zeros(bs, n_actions(env))
    for i = 1 : bs
        q_lo_batch[i, :] = lowfi_values(env, reshape(s_batch, bs, :)[i, :])
        q_lo_p_batch[i, :] = lowfi_values(env, reshape(sp_batch, bs, :)[i, :])
    end
    return q_lo_batch, q_lo_p_batch # normalized q_lo
end

"""
Project specific method
"""
function get_action(graph::TrainGraph, solver::DeepQLearningSolver, env::AbstractEnvironment, o::Array{Float64})
    # If o is empty return :COC
    if all(o[[i * 3 + IND_rho for i in 0 : env.num_sections - 1]] .== 0)
        return :COC
    end

    q_low_fi = vec(lowfi_values(env, o))

    # cannot take a batch of observations
    o = reshape(o, (1, size(o)...))

    if env.nn_policy == true
        # using pure neural network policy
        q_decision = run(graph.sess, graph.q, Dict(graph.s => o) )
    else
        # using corrected VICAS or VICAS
        if env.correction == true
            q_correction = run(graph.sess, graph.q, Dict(graph.s => o) ) # Array{Float32,2}
            # reshape q_correction
            temp_arr = zeros(n_actions(env)) # Array{Float64,1}
            for i in 1:n_actions(env)
                temp_arr[i] = q_correction[1, i]
            end
            q_correction = temp_arr # non-normalized
            # q_correction = normalize_qvals(temp_arr) # Array{Float64, 1}, normalized

            # corrected qvals for decision making
            if solver.weighted_correction_sum == true
                q_decision = env.correction_weight * q_correction + (1 - env.correction_weight) * q_low_fi
            else
                q_decision = q_correction + q_low_fi
            end

            # println("q_correction = ", q_correction)
            # println("q_low_fi = ", q_low_fi)
        else # no correction
            q_decision = q_low_fi
        end
    end

    if env.stochastic_ego_policy == false
        return actions(env)[indmax(q_decision)]
    else
        # softmax policy for ego agent
        q_decision -= maximum(q_decision)
        exp_qvals = exp.(q_decision)
        cat_dist = Categorical(exp_qvals / sum(exp_qvals))
        return actions(env)[rand(cat_dist)]
    end
end

###### RECURRENT VERSION

"""
    RecurrentTrainGraph
type to store all the tensor of the DQN graph including both the Q network
and the training operations
"""
mutable struct RecurrentTrainGraph
    sess::Session
    s::Tensor
    a::Tensor
    sp::Tensor
    r::Tensor
    done_mask::Tensor
    trace_mask::Tensor
    importance_weights::Tensor
    q::Tensor
    hq_in::LSTMStateTuple
    hq_out::LSTMStateTuple
    qp::Tensor
    hqp_in::LSTMStateTuple
    hqp_out::LSTMStateTuple
    target_q::Tensor
    target_hq_in::LSTMStateTuple
    target_hq_out::LSTMStateTuple
    loss::Tensor
    td_errors::Tensor
    train_op::Tensor
    grad_norm::Tensor
    update_op::Tensor
    lstm_policy::LSTMPolicy
end

# placeholder
function build_placeholders(env::AbstractEnvironment, trace_length::Int64)
    obs_dim = obs_dimensions(env)
    n_outs = n_actions(env)
    # bs x trace_length x dim
    s = placeholder(Float32, shape=[-1, trace_length, obs_dim...])
    a = placeholder(Int32, shape=[-1, trace_length])
    sp = placeholder(Float32, shape=[-1, trace_length, obs_dim...])
    r = placeholder(Float32, shape=[-1, trace_length])
    done_mask = placeholder(Bool, shape=[-1, trace_length])
    trace_mask = placeholder(Int32, shape=[-1, trace_length])
    w = placeholder(Float32, shape=[-1])
    return s, a, sp, r, done_mask, trace_mask, w
end

function build_loss(env::AbstractEnvironment,
                    q::Tensor,
                    target_q::Tensor,
                    a::Tensor,
                    r::Tensor,
                    done_mask::Tensor,
                    trace_mask::Tensor,
                    importance_weights::Tensor)
    loss, td_errors = nothing, nothing
    variable_scope("loss") do
        # flatten time dim
        trace_length = 6
        time_mask = reshape(trace_mask, (-1))
        flat_a = reshape(a, (-1))
        flat_r = reshape(r, (-1))
        flat_done_mask = reshape(done_mask, (-1))
        term = cast(flat_done_mask, Float32)
        A = one_hot(flat_a, n_actions(env))
        q_sa = sum(A.*q, 2)
        q_samp = flat_r + (1 - term).*discount(env.problem).*maximum(target_q, 2)
        td_errors = time_mask.*(q_sa - q_samp)
        errors = huber_loss(td_errors)
        loss = sum(importance_weights.*errors)/sum(time_mask)
    end
    return loss, td_errors
end

function build_doubleq_loss(env::AbstractEnvironment,
                            q::Tensor,
                            target_q::Tensor,
                            qp::Tensor,
                            a::Tensor,
                            r::Tensor,
                            done_mask::Tensor,
                            trace_mask::Tensor,
                            importance_weights::Tensor)
    loss, td_errors = nothing, nothing
    variable_scope("loss") do
        # flatten time dim
        time_mask = reshape(trace_mask, (-1))
        flat_a = reshape(a, (-1))
        flat_r = reshape(r, (-1))
        flat_done_mask = reshape(done_mask, (-1))
        term = cast(flat_done_mask, Float32)
        A = one_hot(flat_a, n_actions(env))
        q_sa = sum(A.*q, 2)
        best_a = indmax(qp, 2)
        best_A = one_hot(best_a, n_actions(env))
        target_q_best = sum(best_A.*target_q, 2)
        q_samp = flat_r + (1 - term).*discount(env.problem).*target_q_best
        td_errors = time_mask.*(q_sa - q_samp)
        errors = huber_loss(td_errors)
        loss = sum(importance_weights.*errors)/sum(time_mask)
    end
    return loss, td_errors
end

function build_graph(solver::DeepRecurrentQLearningSolver, env::AbstractEnvironment, graph=Graph())
    sess = init_session(graph)
    s, a, sp, r, done_mask, trace_mask, w = build_placeholders(env, solver.trace_length)
    q, hq_in, hq_out = build_recurrent_q_network(s,
                                                 solver.arch.convs,
                                                 solver.arch.fc_in,
                                                 solver.arch.fc_out,
                                                 solver.arch.lstm_size,
                                                 n_actions(env),
                                                 final_activation = identity,
                                                 scope= Q_SCOPE,
                                                 dueling = solver.dueling)
    qp, hqp_in, hqp_out = build_recurrent_q_network(sp,
                                                    solver.arch.convs,
                                                    solver.arch.fc_in,
                                                    solver.arch.fc_out,
                                                    solver.arch.lstm_size,
                                                    n_actions(env),
                                                    final_activation = identity,
                                                    scope= Q_SCOPE,
                                                    reuse=true,
                                                    dueling = solver.dueling)
    target_q, target_hq_in, target_hq_out = build_recurrent_q_network(sp,
                                                                      solver.arch.convs,
                                                                      solver.arch.fc_in,
                                                                      solver.arch.fc_out,
                                                                      solver.arch.lstm_size,
                                                                      n_actions(env),
                                                                      final_activation = identity,
                                                                      scope= TARGET_Q_SCOPE,
                                                                      reuse=false,
                                                                      dueling = solver.dueling)
    if solver.double_q
        loss, td_errors = build_doubleq_loss(env, q, target_q, qp, a, r, done_mask, trace_mask, w)
    else
        loss, td_errors = build_loss(env, q, target_q, a, r, done_mask, trace_mask, w)
    end
    train_op, grad_norm = build_train_op(loss,
                                         lr=solver.lr,
                                         grad_clip=solver.grad_clip,
                                         clip_val=solver.clip_val)
    update_op = build_update_target_op(Q_SCOPE, TARGET_Q_SCOPE)
    lstm_policy = LSTMPolicy(sess, env, solver.arch, solver.dueling)
    return RecurrentTrainGraph(sess, s, a, sp, r, done_mask, trace_mask, w,
                               q, hq_in, hq_out,
                               qp, hqp_in, hqp_out,
                               target_q, target_hq_in, target_hq_out,
                               loss, td_errors, train_op, grad_norm, update_op, lstm_policy)
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
"""
    bdp(session::Session, security::AbstractString, fields;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        )

Given a single field name or vector of field names
at the `fields` argument,
return a single named tuple with the result of a
`ReferenceDataRequest` request.

See [`ErrorHandling`](@ref) for `error_handling` argument behavior.

For bulk data, `bds` method should be used instead.

# Example

```julia
julia> BLPData.bdp(session, "PETR4 BS Equity", "PX_LAST")
(PX_LAST = 15.95,)

julia> BLPData.bdp(session, "PETR4 BS Equity", ["PX_LAST", "VOLUME"])
(PX_LAST = 15.95, VOLUME = 1.601771e8)
```
"""
function bdp(session::Session, security::AbstractString, fields::Vector{T};
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T<:AbstractString}

    bdp_result = bdp(session, [security], fields, options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds, error_handling=error_handling)
    return bdp_result[security]
end

function bdp(session::Session, security::AbstractString, fields::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        )

    bdp(session, security, [fields], options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds, error_handling=error_handling)
end

function bdp(session::Session, securities::Vector{T1}, fields::Vector{T2};
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T1<:AbstractString, T2<:AbstractString}

    @assert !isempty(fields) "Fields vector should not be empty."

    queue, corr_id = send_request(session, "//blp/refdata", "ReferenceDataRequest") do req
        append!(req["securities"], securities)
        append!(req["fields"], fields)

        if options != nothing
            for (k, v) in options
                req[k] = v
            end
        end
    end

    local result::Dict{String, Any} = Dict()

    for_each_response_message_element(queue, corr_id, timeout_milliseconds=timeout_milliseconds, verbose=verbose) do element
        parse_reference_data_response_into!(element, result, securities, fields, error_handling)
    end

    if length(result) != length(securities)
        throw(BLPUnknownException("Expected to return information about $(length(securities)) securities. Got $(length(result))."))
    end

    return result
end

function bdp(session::Session, securities::Vector{T}, field::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0),
            error_handling::ErrorHandling=Unwrap()
        ) where {T<:AbstractString}

    bdp(session, securities, [field], options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds, error_handling=error_handling)
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using Test
using PETScX

@testset "options tests" begin
    kw_opts = (
        ksp_monitor = nothing,
        ksp_view = nothing,
        da_grid_x = 100,
        da_grid_y = 100,
        pc_type = "mg",
        pc_mg_levels = 1,
        mg_levels_0_pc_type = "ilu",
    )
    for petsclib in PETScX.petsclibs
        PETScX.Initialize(petsclib)

        opts = PETScX.Options(petsclib; kw_opts...)

        # Check that all the keys got added
        for (key, val) in pairs(kw_opts)
            key = string(key)
            val = isnothing(val) ? "" : string(val)
            @test val == opts[key]
        end

        # Check that we throw when a bad key is asked for
        @test_throws KeyError opts["bad key"]

        # try to insert some new keys
        opts["new_opt"] = 1
        opts["nothing_opt"] = nothing
        @test "" == opts["nothing_opt"]
        @test "1" == opts["new_opt"]

        # Check that viewer is working
        _stdout = stdout
        (rd, wr) = redirect_stdout()
        @show opts

        @test readline(rd) == "opts = #PETSc Option Table entries:"
        @test readline(rd) == "-da_grid_x 100"
        @test readline(rd) == "-da_grid_y 100"
        @test readline(rd) == "-ksp_monitor"
        @test readline(rd) == "-ksp_view"
        @test readline(rd) == "-mg_levels_0_pc_type ilu"
        @test readline(rd) == "-new_opt 1"
        @test readline(rd) == "-nothing_opt"
        @test readline(rd) == "-pc_mg_levels 1"
        @test readline(rd) == "-pc_type mg"
        @test readline(rd) == "#End of PETSc Option Table entries"
        @test readline(rd) == ""

        glo_opts = PETScX.GlobalOptions(petsclib)
        show(stdout, "text/plain", glo_opts)
        @test readline(rd) == "#No PETSc Option Table entries"

        # Try to set some options and check that they are set
        PETScX.with(opts) do
            show(stdout, "text/plain", glo_opts)
        end
        @test readline(rd) == "#PETSc Option Table entries:"
        @test readline(rd) == "-da_grid_x 100"
        @test readline(rd) == "-da_grid_y 100"
        @test readline(rd) == "-ksp_monitor"
        @test readline(rd) == "-ksp_view"
        @test readline(rd) == "-mg_levels_0_pc_type ilu"
        @test readline(rd) == "-new_opt 1"
        @test readline(rd) == "-nothing_opt"
        @test readline(rd) == "-pc_mg_levels 1"
        @test readline(rd) == "-pc_type mg"
        @test readline(rd) == "#End of PETSc Option Table entries"

        show(stdout, "text/plain", glo_opts)
        @test readline(rd) == "#No PETSc Option Table entries"

        redirect_stdout(_stdout)

        PETScX.Finalize(petsclib)
    end
end


text:
using ITensors
using ITensorInfiniteMPS
using Random

# Helper function to make an MPO
import ITensors: op
op(::OpName"Zero", ::SiteType"S=1/2", s::Index) = ITensor(s', dag(s))

function tfi_mpo(s, l, r; J=1.0, h)
  dₗ = 3 # The link dimension of the TFI
  Hmat = fill(op("Zero", s), dₗ, dₗ)
  Hmat[1, 1] = op("Id", s)
  Hmat[2, 2] = op("Id", s)
  Hmat[3, 3] = op("Id", s)
  Hmat[2, 1] = -J * op("X", s)
  Hmat[3, 2] = -J * op("X", s)
  Hmat[3, 1] = -h * op("Z", s)
  H = ITensor()
  for i in 1:dₗ, j in 1:dₗ
    H += Hmat[i, j] * setelt(l => i) * setelt(r => j)
  end
  return H
end

Random.seed!(1234)

N = 1
s = siteinds("S=1/2", N)
χ = 6
@assert iseven(χ)
if any(hasqns, s)
  space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
else
  space = χ
end
#ψ = InfiniteMPS(ComplexF64, s; space = space)
ψ = InfiniteMPS(s; space=space)
randn!.(ψ)

# Use a finite MPO to create the infinite MPO
H = InfiniteMPO(N)

# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
J = 1.0
h = 2.0

# For the finite unit cell
s¹ = addtags.(s, "c=1")
l¹ = [Index(3, "Link,l=$n,c=1") for n in 1:N]
for n in 1:N
  # TODO: use a CelledVector here to make the code logic simpler
  l¹ₙ = n == 1 ? replacetags(dag(l¹[N]), "c=1" => "c=0") : l¹[n - 1]
  r¹ₙ = l¹[n]
  H[n] = tfi_mpo(s¹[n], l¹ₙ, r¹ₙ; J=J, h=h)
end

ψf = vumps(H, ψ; nsweeps=10)

s⃗ = [uniqueind(ψ[n], ψ[n - 1], ψ[n + 1]) for n in 1:2]
h = -J * op("X", s⃗, 1) * op("X", s⃗, 2) - h * op("Z", s⃗, 1) * op("Id", s⃗, 2)
#ψ2 = ψf.AL[1] * ψf.C[1] * ψf.AR[2]
ψ2 = ψf.AL[1] * ψf.AL[2] * ψf.C[2]
@show ψ2 * h * prime(dag(ψ2), "Site")


text:
const MPIInteger = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64}
const MPIFloatingPoint = Union{Float32, Float64}
const MPIComplex = Union{ComplexF32, ComplexF64}

const MPIDatatype = Union{Char,
                            Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64,
                            UInt64,
                            Float32, Float64, ComplexF32, ComplexF64}
MPIBuffertype{T} = Union{Ptr{T}, Array{T}, SubArray{T}, Ref{T}}

MPIBuffertypeOrConst{T} = Union{MPIBuffertype{T}, SentinelPtr}

Base.cconvert(::Type{MPIPtr}, x::Union{Ptr{T}, Array{T}, Ref{T}}) where T = Base.cconvert(Ptr{T}, x)
function Base.cconvert(::Type{MPIPtr}, x::SubArray{T}) where T
    Base.cconvert(Ptr{T}, x)
end
function Base.unsafe_convert(::Type{MPIPtr}, x::MPIBuffertype{T}) where T
    ptr = Base.unsafe_convert(Ptr{T}, x)
    reinterpret(MPIPtr, ptr)
end

function Base.cconvert(::Type{MPIPtr}, x::String)
    x
end
function Base.unsafe_convert(::Type{MPIPtr}, x::String)
    reinterpret(MPIPtr, pointer(x))
end

function Base.cconvert(::Type{MPIPtr}, ::Nothing)
    reinterpret(MPIPtr, C_NULL)
end

macro assert_minlength(buffer, count)
    quote
        if $(esc(buffer)) isa AbstractArray
            @assert length($(esc(buffer))) >= $(esc(count))
        end
    end
end

"""
    MPI.MPIPtr

A pointer to an MPI buffer. This type is used only as part of the implicit conversion in
`ccall`: a Julia object can be passed to MPI by defining methods for
`Base.cconvert(::Type{MPIPtr}, ...)`/`Base.unsafe_convert(::Type{MPIPtr}, ...)`.

Currently supported are:
 - `Ptr`
 - `Ref`
 - `Array`
 - `SubArray`
 - `CUDA.CuArray` if CUDA.jl is loaded.

Additionally, certain sentinel values can be used, e.g. `MPI_IN_PLACE` or `MPI_BOTTOM`.
"""
MPIPtr

# MPI_IN_PLACE

struct InPlace
end
Base.cconvert(::Type{MPIPtr}, ::InPlace) = MPI_IN_PLACE

"""
    MPI.IN_PLACE

A sentinel value that can be passed as a buffer argument for certain collective operations
to use the same buffer for send and receive operations.

- [`Scatter!`](@ref) and [`Scatterv!`](@ref): can be used as the `recvbuf` argument on the
  root process.

- [`Gather!`](@ref) and [`Gatherv!`](@ref): can be used as the `sendbuf` argument on the
  root process.

- [`Allgather!`](@ref), [`Allgatherv!`](@ref), [`Alltoall!`](@ref) and
  [`Alltoallv!`](@ref): can be used as the `sendbuf` argument on all processes.

- [`Reduce!`](@ref) (root only), [`Allreduce!`](@ref), [`Scan!`](@ref) and
  [`Exscan!`](@ref): can be used as `sendbuf` argument.

"""
const IN_PLACE = InPlace()

# TODO: MPI_BOTTOM

"""
    MPI.Buffer

An MPI buffer for communication with a single rank. It is used for point-to-point
communication and some collective operations.

# Fields
$(DocStringExtensions.FIELDS)

# Usage

    Buffer(data, count::Integer, datatype::Datatype)

Generic constructor.

    Buffer(data)

Construct a `Buffer` backed by `data`, automatically determining the appropriate `count`
and `datatype`. Methods are provided for

 - `Ref`
 - `Array`
 - `CUDA.CuArray` if CUDA.jl is loaded
 - `SubArray`s of an `Array` or `CUDA.CuArray` where the layout is contiguous, sequential or
   blocked.

# See also

- [`Buffer_send`](@ref)

"""
struct Buffer{A}
    """a Julia object referencing a region of memory to be used for communication. It is
    required that the object can be `cconvert`ed to an [`MPIPtr`](@ref)."""
    data::A

    """the number of elements of `datatype` in the buffer. Note that this may not
    correspond to the number of elements in the array if derived types are used."""
    count::Cint

    """the [`MPI.Datatype`](@ref) stored in the buffer."""
    datatype::Datatype
end
Buffer(buf::Buffer) = buf
Buffer(data, count::Integer, datatype::Datatype) = Buffer(data, Cint(count), datatype)

function Buffer(arr::Array)
    Buffer(arr, Cint(length(arr)), Datatype(eltype(arr)))
end
function Buffer(ref::Ref)
    Buffer(ref, Cint(1), Datatype(eltype(ref)))
end

# SubArray
function Buffer(sub::Base.FastContiguousSubArray)
    Buffer(sub, Cint(length(sub)), Datatype(eltype(sub)))
end
function Buffer(sub::Base.FastSubArray)
    datatype = Types.create_vector(length(sub), 1, sub.stride1,
                                   Datatype(eltype(sub); commit=false))
    Types.commit!(datatype)
    Buffer(sub, Cint(1), datatype)
end
function Buffer(sub::SubArray{T,N,P,I,false}) where {T,N,P,I<:Tuple{Vararg{Union{Base.ScalarIndex, Base.Slice, AbstractUnitRange}}}}
    datatype = Types.create_subarray(size(parent(sub)),
                                     map(length, sub.indices),
                                     map(i -> first(i)-1, sub.indices),
                                     Datatype(eltype(sub), commit=false))
    Types.commit!(datatype)
    Buffer(parent(sub), Cint(1), datatype)
end

Buffer(::InPlace) = Buffer(IN_PLACE, 0, DATATYPE_NULL)
Buffer(::Nothing) = Buffer(nothing, 0, DATATYPE_NULL)

"""
    Buffer_send(data)

Construct a [`Buffer`](@ref) object for a send operation from `data`, allowing cases where
`isbits(data)`.
"""
Buffer_send(data) = isbits(data) ? Buffer(Ref(data)) : Buffer(data)
Buffer_send(str::String) = Buffer(str, sizeof(str), MPI.CHAR)
Buffer_send(::InPlace) = Buffer(InPlace())
Buffer_send(::Nothing) = Buffer(nothing)

"""
    MPI.UBuffer

An MPI buffer for chunked collective communication, where all chunks are of uniform size.

# Fields
$(DocStringExtensions.FIELDS)

# Usage

    UBuffer(data, count::Integer, nchunks::Union{Nothing, Integer}, datatype::Datatype)

Generic constructor.

    UBuffer(data, count::Integer)

Construct a `UBuffer` backed by `data`, where `count` is the number of elements in each chunk.

# See also

- [`VBuffer`](@ref): similar, but supports chunks of non-uniform sizes.
"""
struct UBuffer{A}
    """A Julia object referencing a region of memory to be used for communication. It is
    required that the object can be `cconvert`ed to an [`MPIPtr`](@ref)."""
    data::A

    """The number of elements of `datatype` in each chunk."""
    count::Cint

    """The maximum number of chunks stored in the buffer. This is used only for
    validation, and can be set to `nothing` to disable checks."""
    nchunks::Union{Nothing,Cint}

    """The [`MPI.Datatype`](@ref) stored in the buffer."""
    datatype::Datatype
end
UBuffer(data, count::Integer, nchunks::Union{Integer, Nothing}, datatype::Datatype) =
    UBuffer(data, Cint(count), nchunks isa Integer ? Cint(nchunks) : nothing, datatype)

function UBuffer(arr::AbstractArray, count::Integer)
    @assert stride(arr, 1) == 1
    UBuffer(arr, count, div(length(arr), count), Datatype(eltype(arr)))
end
Base.similar(buf::UBuffer) =
    UBuffer(similar(buf.data), buf.count, buf.nchunks, buf.datatype)

UBuffer(::Nothing) = UBuffer(nothing, 0, nothing, DATATYPE_NULL)
UBuffer(::InPlace) = UBuffer(IN_PLACE, 0, nothing, DATATYPE_NULL)

"""
    MPI.VBuffer

An MPI buffer for chunked collective communication, where chunks can be of different sizes and at different offsets.

# Fields
$(DocStringExtensions.FIELDS)

# Usage

    VBuffer(data, counts[, displs[, datatype]])

Construct a `VBuffer` backed by `data`, where `counts[j]` is the number of elements in the
`j`th chunk, and `displs[j]` is the 0-based displacement. In other words, the `j`th chunk
occurs in indices `displs[j]+1:displs[j]+counts[j]`.

The default value for `displs[j] = sum(counts[1:j-1])`.

# See also

- [`UBuffer`](@ref) when chunks are all of the same size.
"""
struct VBuffer{A}
    """A Julia object referencing a region of memory to be used for communication. It is
    required that the object can be `cconvert`ed to an [`MPIPtr`](@ref)."""
    data::A

    """An array containing the length of each chunk."""
    counts::Vector{Cint}

    """An array containing the (0-based) displacements of each chunk."""
    displs::Vector{Cint}

    """The [`MPI.Datatype`](@ref) stored in the buffer."""
    datatype::Datatype
end
VBuffer(data, counts, displs, datatype::Datatype) =
    VBuffer(data, convert(Vector{Cint}, counts),
                 convert(Vector{Cint}, displs), datatype)
VBuffer(data, counts, displs) =
    VBuffer(data, counts, displs, Datatype(eltype(data)))

function VBuffer(arr::AbstractArray, counts)
    @assert stride(arr,1) == 1
    counts = convert(Vector{Cint}, counts)
    displs = similar(counts)
    d = zero(Cint)
    for i in eachindex(displs)
        displs[i] = d
        d += counts[i]
    end
    @assert length(arr) >= d
    VBuffer(arr, counts, displs, Datatype(eltype(arr)))
end

VBuffer(::Nothing) = VBuffer(nothing, Cint[], Cint[], DATATYPE_NULL)
VBuffer(::InPlace) = VBuffer(IN_PLACE, Cint[], Cint[], DATATYPE_NULL)

"""
    MPI.RBuffer

An MPI buffer for reduction operations ([`MPI.Reduce!`](@ref), [`MPI.Allreduce!`](@ref), [`MPI.Scan!`](@ref), [`MPI.Exscan!`](@ref)).

# Fields
$(DocStringExtensions.FIELDS)

# Usage

    RBuffer(senddata, recvdata[, count, datatype])

Generic constructor.

    RBuffer(senddata, recvdata)

Construct a `Buffer` backed by `senddata` and `recvdata`, automatically determining the
appropriate `count` and `datatype`.

- `senddata` can be [`MPI.IN_PLACE`](@ref)
- `recvdata` can be `nothing` on a non-root node with [`MPI.Reduce!`](@ref)
"""
struct RBuffer{S,R}
    """A Julia object referencing a region of memory to be used for the send buffer. It is
    required that the object can be `cconvert`ed to an [`MPIPtr`](@ref)."""
    senddata::S

    """A Julia object referencing a region of memory to be used for the receive buffer. It is
    required that the object can be `cconvert`ed to an [`MPIPtr`](@ref)."""
    recvdata::R

    """the number of elements of `datatype` in the buffer. Note that this may not
    correspond to the number of elements in the array if derived types are used."""
    count::Cint

    """the [`MPI.Datatype`](@ref) stored in the buffer."""
    datatype::Datatype
end

RBuffer(senddata, recvdata, count::Integer, datatype::Datatype) =
    RBuffer(senddata, recvdata, Cint(count), datatype)

function RBuffer(senddata::AbstractArray{T}, recvdata::AbstractArray{T}) where {T}
    @assert (count = length(senddata)) == length(recvdata)
    @assert stride(senddata,1) == stride(recvdata,1) == 1
    RBuffer(senddata, recvdata, count, Datatype(T))
end
function RBuffer(::InPlace, recvdata::AbstractArray{T}) where {T}
    count = length(recvdata)
    @assert stride(recvdata,1) == 1
    RBuffer(IN_PLACE, recvdata, count, Datatype(T))
end
function RBuffer(senddata::AbstractArray{T}, recvdata::Nothing) where {T}
    count = length(senddata)
    @assert stride(senddata,1) == 1
    RBuffer(senddata, nothing, count, Datatype(T))
end

function RBuffer(senddata::Ref{T}, recvdata::Ref{T}) where {T}
    RBuffer(senddata, recvdata, 1, Datatype(T))
end
function RBuffer(senddata::InPlace, recvdata::Ref{T}) where {T}
    RBuffer(IN_PLACE, recvdata, 1, Datatype(T))
end
function RBuffer(senddata::Ref{T}, recvdata::Nothing) where {T}
    RBuffer(senddata, nothing, 1, Datatype(T))
end

Base.eltype(rbuf::RBuffer) = eltype(rbuf.senddata)
Base.eltype(rbuf::RBuffer{InPlace}) = eltype(rbuf.recvdata)


text:
abstract type AbstractBSONCodec end

struct MinimalCodec <: AbstractBSONCodec end

struct BSONCodec <: AbstractBSONCodec
    dicttype::Type
    as_symbol::Bool
    as_string::Bool
end

function BSONCodec(dicttype = Dict; as_symbol = true, as_string = true)
    BSONCodec(dicttype, as_symbol, as_string)
end

struct JSONCodec <: AbstractBSONCodec
    dicttype::Type
    as_symbol::Bool
    as_string::Bool
end

function JSONCodec(dicttype = Dict; as_symbol = true, as_string = true)
    JSONCodec(dicttype, as_symbol, as_string)
end

const _codec = Dict(
                    :minimal => MinimalCodec(),
                    :bson => BSONCodec(),
                    :json => JSONCodec()
                   )

struct Element
    key::String
    value::Any

    function Element(k::AbstractString, v)
        ('\0' in k) && throw(ArgumentError("Element key cannot have embbeded null"))
        new(k, v)
    end
end

struct Maxkey end
struct Minkey end

const maxkey = Maxkey()
const minkey = Minkey()

struct ObjectId
    value::String

    function ObjectId(bytes::Vector{UInt8})
        length(bytes) == 12 || throw(ArgumentError("ObjectId buffer shall be 12 bytes"))
        new(String(bytes))
    end
end

function _generate_oid_time(buf::Vector{UInt8})
    io = IOBuffer()
    write(io, hton(UInt32(div(value(now(UTC)) - UNIXEPOCH, 1000))))
    append!(buf, take!(io))
end

function _generate_oid_counter(buf::Vector{UInt8})
    global _oid_context

    counter = UInt32(_oid_context[:counter] + 1)
    _oid_context[:counter] = counter

    io = IOBuffer()
    write(io, hton(counter))
    append!(buf, take!(io)[2:end])
end

function ObjectId()

    ret = Vector{UInt8}()
    sizehint!(ret, 12)
    _generate_oid_time(ret)
    append!(ret, _oid_context[:rand])
    _generate_oid_counter(ret)

    ObjectId(ret)
end

ObjectId(str::AbstractString) = ObjectId(hex2bytes(str))

string(oid::ObjectId) = bytes2hex(Vector{UInt8}(oid))

show(io::IO, oid::ObjectId) = write(io, "ObjectId($(string(oid)))")

Vector{UInt8}(oid::ObjectId) = Vector{UInt8}(oid.value)

struct Binary
    subtype::UInt8
    bytes::Vector{UInt8}

    function Binary(t::UInt8, b::Vector{UInt8})
        t <= 0x05 || t >= 0x80 || throw(ArgumentError("Invalid Binary subtype $t"))
        new(t, b)
    end
end

struct BSONSymbol
    value::String
end

struct Timestamp
    value::UInt64
end

struct Decimal128
    value::UInt128
end

struct BSONRegex
    pattern::String
    options::String

    function BSONRegex(p::AbstractString, o::AbstractString)
        ('\0' in p) && throw(ArgumentError("BSON regex pattern cannot have embbeded null"))
        ('\0' in o) && throw(ArgumentError("BSON regex options cannot have embbeded null"))
        new(p, o)
    end
end

struct DBPointer
    ref::String
    id::ObjectId
end

struct Document
    elist::Vector{Element}
end

function getindex(doc::Document, key::AbstractString)
    for kv in doc.elist
        if kv.key == key
            return kv.value
        end
    end
    throw(KeyError("$key not found"))
end

struct BSONArray
    elist::Vector{Element}
end

getindex(doc::Document, key::Integer) = doc.elist[key].value

struct Code
    code::String
end

struct CodeWithScope
    code::String
    scope::Document
end

const bson_type = Dict(
                       0x01 => Float64,
                       0x02 => String,
                       0x03 => Document,
                       0x04 => BSONArray,
                       0x05 => Binary,
                       0x06 => Missing, # undefined
                       0x07 => ObjectId,
                       0x08 => Bool,
                       0x09 => DateTime,
                       0x0A => Nothing, # null
                       0x0B => BSONRegex,
                       0x0C => DBPointer,
                       0x0D => Code,
                       0x0E => BSONSymbol,
                       0x0F => CodeWithScope,
                       0x10 => Int32,
                       0x11 => Timestamp,
                       0x12 => Int64,
                       0x13 => Decimal128,
                       0x7F => Maxkey,
                       0xFF => Minkey,
                      )

for (k, v) in bson_type
    @eval bson_code(::$v) = $k
end

const BSONType = Union{
                       Float64,
                       String,
                       Document,
                       BSONArray,
                       Binary,
                       Missing, # undefined
                       ObjectId,
                       Bool,
                       DateTime,
                       Nothing, # null
                       BSONRegex,
                       DBPointer,
                       Code,
                       BSONSymbol,
                       CodeWithScope,
                       Int32,
                       Timestamp,
                       Int64,
                       Decimal128,
                       Maxkey,
                       Minkey,
                      }


text:
using Test

@testset "doctest" begin
    include("../docs/utils.jl")
    transducers_rm_examples()
    transducers_makedocs(examples=[])
    @test true
end


text:
module LoadModule

# Currently, it is possible for community code to access everything on the system and thus it can easally steal the keys. Ideally it would be cool if the modules could be loaded with some sort of restricted namespace which could access only given folders on the system. But that seems unlikelly to be a focus of julia developers. Instead one could try to start a julia process in the container which does not have a filesystem access at all (except read only access to the community data) and then connect it over a socket and make a wrapper for restriction. Also worth exploring Cassete.jl

### Could be reduced in length
# function loadmodule(m::Symbol)
#     expr = :(import $m)
#     @eval $expr
#     m = @eval $m
#     return m
# end

import Pkg.Types.UUID
import Pkg.Types.Context

function loadmodule(ctx::Context,uuid::UUID)
    name = Symbol(ctx.env.manifest[uuid].name)
    expr = :(import $name)
    @eval $expr
    #return @eval $name
end

end

import Pkg.Types.UUID
importcommunity(uuid::UUID) = LoadModule.loadmodule(Context(),uuid)

addcommunity(spec::PackageSpec) = addcommunity(spec::PackageSpec,CONFIG_DIR)
addcommunity(name) = addcommunity(PackageSpec(name))
addcommunity(;url=nothing) = addcommunity(PackageSpec(url=url))

### Would be usefull for setting up the system for the developer
function uuid(name::AbstractString)
    ctx = Context()
    return ctx.env.project.deps[name]
end

function adduuid!(spec::PackageSpec)
    name = basename(spec.repo.url)[1:end-3]
    spec.uuid = uuid(name)
end

# The question however how does

function addcommunity(spec::PackageSpec,CONFIG_DIR)
    mkpath(CONFIG_DIR * "/communities/")
    withenv("JULIA_PKG_DEVDIR" => CONFIG_DIR * "/communities/") do     # Good enough for the prototype
        Pkg.develop(spec)
    end
end


text:
# this example is transported from Andrew Delong and Anton Osokin's Matlab wrapper. All credit goes to them.
using GCoptimization

h = gco_create(4,3)     # % Create new object with NumSites=4, NumLabels=3

gco_setdatacost(h, [
    0 9 2 0;            # % Sites 1,4 prefer  label 1
    3 0 3 3;            # % Site  2   prefers label 2 (strongly)
    5 9 0 5])           # % Site  3   prefers label 3

gco_setsmoothcost(h, [
    0 1 2;              # %
    1 0 1;              # % Linear (Total Variation) pairwise cost
    2 1 0])             # %

gco_setneighbors(h, [
    0 1 0 0;            # % Sites 1 and 2 connected with weight 1
    0 0 1 0;            # % Sites 2 and 3 connected with weight 1
    0 0 0 2;            # % Sites 3 and 4 connected with weight 2
    0 0 0 0])

gco_expansion(h)        # % Compute optimal labeling via alpha-expansion

gco_getlabeling(h)      # % Optimal labeling is (1,2,1,1)

gco_energy(h)           # % Energy = Data Energy + Smooth Energy + Label Energy


text:
function TriMesh_ProjectDerivatives(x,eConn,u;nodeList=[])
#= """
  TriMesh_ProjectDerivatives - project the derivative of a C0 scalar field
             onto the continuous finite element space (a ZZ projection).

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    dudx1, dudx2, elError, nodeList =
           TriMesh_ProjectDerivatives(x, eConn, u, nodeList)
  ```
                 Multiple scalar fields can be treated at once
                     u = [u1  u2  ... ]

                  in which case, the function returns
                     dudx1 = [∂(s1)/∂(x1)  ∂(s2)/∂(x1)  ... ]
                     dudx2 = [∂(s1)/∂(x1)  ∂(s2)/∂(x1)  ... ]

  Arguments:
  - `x`: Nodal coordinates
  - `eConn`: Element connectivity
  - `u`: Nodal values of scalar quantity

  - `nodeList`: (optional, can be omitted)
            the nodeList structure can be reused for multiple
            projections involving the same mesh, e.g. for
            time-dependent problems. It is generated and
            returned if not provided.

  Output arguments:
  - `dudx1`: Projection of the x-derivatives evaluated at nodes
  - `dudx2`: Projection of the y-derivatives evaluated at nodes

  - `elError`: element error (H1-seminorm for each field)
  - `nodeList`: a structure that contains an element list for each node
""" =#

  nNodes = size(x,1)
  nElements = size(eConn,1)

  nFields = size(u,2)

  #----------------------------------------------------------------------------
  #  For every vertex node, construct a list of elements that share it
  #----------------------------------------------------------------------------
  if typeof(nodeList) == typeof([])
    nodeList = zeros(UInt64,nNodes,15) # the first entry is reserved as a
                                       # counter and
                                       # we assume <= 14 elements per vertex

    for n=1:nElements  # local nodes 1,2, and 3 are vertices
      nodeList[eConn[n,1],1] = nodeList[eConn[n,1],1] + 1
      nodeList[eConn[n,1],nodeList[eConn[n,1],1]+1] = n

      nodeList[eConn[n,2],1] = nodeList[eConn[n,2],1] + 1
      nodeList[eConn[n,2],nodeList[eConn[n,2],1]+1] = n

      nodeList[eConn[n,3],1] = nodeList[eConn[n,3],1] + 1
      nodeList[eConn[n,3],nodeList[eConn[n,3],1]+1] = n
    end
  end

  #----------------------------------------------------------------------------
  #  Compute the least-squares projection over each patch,
  #     least-squares of gradients are performed at element quadrature points
  #----------------------------------------------------------------------------
  nQuadrature = 7
  r,s,w = twodQuadratureRule(nQuadrature)

  d1_p = fill(Inf,nNodes,nFields)
  d2_p = fill(Inf,nNodes,nFields)  # initialize to ridiculous values

  for n=1:nNodes
    if ( nodeList[n,1]>0 ) # this is a vertex node
      nel = nodeList[n,1]
      P   = zeros(nel*nQuadrature,6)
      d1  = zeros(nel*nQuadrature,nFields)
      d2  = zeros(nel*nQuadrature,nFields)

      for j=1:nodeList[n,1] # loop over elements in the patch
        el = nodeList[n,j+1]
        nodesLocal = eConn[el,:]
        xLocal = x[nodesLocal,:]
        uLocal = u[nodesLocal,:]

        x_g, w_g, ϕ, ϕ_x, ϕ_y = twodShape(xLocal, r, s, w)

        d_x = ϕ_x*uLocal   # calculate the finite element derivatives
        d_y = ϕ_y*uLocal   # on this element

        idx = 1+(j-1)*nQuadrature:j*nQuadrature
        P[idx,:] = hcat( ones(nQuadrature,1), x_g[:,1], x_g[:,2], x_g[:,1].^2, x_g[:,1].*x_g[:,2], x_g[:,2].^2 )
        d1[idx,:] = d_x
        d2[idx,:] = d_y

#       for i=1:n_quadrature
#         P  = [ P ; ...
#                1 x_g(i,1) x_g(i,2) x_g(i,1)^2 x_g(i,1)*x_g(i,2) x_g(i,2)^2 ]
#         d1 = [ d1; d_x(i) ]
#         d2 = [ d2; d_y(i) ]
#       end
      end # finished filling least-squares system

      # calculate the polynomial coefficients for each field
      for i=1:nFields
        a_x = P\d1[:,i]
        a_y = P\d2[:,i]

        #-------------------------------------------------------------------78--
        #  Compute projected derivatives at vertex node
        #-----------------------------------------------------------------------
        xp  = x[n,1]
        yp  = x[n,2]
        d1_p[n,i] = a_x[1]      + a_x[2]*xp    + a_x[3]*yp +
                    a_x[4]*xp^2 + a_x[5]*xp*yp + a_x[6]*yp^2
        d2_p[n,i] = a_y[1]      + a_y[2]*xp    + a_y[3]*yp +
                    a_y[4]*xp^2 + a_y[5]*xp*yp + a_y[6]*yp^2

        for j=1:nodeList[n,1]
          el = nodeList[n,j+1]
          nodesLocal = eConn[el,:]
          xLocal = x[nodesLocal,:]

          #---------------------------------------------------------------------
          #  Compute contribution to edge nodes.  These are either on the
          #  boundary or are averaged over two element faces.
          #---------------------------------------------------------------------
          if ( n==nodesLocal[1] )
            x4 = xLocal[4,1]
            y4 = xLocal[4,2]

            x6 = xLocal[6,1]
            y6 = xLocal[6,2]

            if ( d1_p[nodesLocal[4],i] == Inf ) # contribution not rec'd
              d1_p[nodesLocal[4],i] =
                  a_x[1]      + a_x[2]*x4    + a_x[3]*y4 +
                  a_x[4]*x4^2 + a_x[5]*x4*y4 + a_x[6]*y4^2
              d2_p[nodesLocal[4],i] =
                  a_y[1]      + a_y[2]*x4    + a_y[3]*y4 +
                  a_y[4]*x4^2 + a_y[5]*x4*y4 + a_y[6]*y4^2
            else
              d1_p[nodesLocal[4],i] = .5*( d1_p[nodesLocal[4],i] +
                  a_x[1]      + a_x[2]*x4    + a_x[3]*y4 +
                  a_x[4]*x4^2 + a_x[5]*x4*y4 + a_x[6]*y4^2 )
              d2_p[nodesLocal[4],i] = .5*( d2_p[nodesLocal[4],i] +
                  a_y[1]      + a_y[2]*x4    + a_y[3]*y4 +
                  a_y[4]*x4^2 + a_y[5]*x4*y4 + a_y[6]*y4^2 )
            end
            if ( d1_p[nodesLocal[6],i] == Inf )
              d1_p[nodesLocal[6],i] =
                  a_x[1]      + a_x[2]*x6    + a_x[3]*y6 +
                  a_x[4]*x6^2 + a_x[5]*x6*y6 + a_x[6]*y6^2
              d2_p[nodesLocal[6],i] =
                  a_y[1]      + a_y[2]*x6    + a_y[3]*y6 +
                  a_y[4]*x6^2 + a_y[5]*x6*y6 + a_y[6]*y6^2
            else
              d1_p[nodesLocal[6],i] = .5*( d1_p[nodesLocal[6],i] +
                  a_x[1]      + a_x[2]*x6    + a_x[3]*y6 +
                  a_x[4]*x6^2 + a_x[5]*x6*y6 + a_x[6]*y6^2 )
              d2_p[nodesLocal[6],i] = .5*( d2_p[nodesLocal[6],i] +
                  a_y[1]      + a_y[2]*x6    + a_y[3]*y6 +
                  a_y[4]*x6^2 + a_y[5]*x6*y6 + a_y[6]*y6^2 )
            end

          elseif ( n==nodesLocal[2] )
            x4 = xLocal[4,1]
            y4 = xLocal[4,2]

            x5 = xLocal[5,1]
            y5 = xLocal[5,2]

            if ( d1_p[nodesLocal[4],i] == Inf )
              d1_p[nodesLocal[4],i] =
                  a_x[1]      + a_x[2]*x4    + a_x[3]*y4 +
                  a_x[4]*x4^2 + a_x[5]*x4*y4 + a_x[6]*y4^2
              d2_p[nodesLocal[4],i] =
                  a_y[1]      + a_y[2]*x4    + a_y[3]*y4 +
                  a_y[4]*x4^2 + a_y[5]*x4*y4 + a_y[6]*y4^2
            else
              d1_p[nodesLocal[4],i] = .5*( d1_p[nodesLocal[4],i] +
                  a_x[1]      + a_x[2]*x4    + a_x[3]*y4 +
                  a_x[4]*x4^2 + a_x[5]*x4*y4 + a_x[6]*y4^2 )
              d2_p[nodesLocal[4],i] = .5*( d2_p[nodesLocal[4],i] +
                  a_y[1]      + a_y[2]*x4    + a_y[3]*y4 +
                  a_y[4]*x4^2 + a_y[5]*x4*y4 + a_y[6]*y4^2 )
            end
            if ( d1_p[nodesLocal[5],i] == Inf )
              d1_p[nodesLocal[5],i] =
                  a_x[1]      + a_x[2]*x5    + a_x[3]*y5 +
                  a_x[4]*x5^2 + a_x[5]*x5*y5 + a_x[6]*y5^2
              d2_p[nodesLocal[5],i] =
                  a_y[1]      + a_y[2]*x5    + a_y[3]*y5 +
                  a_y[4]*x5^2 + a_y[5]*x5*y5 + a_y[6]*y5^2
            else
              d1_p[nodesLocal[5],i] = .5*( d1_p[nodesLocal[5],i] +
                  a_x[1]      + a_x[2]*x5    + a_x[3]*y5 +
                  a_x[4]*x5^2 + a_x[5]*x5*y5 + a_x[6]*y5^2 )
              d2_p[nodesLocal[5],i] = .5*( d2_p[nodesLocal[5],i] +
                  a_y[1]      + a_y[2]*x5    + a_y[3]*y5 +
                  a_y[4]*x5^2 + a_y[5]*x5*y5 + a_y[6]*y5^2 )
            end

          elseif ( n==nodesLocal[3] )
            x5 = xLocal[5,1]
            y5 = xLocal[5,2]

            x6 = xLocal[6,1]
            y6 = xLocal[6,2]

            if ( d1_p[nodesLocal[5],i] == Inf )
              d1_p[nodesLocal[5],i] =
                  a_x[1]      + a_x[2]*x5    + a_x[3]*y5 +
                  a_x[4]*x5^2 + a_x[5]*x5*y5 + a_x[6]*y5^2
              d2_p[nodesLocal[5],i] =
                  a_y[1]      + a_y[2]*x5    + a_y[3]*y5 +
                  a_y[4]*x5^2 + a_y[5]*x5*y5 + a_y[6]*y5^2
            else
              d1_p[nodesLocal[5],i] = .5*( d1_p[nodesLocal[5],i] +
                  a_x[1]      + a_x[2]*x5    + a_x[3]*y5 +
                  a_x[4]*x5^2 + a_x[5]*x5*y5 + a_x[6]*y5^2 )
              d2_p[nodesLocal[5],i] = .5*( d2_p[nodesLocal[5],i] +
                  a_y[1]      + a_y[2]*x5    + a_y[3]*y5 +
                  a_y[4]*x5^2 + a_y[5]*x5*y5 + a_y[6]*y5^2 )
            end
            if ( d1_p[nodesLocal[6],i] == Inf )
              d1_p[nodesLocal[6],i] =
                  a_x[1]      + a_x[2]*x6    + a_x[3]*y6 +
                  a_x[4]*x6^2 + a_x[5]*x6*y6 + a_x[6]*y6^2
              d2_p[nodesLocal[6],i] =
                  a_y[1]      + a_y[2]*x6    + a_y[3]*y6 +
                  a_y[4]*x6^2 + a_y[5]*x6*y6 + a_y[6]*y6^2
            else
              d1_p[nodesLocal[6],i] = .5*( d1_p[nodesLocal[6],i] +
                  a_x[1]      + a_x[2]*x6    + a_x[3]*y6 +
                  a_x[4]*x6^2 + a_x[5]*x6*y6 + a_x[6]*y6^2 )
              d2_p[nodesLocal[6],i] = .5*( d2_p[nodesLocal[6],i] +
                  a_y[1]      + a_y[2]*x6    + a_y[3]*y6 +
                  a_y[4]*x6^2 + a_y[5]*x6*y6 + a_y[6]*y6^2 )
            end
          else
            fprintf("Problem identifying node at a vertex")
          end

        end
      end
    end
  end
#  Compare d1_p with the true gradient
#d_true=-2*x(:,1).*(1-x(:,2).^2)     % elliptic_2da example
#disp([d1_p d_true])

#  Compare d2_p with the true gradient
#d_true=-2*x(:,2).*(1-x(:,1).^2)
#disp([d2_p d_true])
  computeElementError = false  # later define this based on vargout?
  if ( computeElementError )
    #---------------------------------------------------------------------------
    #  Calculate the H1-seminorm of the error on each element
    #---------------------------------------------------------------------------
    eError = zeros(Float64,nElements,nFields)

    nQuadrature = 7;
    r,s,w   = twodQuadratureRule(nQuadrature)

    for nEl=1:nElements
      nodesLocal       = eConn[nEl,:]
      xLocal           = x[nodesLocal,:]
      uLocal           = u[nodesLocal,:]
      d1Local          = d1_p[nodesLocal,:]
      d2Local          = d2_p[nodesLocal,:]

      x_g,w_g,ϕ,ϕ_x,ϕ_y = twodShape(xLocal,r,s,w)

      # calculate the finite element derivatives and their projections at x_g
      dx_fe           = ϕ_x*uLocal
      dy_fe           = ϕ_y*uLocal

      dx_p            = ϕ*d1Local
      dy_p            = ϕ*d2Local

  # % exact derivatives as a test
  # dx_p = -2*x_g(:,1).*(1-x_g(:,2).^2)
  # dy_p = -2*x_g(:,2).*(1-x_g(:,1).^2)
      for i=1:nFields
        sqError     = (dx_fe[:,i]-dx_p[:,i]).^2 + (dy_fe[:,i]-dy_p[:,i]).^2
        eError[nEl,i] = sqError'*w_g

      end

    eError = sqrt.(eError)
    end
  end # computeElementError conditional

  return d1_p, d2_p, nodeList #, eError
end


text:
using Test

using CSV
readcsvmat(path;header=false, kwargs...) = convert(Matrix{Float64},CSV.read(path;header=header, kwargs...))

@testset "Lasso" begin

include("lasso.jl")
include("gammalasso.jl")
include("fusedlasso.jl")
include("trendfiltering.jl")
include("cross_validation.jl")

end