text:
module Registry

import ..Pkg
using ..Pkg: depots1, printpkgstyle, DEFAULT_IO, isdir_nothrow, pathrepr, pkg_server,
             GitTools, OFFLINE_MODE, UPDATED_REGISTRY_THIS_SESSION
using ..Pkg.PlatformEngines: download_verify_unpack, download
using UUIDs, LibGit2

include("registry_instance.jl")

mutable struct RegistrySpec
    name::Union{String,Nothing}
    uuid::Union{UUID,Nothing}
    url::Union{String,Nothing}
    # the path field can be a local source when adding a registry
    # otherwise it is the path where the registry is installed
    path::Union{String,Nothing}
end
RegistrySpec(name::String) = RegistrySpec(name = name)
RegistrySpec(;name::Union{String,Nothing}=nothing, uuid::Union{String,UUID,Nothing}=nothing,
              url::Union{String,Nothing}=nothing, path::Union{String,Nothing}=nothing) =
    RegistrySpec(name, isa(uuid, String) ? UUID(uuid) : uuid, url, path)

"""
    Pkg.Registry.add(url::String)
    Pkg.Registry.add(registry::RegistrySpec)

Add new package registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.add("General")
Pkg.Registry.add(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
Pkg.Registry.add(RegistrySpec(url = "https://github.com/JuliaRegistries/General.git"))
```
"""
add(reg::Union{String,RegistrySpec}; kwargs...) = add([reg]; kwargs...)
add(regs::Vector{String}; kwargs...) = add(RegistrySpec[RegistrySpec(name = name) for name in regs]; kwargs...)
function add(regs::Vector{RegistrySpec}; io::IO=DEFAULT_IO[])
    if isempty(regs)
        download_default_registries(io, only_if_empty = false)
    else
        download_registries(io, regs)
    end
end

const DEFAULT_REGISTRIES =
    RegistrySpec[RegistrySpec(name = "General",
                              uuid = UUID("23338594-aafe-5451-b93e-139f81909106"),
                              url = "https://github.com/JuliaRegistries/General.git")]

# Use the pattern
#
# registry_urls = nothing
# for ...
#     url, registry_urls = pkg_server_registry_url(uuid, registry_urls)
# end
#
# to query the pkg server at most once for registries.
pkg_server_registry_url(uuid::UUID, ::Nothing) =
    pkg_server_registry_url(uuid, pkg_server_registry_urls())

pkg_server_registry_url(uuid::UUID, registry_urls::Dict{UUID, String}) =
    get(registry_urls, uuid, nothing), registry_urls

pkg_server_registry_url(::Nothing, registry_urls) = nothing, registry_urls

function pkg_server_registry_urls()
    registry_urls = Dict{UUID, String}()
    server = pkg_server()
    server === nothing && return registry_urls
    tmp_path = tempname()
    download_ok = false
    try
        download("$server/registries", tmp_path, verbose=false)
        download_ok = true
    catch err
        @warn "could not download $server/registries" exception=err
    end
    download_ok || return registry_urls
    open(tmp_path) do io
        for line in eachline(io)
            if (m = match(r"^/registry/([^/]+)/([^/]+)$", line)) !== nothing
                uuid = UUID(m.captures[1])
                hash = String(m.captures[2])
                registry_urls[uuid] = "$server/registry/$uuid/$hash"
            end
        end
    end
    Base.rm(tmp_path, force=true)
    return registry_urls
end

pkg_server_url_hash(url::String) = Base.SHA1(split(url, '/')[end])

function download_default_registries(io::IO; only_if_empty::Bool = true)
    installed_registries = reachable_registries()
    # Only clone if there are no installed registries, unless called
    # with false keyword argument.
    if isempty(installed_registries) || !only_if_empty
        printpkgstyle(io, :Installing, "known registries into $(pathrepr(depots1()))")
        registries = copy(DEFAULT_REGISTRIES)
        for uuid in keys(pkg_server_registry_urls())
            if !(uuid in (reg.uuid for reg in registries))
                push!(registries, RegistrySpec(uuid = uuid))
            end
        end
        filter!(reg -> !(reg.uuid in installed_registries), registries)
        download_registries(io, registries)
        return true
    end
    return false
end

# Hacky way to make e.g. `registry add General` work.
function populate_known_registries_with_urls!(registries::Vector{RegistrySpec})
    known_registries = DEFAULT_REGISTRIES # TODO: Some way to add stuff here?
    for reg in registries, known in known_registries
        if reg.uuid !== nothing
            if reg.uuid === known.uuid
                reg.url = known.url
            end
        elseif reg.name !== nothing
            if reg.name == known.name
                named_regs = filter(r -> r.name == reg.name, known_registries)
                if !all(r -> r.uuid == first(named_regs).uuid, named_regs)
                    Pkg.Types.pkgerror("multiple registries with name `$(reg.name)`, please specify with uuid.")
                end
                reg.url = known.url
                reg.uuid = known.uuid
            end
        end
    end
end

registry_use_pkg_server() =
    !Sys.iswindows() || haskey(ENV, "JULIA_PKG_SERVER")

function download_registries(io::IO, regs::Vector{RegistrySpec}, depot::String=depots1())
    populate_known_registries_with_urls!(regs)
    registry_urls = nothing
    for reg in regs
        if reg.path !== nothing && reg.url !== nothing
            Pkg.Types.pkgerror("ambiguous registry specification; both url and path is set.")
        end
        # clone to tmpdir first
        mktempdir() do tmp
            url, registry_urls = pkg_server_registry_url(reg.uuid, registry_urls)
            if url !== nothing && registry_use_pkg_server()
                # download from Pkg server
                try
                    download_verify_unpack(url, nothing, tmp, ignore_existence = true)
                catch err
                    Pkg.Types.pkgerror("could not download $url")
                end
                tree_info_file = joinpath(tmp, ".tree_info.toml")
                hash = pkg_server_url_hash(url)
                write(tree_info_file, "git-tree-sha1 = " * repr(string(hash)))
            elseif reg.path !== nothing # copy from local source
                printpkgstyle(io, :Copying, "registry from `$(Base.contractuser(reg.path))`")
                mv(reg.path, tmp; force=true)
            elseif reg.url !== nothing # clone from url
                repo = GitTools.clone(io, reg.url, tmp; header = "registry from $(repr(reg.url))")
                LibGit2.close(repo)
            else
                Pkg.Types.pkgerror("no path or url specified for registry")
            end
            # verify that the clone looks like a registry
            if !isfile(joinpath(tmp, "Registry.toml"))
                Pkg.Types.pkgerror("no `Registry.toml` file in cloned registry.")
            end
            registry = Registry.RegistryInstance(tmp; parse_packages=false)
            # copy to `depot`
            regpath = joinpath(depot, "registries", registry.name)
            ispath(dirname(regpath)) || mkpath(dirname(regpath))
            if isfile(joinpath(regpath, "Registry.toml"))
                existing_registry = Registry.RegistryInstance(regpath; parse_packages=false)
                if registry.uuid == existing_registry.uuid
                    println(io,
                            "registry `$(registry.name)` already exist in `$(Base.contractuser(regpath))`.")
                else
                    throw(Pkg.Types.PkgError("registry `$(registry.name)=\"$(registry.uuid)\"` conflicts with " *
                        "existing registry `$(existing_registry.name)=\"$(existing_registry.uuid)\"`. " *
                        "To install it you can clone it manually into e.g. " *
                        "`$(Base.contractuser(joinpath(depot, "registries", registry.name*"-2")))`."))
                end
            else
                # if the dir doesn't exist, or exists but doesn't contain a Registry.toml
                mv(tmp, regpath, force=true)
                printpkgstyle(io, :Added, "registry `$(registry.name)` to `$(Base.contractuser(regpath))`")
            end
        end
    end
    return nothing
end

"""
    Pkg.Registry.rm(registry::String)
    Pkg.Registry.rm(registry::RegistrySpec)

Remove registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.rm("General")
Pkg.Registry.rm(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
rm(reg::Union{String,RegistrySpec}; kwargs...) = rm([reg]; kwargs...)
rm(regs::Vector{String}; kwargs...) = rm([RegistrySpec(name = name) for name in regs]; kwargs...)
function rm(regs::Vector{RegistrySpec}; io::IO=DEFAULT_IO[])
    for registry in find_installed_registries(io, regs)
        printpkgstyle(io, :Removing, "registry `$(registry.name)` from $(Base.contractuser(registry.path))")
        Base.rm(registry.path; force=true, recursive=true)
    end
    return nothing
end

# Search for the input registries among installed ones
function find_installed_registries(io::IO,
                                   needles::Union{Vector{Registry.RegistryInstance}, Vector{RegistrySpec}})
    haystack = reachable_registries(; parse_packages=false)
    output = Registry.RegistryInstance[]
    for needle in needles
        if needle.name === nothing && needle.uuid === nothing
            Pkg.Types.pkgerror("no name or uuid specified for registry.")
        end
        found = false
        for candidate in haystack
            if needle.uuid !== nothing
                if needle.uuid == candidate.uuid
                    push!(output, candidate)
                    found = true
                end
            elseif needle.name !== nothing
                if needle.name == candidate.name
                    named_regs = filter(r -> r.name == needle.name, haystack)
                    if !all(r -> r.uuid == first(named_regs).uuid, named_regs)
                        Pkg.Types.pkgerror("multiple registries with name `$(needle.name)`, please specify with uuid.")
                    end
                    push!(output, candidate)
                    found = true
                end
            end
        end
        if !found
            println(io, "registry `$(needle.name === nothing ? needle.uuid :
                                         needle.uuid === nothing ? needle.name :
                                         "$(needle.name)=$(needle.uuid)")` not found.")
        end
    end
    return output
end

"""
    Pkg.Registry.update()
    Pkg.Registry.update(registry::RegistrySpec)
    Pkg.Registry.update(registry::Vector{RegistrySpec})

Update registries. If no registries are given, update
all available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.update()
Pkg.Registry.update("General")
Pkg.Registry.update(RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"))
```
"""
update(reg::Union{String,RegistrySpec}; kwargs...) = update([reg]; kwargs...)
update(regs::Vector{String}; kwargs...) = update([RegistrySpec(name = name) for name in regs]; kwargs...)
function update(regs::Vector{RegistrySpec} = RegistrySpec[]; io::IO=DEFAULT_IO[], force::Bool=true)
    OFFLINE_MODE[] && return
    !force && UPDATED_REGISTRY_THIS_SESSION[] && return

    isempty(regs) && (regs = reachable_registries(; depots=depots1()))
    errors = Tuple{String, String}[]
    registry_urls = nothing
    for reg in unique(r -> r.uuid, find_installed_registries(io, regs); seen=Set{UUID}())
        let reg=reg
            regpath = pathrepr(reg.path)
            if reg.tree_info !== nothing
                printpkgstyle(io, :Updating, "registry at " * regpath)
                old_hash = reg.tree_info
                url, registry_urls = pkg_server_registry_url(reg.uuid, registry_urls)
                if url !== nothing && (new_hash = pkg_server_url_hash(url)) != old_hash
                    let new_hash = new_hash
                        # TODO: update faster by using a diff, if available
                        mktempdir() do tmp
                            try
                                download_verify_unpack(url, nothing, tmp, ignore_existence = true)
                            catch err
                                @error "could not download $url" exception=err
                            end
                            tree_info_file = joinpath(tmp, ".tree_info.toml")
                            hash = pkg_server_url_hash(url)
                            write(tree_info_file, "git-tree-sha1 = " * repr(string(new_hash)))
                            mv(tmp, reg.path, force=true)
                        end
                    end
                end
            elseif isdir(joinpath(reg.path, ".git"))
                printpkgstyle(io, :Updating, "registry at " * regpath)
                LibGit2.with(LibGit2.GitRepo(reg.path)) do repo
                    if LibGit2.isdirty(repo)
                        push!(errors, (regpath, "registry dirty"))
                        @goto done
                    end
                    if !LibGit2.isattached(repo)
                        push!(errors, (regpath, "registry detached"))
                        @goto done
                    end
                    if !("origin" in LibGit2.remotes(repo))
                        push!(errors, (regpath, "origin not in the list of remotes"))
                        @goto done
                    end
                    branch = LibGit2.headname(repo)
                    try
                        GitTools.fetch(io, repo; refspecs=["+refs/heads/$branch:refs/remotes/origin/$branch"])
                    catch e
                        e isa Pkg.Types.PkgError || rethrow()
                        push!(errors, (reg.path, "failed to fetch from repo"))
                        @goto done
                    end
                    ff_succeeded = try
                        LibGit2.merge!(repo; branch="refs/remotes/origin/$branch", fastforward=true)
                    catch e
                        e isa LibGit2.GitError && e.code == LibGit2.Error.ENOTFOUND || rethrow()
                        push!(errors, (reg.path, "branch origin/$branch not found"))
                        @goto done
                    end

                    if !ff_succeeded
                        try LibGit2.rebase!(repo, "origin/$branch")
                        catch e
                            e isa LibGit2.GitError || rethrow()
                            push!(errors, (reg.path, "registry failed to rebase on origin/$branch"))
                            @goto done
                        end
                    end
                    @label done
                end
            end
        end
    end
    if !isempty(errors)
        warn_str = "Some registries failed to update:"
        for (reg, err) in errors
            warn_str *= "\n    — $reg — $err"
        end
        @error warn_str
    end
    UPDATED_REGISTRY_THIS_SESSION[] = true
    return
end

"""
    Pkg.Registry.status()

Display information about available registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

# Examples
```julia
Pkg.Registry.status()
```
"""
function status(io::IO=DEFAULT_IO[])
    regs = reachable_registries()
    regs = unique(r -> r.uuid, regs; seen=Set{Union{UUID,Nothing}}())
    printpkgstyle(io, Symbol("Registry Status"), "")
    if isempty(regs)
        println(io, "  (no registries found)")
    else
        for reg in regs
            printstyled(io, " [$(string(reg.uuid)[1:8])]"; color = :light_black)
            print(io, " $(reg.name)")
            reg.url === nothing || print(io, " ($(reg.url))")
            println(io)
        end
    end
end

end # module


text:
@testset "femesh.jl" begin

end

ext: jl


text:
DataSubArray2{T,N,P<:DataValueArray,IV,LD} = SubArray{T,N,P,IV,LD}

@inline function Base.isnull(V::DataSubArray2, I::Int...)
    @boundscheck checkbounds(V, I...)
    @inbounds return V.parent.isnull[Base.reindex(V, V.indexes, I)...]
end

@inline function Base.values(V::DataSubArray2, I::Int...)
    @boundscheck checkbounds(V, I...)
    @inbounds return V.parent.values[Base.reindex(V, V.indexes, I)...]
end

FastDataSubArray2{T,N,P<:DataValueArray,IV} = SubArray{T,N,P,IV,true}

@inline function Base.isnull(V::FastDataSubArray2, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.isnull[V.first_index + V.stride1*i-1]
end

@inline function Base.values(V::FastDataSubArray2, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.values[V.first_index + V.stride1*i-1]
end

# We can avoid a multiplication if the first parent index is a Colon or UnitRange
FastDataContiguousSubArray2{T,N,P<:DataValueArray,I<:Tuple{Union{Colon, UnitRange}, Vararg{Any}}} = SubArray{T,N,P,I,true}

@inline function Base.isnull(V::FastDataContiguousSubArray2, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.isnull[V.first_index + i - 1]
end

@inline function Base.values(V::FastDataContiguousSubArray2, i::Int)
    @boundscheck checkbounds(V, i)
    @inbounds return V.parent.values[V.first_index + i - 1]
end

@generated function anynull{T, N, U<:DataValueArray}(S::SubArray{T, N, U})
    return quote
        isnull = view(S.parent.isnull, S.indexes...)
        @nloops $N i S begin
            (@nref $N isnull i) && (return true)
        end
        return false
    end
end


text:
# Julia program to Left Rotate a number by a specific bits.

# Left Rotate 'cnt' number of bits of the given number 'n'
function left_rotate_bits(n, cnt)
    cnt = cnt % 31
    while(cnt > 0)
        # Store the current MSB in a temporary variable
        msb = (n >> 31) & 1;
        # Left rotate the given number by one bit
        n = (n<<1);
        # Set the dropped MSB as the new LSB
        n = n | msb;
        # Decrement cnt
        cnt = cnt - 1
    end
    return n
end

print("Enter the number? ")
num = readline()
num = parse(Int, num)
print("How many bits do you want to rotate? ")
cnt = readline()
cnt = parse(Int, cnt)
left = left_rotate_bits(num, cnt)
print("The Left-rotated number is: $left")

"""
Time Complexity: O(n)
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

Enter the number? 39
How many bits do you want to rotate? 17
The Left-rotated number is: 5111808
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
using SemidefiniteOptInterface
const SDOI = SemidefiniteOptInterface

using Compat
using Compat.Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges

include("sdpa.jl")

const MOIU = MOI.Utilities
MOIU.@model(SDModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.PositiveSemidefiniteConeTriangle),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            ())

mock = SDOI.MockSDOptimizer{Float64}()
mock_optimizer = SDOI.SDOIOptimizer(mock, Float64)
@testset "supports_allocate_load" begin
    @test MOIU.supports_allocate_load(mock_optimizer, false)
    @test !MOIU.supports_allocate_load(mock_optimizer, true)
end
cached = MOIU.CachingOptimizer(SDModelData{Float64}(), mock_optimizer)
bridged = MOIB.full_bridge_optimizer(cached, Float64)
@testset "SolverName" begin
    @test MOI.get(mock,           MOI.SolverName()) == "MockSD"
    @test MOI.get(mock_optimizer, MOI.SolverName()) == "MockSD"
    @test MOI.get(cached,         MOI.SolverName()) == "MockSD"
    @test MOI.get(bridged,        MOI.SolverName()) == "MockSD"
end
config = MOIT.TestConfig(atol=1e-4, rtol=1e-4)

@testset "Unit" begin
    include("unit.jl")
end
@testset "MOI Continuous Linear" begin
    include("contlinear.jl")
end
@testset "MOI Continuous Conic" begin
    include("contconic.jl")
end


text:
struct MACAddress
    value::NTuple{6, UInt8}
end

function Base.show(io::IO, x::MACAddress)
    @printf(io, "%02x:%02x:%02x:%02x:%02x:%02x", x.value...)
end

isunicast(x::MACAddress) = (x.value[1] & 0x1) == 0
ismulticast(x::MACAddress) = (x.value[1] & 0x1) != 0

isuniversal(x::MACAddress) = (x.value[1] & 0x2) == 0
islocal(x::MACAddress) = (x.value[1] & 0x2) != 0

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols


text:
"""
    BnBModel

Stores attributes of stack used to solve BnB problem. Has the following fields:
* `Init_Box::Vector{Interval{Float64}}`:        stores initial interval box used
* `box::Vector{Vector{Interval{Float64}}}`      interval box storage stack
* `Init_Integer::Vector{Vector{Int64}}`         initial integer range
* `integers::Vector{Vector{Vector{Int64}}}`     integer range storage stack
* `LBD::Vector{Float64}`:                       lower bounds associated with each stack item
* `UBD::Vector{Float64}`:                       Upper bounds associated with each stack item
* `id::Vector{Int64}`:                          Node ID for each stack item
* `pos::Vector{Int64}`:                         Position in BnB Tree for each stack item
* `LBDg::Float64`:                              Global Lower Bound
* `UBDg::Float64`:                              Global Upper Bound
* `LBDg_hist::Vector{Float64}`:                 Value history LBD problem
* `UBDg_hist::Vector{Float64}`:                 Value history UBD problem
* `LBDgtime::Vector{Float64}`:                  Run time history LBD problem
* `UBDgtime::Vector{Float64}`:                  Run time history UBD problem
* `Pretime::Vector{Float64}`:                   Run time history preprocessing
* `Posttime::Vector{Float64}`:                  Run time history postprocessing
* `max_id::Int64`:                              Max node used
* `pstar::Vector{Interval{Float64}}`:           IntervalBox with solution
* `soln::Vector{Float64}`:                      Storage for solution
* `soln_val::Float64`:                          Solution value found
* `first_fnd::Bool`:                            Has a solution been found
* `feas_fnd::Bool`:                             Has a feasible point been found
* `first_num::Int64`:                           Iteration at which first solution found
* `lbcnt::Int64`:                               number of lower bounding problems solved
* `ubcnt::Int64`:                               number of upper bounding problems solved
"""
mutable struct BnBModel{V}
  Init_Box::Vector{V}
  box::Vector{Vector{V}} # interval box stack
  Init_Integer::Vector{Vector{Int64}}
  integers::Vector{Vector{Vector{Int64}}}
  LBD::Vector{Float64} # lower bounds associated with each stack
  UBD::Vector{Float64} # Upper bounds associated with each stack
  id::Vector{Int64}  # Node ID
  pos::Vector{Int64} # Position in BnB Tree
  LBDg::Float64 # Global Lower Bound
  UBDg::Float64 # Global Upper Bound
  LBDg_hist::Vector{Float64} # Value history LBD problem
  UBDg_hist::Vector{Float64} # Value history UBD problem
  LBDgtime::Vector{Float64} # Run time history LBD problem
  UBDgtime::Vector{Float64} # Run time history UBD problem
  Pretime::Vector{Float64} # Run time history LBD problem
  Posttime::Vector{Float64} # Run time history UBD problem
  max_id::Int64             # Max node used
  pstar::Vector{V}
  soln::Vector{Float64}
  soln_val::Float64
  first_fnd::Bool
  feas_fnd::Bool
  first_num::Int64
  lbcnt::Int64
  ubcnt::Int64
  lastgap::Float64
end

"""
    BnBModel(X::Vector{Interval{Float64}})

Initializes a `BnBModel` with `.Init_Box` = `X` and `.box` = `[X]`.
"""
BnBModel(X::Vector{V}) where {V} = BnBModel{V}(deepcopy(X),
                                            [deepcopy(X)],
                                            [[1]],
                                            [[[1]]],
                                            [-Inf],
                                            [Inf],
                                            [1],
                                            [1],
                                            -Inf,
                                            Inf,
                                            [-Inf],
                                            [Inf],
                                            [0.0],
                                            [0.0],
                                            [0.0],
                                            [0,0],
                                            1,
                                            deepcopy(X),
                                            [0.0],
                                            Inf,
                                            false,
                                            false,
                                            -1,
                                            0,
                                            0,
                                            Inf)
#BnBModel() = BnBModel{MCInterval{Float64}}([MCInterval(0.0,1.0)])


text:
using DataFrames, CSV, WAV, DSP, Images

class_to_i = Dict(
    "siren" => 0,
    "street_music" => 1,
    "drilling" => 2,
    "dog_bark" => 3,
    "children_playing" => 4,
    "gun_shot" => 5,
    "engine_idling" => 6,
    "air_conditioner" => 7,
    "jackhammer" => 8,
    "car_horn" => 9
)

# Read label dictionary : (sound id) => class
function gen_label_dict(labels_filename)
    # Read class data
    df = CSV.read(labels_filename, header=true)
    labels_dict = Dict()
    for i in 1:size(df, 1)
        labels_dict[df[i,1]] = df[i,2]
    end
    delete!(labels_dict, 8542) # Bad file
    delete!(labels_dict, 3190) # Spectrogram has 0 columns (weird)
    return labels_dict
end

function simple_labels(labels_dict)
    result = Array{Int}(undef, 0)
    for id in sort(collect(keys(labels_dict)))
        push!(result, class_to_i[labels_dict[id]])
    end
    return result
end

function create_spectrogram(wavfile)
    signal, sampleRate = wavread(wavfile)
    spec = spectrogram(signal[:,1], convert(Int, round(25e-3*sampleRate)),
                        convert(Int, round(10e-3*sampleRate)); window=hanning)
    t = time(spec)
    f = freq(spec)
    return log10.(power(spec))
end

function create_spectrogram_base(labels, n=200, m=100)
    spectres = Array{Any}(undef, 0)
    iter = 0
    for id in sort(collect(keys(labels)))
        spectr = create_spectrogram("train/Train/$id.wav")
        if size(spectr, 2) == 0
            @info("$id.wav spectrogram has 0 columns. Skipping...")
        else
            push!(spectres, imresize(sigmoid.(spectr), (n, m)))
        end
        iter += 1
        if iter % 300 == 0
            @info("Created $iter spectrograms")
        end
    end
    return spectres
end


text:
using TimeZones: Class

@testset "istimezone" begin
    # Invalidate the cache to ensure that `istimezone` works for non-loaded time zones.
    empty!(TimeZones.TIME_ZONE_CACHE)

    @test istimezone("Europe/Warsaw")
    @test istimezone("UTC+02")
    @test !istimezone("Europe/Camelot")
end

# Deserialization can cause us to have two immutables that are not using the same memory
@test TimeZone("Europe/Warsaw") === TimeZone("Europe/Warsaw")
@test tz"Africa/Nairobi" === TimeZone("Africa/Nairobi")

@testset "etcetera" begin
    # Note: In previous versions of TimeZones.jl the "etcetera" source file was not parsed
    # by default.
    @test !istimezone("Etc/GMT")
    @test !istimezone("Etc/GMT+12")
    @test !istimezone("Etc/GMT-14")

    @test istimezone("Etc/GMT", Class(:LEGACY))
    @test istimezone("Etc/GMT+12", Class(:LEGACY))
    @test istimezone("Etc/GMT-14", Class(:LEGACY))

    @test TimeZone("Etc/GMT", Class(:LEGACY)) == FixedTimeZone("Etc/GMT", 0)
    @test TimeZone("Etc/GMT+12", Class(:LEGACY)) == FixedTimeZone("Etc/GMT+12", -12 * 3600)
    @test TimeZone("Etc/GMT-14", Class(:LEGACY)) == FixedTimeZone("Etc/GMT-14", 14 * 3600)
end


text:
using MolecularIntegrals, Profile
bfs = build_basis(ethane,"cc-pvdz") #"6-31G")
@profview MolecularIntegrals.all_twoe_ints_rys(bfs)


text:
export generate

# TODO
# load triplet from flake or drv

# generate(mod::Module, args...; kwargs...) = generate(pkgdir(mod), args...; kwargs...)
# function generate(path::String; platforms::Vector{AbstractPlatform}=[HostPlatform()], include_lazy::Bool=false)
function generate(
    env_path::String = pwd(),
    platforms::Vector{<:AbstractPlatform}=[HostPlatform()],
    include_lazy::Bool=false,
    out_path::AbstractString="Depot.json",
    registries = Pkg.Registry.reachable_registries()
)
    ctx = Context(; registries, env=EnvCache(Types.projectfile_path(env_path)))
    pkgs = collect_packages(ctx, platforms, include_lazy)
    depot = resolve_depot(generate_depot(pkgs, ctx), ctx)
    # depot = generate_depot(pkgs, ctx)

    open(out_path, "w") do io_out
        if success(`prettier --help`)
            open(pipeline(`prettier --parser json`, io_out), write=true) do io_in
                JSON3.write(io_in, depot)
            end
        else
            JSON3.write(io_out, depot)
        end
    end
end

function printtypes(x::AbstractDict, indent="")
    for (k, v) in x
        print(indent, k, "=>")
        printtypes(v, indent * "  ")
    end
end
printtypes(x, indent) = println(indent, typeof(x))


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

include("distribution_functions.jl")
include("bat_sampler.jl")
include("HistogramAsUvDistribution.jl")
include("HistogramAsMvDistribution.jl")


text:
array = [1,2,3]
rand(array)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using TrajectoryOptimization
using StaticArrays
using LinearAlgebra
using Statistics
using BenchmarkTools
using ForwardDiff
using Test
using SparseArrays
const TO = TrajectoryOptimization

function make_mutable(conSet::ConstraintSet)
    n,m,P = size(conSet)
    N =  length(conSet.p)
    conSet_ = ConstraintSet(n,m,N)
    for con in conSet
        add_constraint!(conSet_, con.con, con.inds, mutable=true)
    end
    return conSet_
end
function make_mutable(prob::Problem)
    return Problem(prob, constraints = make_mutable(prob.constraints))
end

# Test entire solve
prob,opts = Problems.DubinsCar(:escape)
solver0 = iLQRSolver(copy(prob), opts)
solver1 = iLQRSolver2(copy(make_mutable(prob)), opts)
solver0 = ALTROSolver(copy(prob), opts, infeasible=true, R_inf=0.1)
solver1 = ALTROSolver(copy(make_mutable(prob)), opts, solver_uncon=iLQRSolver2, infeasible=true, R_inf=0.1)
# solve!(solver0)
# solve!(solver1)
# solver0.opts.projected_newton = false
# solver1.opts.projected_newton = false
b1 = benchmark_solve!(solver0)
b2 = benchmark_solve!(solver1)
judge(median(b2),median(b1))
@show median(b2).allocs
iterations(solver0)
iterations(solver1)
@show cost(solver0) ≈ cost(solver1)

model = get_model(solver1)
n,m = size(model)
x,u = rand(model)
z = KnotPoint(x,u,0.1)
∇f = zeros(n,n+m+1)
@btime TO.discrete_jacobian!($RK3, $∇f, $model, $z)
∇f = SizedMatrix{n,n+m+1}(∇f)
@btime TO.discrete_jacobian!($RK3, $∇f, $model, $z)

altro0 = ALTROSolver(copy(prob), opts)#, infeasible=true, R_inf=0.1)
altro1 = ALTROSolver(copy(make_mutable(prob)), opts, solver_uncon=iLQRSolver2)#, infeasible=true, R_inf=0.1)
al0 = altro0.solver_al
al1 = altro1.solver_al
initialize!(al0)
initialize!(al1)
solver0 = al0.solver_uncon
solver1 = al1.solver_uncon
solve!(solver0)
solve!(solver1)
iterations(solver0)
iterations(solver1)

# Step through
altro0 = ALTROSolver(copy(prob), opts)#, infeasible=true, R_inf=0.1)
altro1 = ALTROSolver(copy(make_mutable(prob)), opts, solver_uncon=iLQRSolver2)#, infeasible=true, R_inf=0.1)
al0 = altro0.solver_al
al1 = altro1.solver_al
solver0 = al0.solver_uncon
solver1 = al1.solver_uncon
rollout!(solver0)
rollout!(solver1)
states(solver0) ≈ states(solver1)
cost(solver0) ≈ cost(solver1)
TO.state_diff_jacobian!(solver0.G, solver0.model, solver0.Z)
TO.state_diff_jacobian!(solver1.G, solver1.model, solver1.Z)
# solver0.G ≈ solver1.G
TO.cost_expansion!(solver0.Q, solver0.G, solver0.obj, solver0.model, solver0.Z)
TO.cost_expansion!(solver1.Q, solver1.obj, solver1.Z)
TO.error_expansion!(solver1.Q, solver1.model, solver1.Z, solver1.G)

Q = [TO.error_expansion(Q, solver1.model) for Q in solver1.Q]
Qxx = [Q.xx for Q in Q]
Quu = [Q.uu for Q in Q]
Qux = [Q.ux for Q in Q]
Qx  = [Q.x for Q in Q]
Qu  = [Q.u for Q in Q]

Qxx ≈ solver0.Q.xx
Quu ≈ solver0.Q.uu
Qux ≈ solver0.Q.ux
Qx ≈ solver0.Q.x
Qu ≈ solver0.Q.u
TO.discrete_jacobian!(solver0.∇F, solver0.model, solver0.Z)
TO.dynamics_expansion!(solver1.D, solver1.model, solver1.Z)
∇F = [D.∇f for D in solver1.D]
∇F ≈ solver0.∇F

# step through backward pass
ΔV0 = backwardpass!(solver0)
ΔV1 = backwardpass!(solver1)
ΔV2 = TO.static_backwardpass!(solver1)
ΔV0 ≈ ΔV1 ≈ ΔV2
solver0.K ≈ solver1.K
solver0.d ≈ solver1.d

J0 = TO.forwardpass!(solver0, ΔV0, cost(solver0))
J1 = forwardpass!(solver1, ΔV1, cost(solver1))
J0 ≈ J1
states(solver0.Z̄) ≈ states(solver1.Z̄)

TO.copy_trajectories!(solver0)
TO.copy_trajectories!(solver1)

# Compare timing
al0 = AugmentedLagrangianSolver(copy(prob), opts)
al1 = AugmentedLagrangianSolver(copy(make_mutable(prob)), opts, solver_uncon=iLQRSolver2)
solver0 = al0.solver_uncon
solver1 = al1.solver_uncon
@btime rollout!($solver0)
@btime rollout!($solver1)
@btime TO.cost!($solver0.obj, $solver0.Z)
@btime TO.cost!($solver1.obj, $solver1.Z)
@btime TO.state_diff_jacobian!($solver0.G, $solver0.model, $solver0.Z)
@btime TO.state_diff_jacobian!($solver1.G, $solver1.model, $solver1.Z)
@btime TO.cost_expansion!($solver0.Q, $solver0.G, $solver0.obj, $solver0.model, $solver0.Z)
@btime TO.cost_expansion!($solver1.Q, $solver1.obj, $solver1.Z)
@btime TO.error_expansion!($solver1.Q, $solver1.model, $solver1.Z, $solver1.G)
@btime TO.discrete_jacobian!($solver0.∇F, $solver0.model, $solver0.Z)
@btime TO.dynamics_expansion!($solver1.D, $solver1.model, $solver1.Z)
@btime TO.error_expansion!($solver1.D, $solver1.model, $solver1.G)

@btime ΔV0 = backwardpass!($solver0)
@btime ΔV1 = backwardpass!($solver1)
@btime ΔV1 = TO.static_backwardpass!($solver1)

@btime TO.forwardpass!($solver0, $ΔV0, cost($solver0))
@btime TO.forwardpass!($solver1, $ΔV1, cost($solver1))
@btime rollout!($solver0, 1.0)
@btime rollout!($solver1, 1.0)

E1 = TO.GeneralExpansion{Float64}(SizedArray, 4,4,1)
E2 = TO.GeneralExpansion{Float64}(SizedArray, 4,4,1)
E1 = TO.SizedExpansion{Float64}(4,4,1)
E2 = TO.SizedExpansion{Float64}(4,4,1)

@btime TO._error_expansion!($E1,$E2,$G)
@btime TO.error_expansion!($solver.E, $solver.Q[1], $solver.model, $solver.Z[1], $(solver.G[1]))

struct MyType{T}


text:
struct Gene{T}
    param::Symbol
    val::T
    pow::Int
end

struct Genotype
    a::Vector{Gene}
    b::Vector{Gene}
end


text:
#=
  A/B test simple

  From ~/stan/ab_test_simple.stan
  """
                  mean      se_mean         sd         2.5%          25%
  rateA       0.09061974 0.0012142080 0.06150084   0.01159885   0.04405893
  rateB       0.07498042 0.0009625318 0.05068273   0.01016666   0.03781876
  rate_diff  -0.01563931 0.0015525787 0.07993122  -0.18981938  -0.05958019
  lp__      -14.93856837 0.0297626673 1.09616975 -17.86292780 -15.35567844
  """

  This Turing model
  """
  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

       rateA    0.0934    0.0612     0.0003    0.0014   1715.9223    1.0023      142.5185
       rateB    0.0763    0.0494     0.0002    0.0011   1683.9340    1.0002      139.8616
        diff   -0.0171    0.0784     0.0004    0.0016   2026.4826    1.0015      168.3125
     diffGt0    0.4111    0.4920     0.0025    0.0137   1085.2272    1.0029       90.1352

  Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       rateA    0.0126    0.0489    0.0820    0.1255    0.2451
       rateB    0.0094    0.0402    0.0671    0.1034    0.1997
        diff   -0.1845   -0.0619   -0.0147    0.0322    0.1342
     diffGt0    0.0000    0.0000    0.0000    1.0000    1.0000

  Credible interval for diff with mass 0.9: (-0.139303..0.114952)
  """

  See ~/stan/ab_test_simple.R and ~/stan/ab_test_simple.stan
      ~/blog/ab_test_simple.blog
      ~/webppl/ab_test_simple.wppl
=#

using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function ab_test_simple(nA=20,nB=25)
    # nA = 16 # Number of sent mail
    # nB = 16

    rateA ~ Beta(1,1) # priors
    rateB ~ Beta(1,1)

    1 ~ Binomial(nA,rateA) # likelihood
    1 ~ Binomial(nB,rateB)

    diff ~ Dirac(rateB - rateA)
    diffGt0 ~ Dirac(diff > 0)
end

nA = 20 # number of trial (sent mail) for A
nB = 25 # number of trial (sent mail) for B
model = ab_test_simple(nA,nB)

num_chns = 4
# chns = sample(model, Prior(), 1000)
# chns = sample(model, MH(), 40_000)

# chns = sample(model, PG(20), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)

display(chns)
# display(plot(chns))

credible_interval(chns,"diff",0.9)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module ColorTypes

using FixedPointNumbers
import FixedPointNumbers: floattype
using Base: @pure

const Fractional = Union{AbstractFloat, FixedPoint}

import Base: ==, isapprox, hash, convert, eltype, length, show, oneunit, zero, reinterpret, getindex
using Random
import Random: rand

## Types
export Fractional

export Colorant
export Color, TransparentColor, AlphaColor, ColorAlpha, AbstractRGB
export AbstractGray, Color3, TransparentGray, Transparent3, TransparentRGB, ColorantNormed
export AbstractAGray, AbstractGrayA, AbstractARGB, AbstractRGBA

export RGB, BGR, XRGB, RGBX
export HSV, HSB, HSL, HSI
export XYZ, xyY, LMS, Lab, LCHab, Luv, LCHuv
export DIN99, DIN99d, DIN99o
export YIQ, YCbCr

export Gray

export RGB24, ARGB32, Gray24, AGray32

# Note: the parametric TransparentColorColors are exported
# algorithmically, see `@make_alpha` in types.jl.

## Functions
export base_color_type, base_colorant_type, ccolor, color, color_type
export alphacolor, coloralpha
export alpha, red, green, blue, gray   # accessor functions that generalize to RGB24, etc.
export chroma, hue
export comp1, comp2, comp3
export mapc, reducec, mapreducec, gamutmax, gamutmin

include("types.jl")
include("traits.jl")
include("conversions.jl")
include("show.jl")
include("operations.jl")

Base.@deprecate_binding RGB1 XRGB
Base.@deprecate_binding RGB4 RGBX

"""
ColorTypes summary:

Type hierarchy:
```
                          Colorant{T,N}
             Color{T,N}                    TransparentColor{C,T,N}
     AbstractRGB{T}                  AlphaColor{C,T,N}  ColorAlpha{C,T,N}
```

Concrete types:
- `RGB`, `BGR`, `XRGB`, `RGBX`, `RGB24` are all subtypes of `AbstractRGB`

- `HSV`, `HSL`, `HSI`, `XYZ`, `xyY`, `Lab`, `LCHab`, `Luv`, `LCHuv`,
  `DIN99`, `DIN99d`, `DIN99o`, `LMS`, `YIQ`, `YCbCR` are subtypes of
  `Color{T,3}`

- Alpha-channel analogs in such as `ARGB` and `RGBA` for most of those
  types (with a few exceptions like `RGB24`, which has `ARGB32`)

- Grayscale types `Gray` and `Gray24` (subtypes of `Color{T,1}`), and
  the corresponding transparent types `AGray`, `GrayA`, and `AGray32`

- Trait functions `eltype`, `length`, `alphacolor`, `coloralpha`,
  `color_type`, `base_color_type`, `base_colorant_type`, `ccolor`

- Getters `red`, `green`, `blue`, `alpha`, `gray`, `comp1`, `comp2`, `comp3`

Use `?` to get more information about specific types or functions.
""" ColorTypes

if VERSION >= v"1.1" # work around https://github.com/JuliaLang/julia/issues/34121
    include("precompile.jl")
    _precompile_()
end

end # module


text:
using SpinGlassNetworks
using SpinGlassDynamics

using Test

my_tests = []

push!(my_tests,
    "simpleCIM.jl",
    "chimeraCIM.jl",
)

for my_test in my_tests
    include(my_test)
end


text:
using Documenter, ACME

makedocs(
    modules = [ACME],
    format = :html,
    sitename = "ACME.jl",
    pages = Any[
        "Home" => "index.md",
        "gettingstarted.md",
        "ug.md",
        "elements.md",
    ],
)

deploydocs(
    target = "build",
    repo = "github.com/HSU-ANT/ACME.jl.git",
    latest = "develop",
    deps = nothing,
    make = nothing,
    julia = "0.6",
)


text:
using WiltonInts84

using Test
using StaticArrays

include("num_quad.jl")

if !@isdefined record
    record = false
end

p1 = SVector(0.4,0.7000000000000001,0.0)
p2 = SVector(0.5,0.7000000000000001,0.0)
p3 = SVector(0.5,0.6000000000000001,0.0)
z  = SVector(0.0,0.0,1.0)

N = 2
T = eltype(Float64)
P = SVector{3,T}
J = Vector{SVector{N+3,T}}()
L = Vector{SVector{N+3,P}}()
I = similar(J)
K = similar(L)

h = 1.2

# circle through triangle vertex
c = SVector(0.020000000000000004,0.06,0.0)
r, R = 0.0, 0.8
q = contour(p1,p2,p3,c,r,R)
@test length(q.segments) == 0
@test length(q.arcs) == 0
@test length(q.circles) == 0
A,B = wiltonints(q,c, Val{N})
push!(I,SVector(A)); push!(K,SVector(B))
if record
    P,Q = dblquadints1(p1,p2,p3,c,Val{N},r,R)
    push!(J,P); push!(L,Q)
end


text:
@testset "Multiple Drivers" begin
    println("VITUIX MULTIPLE DRIVER TESTS")

    Hs = load_vituix_FIRs((joinpath(test_dir, "FIR", wav) for wav in ["T.wav", "M.wav", "W.wav"])...; Δt=-371.519ms, names=[:FIR_T, :FIR_M, :FIR_W])

    template_hor = Regex(raw"VituixCAD_PolarFR hor (.*)\.frd")
    template_ver = Regex(raw"VituixCAD_PolarFR ver (.*)\.frd")

    matchFIR = true

    test_set = [
            (src=[
                Dict(:name => "T", :spinorama => "T/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_T], :y => 0.0)),
                Dict(:name => "M", :spinorama => "M/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_M], :y => -0.118)),
                Dict(:name => "W", :spinorama => "W/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_W], :y => -0.348)),
                ], target="total/polar", freqrange=(0Hz, 20kHz), T=Dict(), max_mag=14, mae_mag=0.2, max_phi=180.0, mae_phi=2.5),
            (src=[
                Dict(:name => "T", :spinorama => "T/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_T], :r => 10.0, :t => 15.0)),
                Dict(:name => "M", :spinorama => "M/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_M], :x => 0.095, :y => -0.118, :z => 0.150, :r => 12.0, :t => -15.0)),
                Dict(:name => "W", :spinorama => "W/polar_noXO", :transform => Dict(:H => Hs[:, :FIR_W], :x => -0.15, :y => -0.348, :z => -0.035, :r => -90.0, :t => 3.0)),
                ], target="total_rotated_tilted2/polar", freqrange=(0Hz, 20kHz), T=Dict(), max_mag=21.0, mae_mag=0.7, max_phi=180.0, mae_phi=10.5),
        ]

    for matchFIR in [true, false]
        println("match FIR: $matchFIR")
        for e in test_set
            sources = String[]
            drivers = []
            for s in e.src
                push!(sources, s[:spinorama])
                push!(drivers, Dict(
                    :spinorama => load_spinorama_from_FRDs(joinpath(test_dir, s[:spinorama]), template_hor=template_hor, template_ver=template_ver),
                    :transform => s[:transform]
                    ))
            end
            println("  \"$(sources)\" -> \"$(e.target)\"")

            target_spinorama = load_spinorama_from_FRDs(joinpath(test_dir, e.target), template_hor=template_hor, template_ver=template_ver)

            print("\tA:")
            @test test_transform(drivers, target_spinorama, e.freqrange, compare_magdb, e.max_mag, e.mae_mag, matchFIR)
            print("\tΦ:")
            @test test_transform(drivers, target_spinorama, e.freqrange, compare_phasedeg, e.max_phi, e.mae_phi, matchFIR)
        end
    end
end


text:
# ---
# title: 926. Flip String to Monotone Increasing
# id: problem926
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/flip-string-to-monotone-increasing/description/>
# hidden: true
# ---
#
# A string of `'0'`s and `'1'`s is _monotone increasing_ if it consists of some
# number of `'0'`s (possibly 0), followed by some number of `'1'`s (also
# possibly 0.)
#
# We are given a string `S` of `'0'`s and `'1'`s, and we may flip any `'0'` to a
# `'1'` or a `'1'` to a `'0'`.
#
# Return the minimum number of flips to make `S` monotone increasing.
#
#
#
# **Example 1:**
#
#
#
#     Input: "00110"
#     Output: 1
#     Explanation: We flip the last digit to get 00111.
#
#
# **Example 2:**
#
#
#
#     Input: "010110"
#     Output: 2
#     Explanation: We flip to get 011111, or alternatively 000111.
#
#
# **Example 3:**
#
#
#
#     Input: "00011000"
#     Output: 2
#     Explanation: We flip to get 00000000.
#
#
#
#
# **Note:**
#
#   1. `1 <= S.length <= 20000`
#   2. `S` only consists of `'0'` and `'1'` characters.
#
#
## @lc code=start
using LeetCode

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


text:
# Unit testing
#=
TODO: Add more comparison tests?
TODO: Add bootstrapStatistic tests for p-dimensional statistics
TODO: Add comparison tests for statistics?
=#
using BlockBootstrap
using Test
using Random
using StatsBase
include("getTestData.jl")
include("testFunctions.jl")

function main()
    println("Starting Bootstrap tests")
    testAllCases()
    println("Tests finished")
end

main()


text:
using Pkg.Artifacts

const srm470_k412 = parse(Material, "(0.4535±0.0020)*SiO2+(0.1933±0.0020)*MgO+(0.1525±0.0020)*CaO+(0.0927±0.0020)*Al2O3+(0.0996±0.0020)*FeO",
    name="SRM-470 K412", density = 3.45, description="https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication260-74.pdf",
    pedigree="NIST SRM-470", conductivity=:Insulator)
const srm470_k411 = parse(Material,"(0.5430±0.0020)*SiO2+(0.1467±0.0020)*MgO+(0.1547±0.0020)*CaO+(0.1443±0.0020)*FeO",
    name="SRM-470 K411", pedigree="SRM-470", conductivity=:Insulator,
    description="https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nbsspecialpublication260-74.pdf")

# Materials from the "Mengason Mineral Mount I" SPI # 1025-AB
const mmm_albite = parse(
    Material,
    "0.1168*Na₂O+0.6876*SiO₂+0.1976*Al₂O₃+0.0017*K₂O+0.0023*CaO",
    name = "MMM1 Albite",
)
const mmm_almadine = parse(
    Material,
    "0.3915*SiO₂+0.0006*TiO₂+0.2271*Al₂O₃+0.2387*FeO+0.0055*MnO+0.1015*MgO+0.0401*CaO",
    name = "MMM1 Almadine",
)
const mmm_andradite = parse(
    Material,
    "0.002*MgO+0.0698*Al₂O₃+0.3423*SiO₂+0.3301*CaO+0.0023*TiO2₂+0.005*MnO+0.0222*FeO+0.2262*Fe₂O₃",
    name = "MMM1 Andradite",
)
const mmm_anhydrite =
    parse(Material, "0.0061*SrO+0.5834*SO₃+0.4113*CaO", name = "MMM1 Anhydrite")
const mmm_arsenopyrite =
    parse(Material, "0.3416*Fe+0.0002*Co+0.1928*S+0.4610*As", name = "MMM1 Arsenopyrite")
const mmm_barite = parse(Material, "0.3521*SO₃+0.6494*BaO", name = "MMM1 Barite")
const mmm_benitoite = parse(
    Material,
    "0.4534*SiO₂+0.1910*TiO₂+0.3494*BaO+0.0010*Na₂O",
    name = "MMM1 Benitotie",
)
const mmm_biotite = parse(
    Material,
    "0.3831*SiO₂+0.1598*Al₂O₃+0.1003*FeO+0.1967*MgO+0.1006*K₂O+0.0157*TiO₂+0.0385*H₂O",
    name = "MMM1 Biotite",
)
const mmm_bustamite = parse(
    Material,
    "0.4844*SiO₂+0.0023*MgO+0.2415*MnO+0.0813*FeO+0.1905*CaO+0.0027*ZnO",
    name = "MMM1 Bustamite",
)
const mmm_calcite =
    parse(Material, "0.5631*CaO+0.0005*SrO+0.4387*CO₂", name = "MMM1 Calcite")
const mmm_celestite =
    parse(Material, "0.5601*SrO+0.0009*CaO+0.4339*SO₃+0.0008*SiO2", name = "MMM1 Celestite")
const mmm_calcopyrite =
    parse(Material, "0.3048*Fe+0.3419*Cu+0.3514*S", name = "MMM1 Calcopyrite")
const mmm_chlorite = parse(
    Material,
    "0.3216*SiO₂+0.1487*Al₂O₃+0.0341*FeO+0.3355*MgO+0.0225*Cr₂O₃+0.0022*NiO+0.1376*H₂O",
    name = "MMM1 Chlorite",
)
const mmm_chromite = parse(
    Material,
    "0.0014*SiO₂+0.1678*MgO+0.2356*Al₂O₃+0.0016*NiO+0.0016*V₂O₃+0.0004*ZnO+0.0012*TiO₂+0.4542*Cr₂O₃+0.1291*FeO",
    name = "MMM1 Chromite",
)
const mmm_cinnabar = parse(Material, "0.8529*Hg+0.1350*S", name = "MMM1 Cinnabar")
const mmm_crocoite =
    parse(Material, "0.6914*PbO+0.3097*CrO₃+0.0009*SiO₂", name = "MMM1 Crocoite")
const mmm_cryolite = parse(Material, "0.6*NaF+0.4*AlF3", name = "MMM1 Cryolite")
const mmm_chromiumoxide = parse(Material, "Cr₂O₃", name = "MMM1 Chromium Oxide")
const mmm_diopside = parse(
    Material,
    "0.0041*NaO+0.5534*SiO₂+0.1776*MgO+0.0062*Al₂O₃+0.0083*FeO+0.2480*CaO",
    name = "MMM1 Diopside",
)
const mmm_dolomite =
    parse(Material, "0.3039*CaO+0.2166*MgO+0.0015*FeO+0.4768*CO₂", name = "MMM1 Dolomite")
const mmm_fluorapatite = parse(
    Material,
    "0.0299*F+0.0004*SiO₂+0.0028*SrO+0.4087*P₂O₅+0.0014*La₂O₃+0.0029*Ce₂O₃+0.0039*SO₃+0.5468*CaO+0.0006*Y₂O₃+0.0014*Nd₂O₃+0.0005*Pr₂O₃+0.0020*H₂O",
    name = "MMM1 Fluorapatite",
)
const mmm_galena = parse(Material, "0.8634*Pb+0.1351*S", name = "MMM1 Galena")
const mmm_hematite = parse(Material, "Fe2O3", name = "MMM1 Hematite")
const mmm_jadeite = parse(
    Material,
    "0.1069*Na+0.4713*O+0.0043*Mg+0.1202*Al+0.2776*S+0.0123*Ca+0.0072*Fe",
    name = "MMM1 Jadeite",
) # Not whats in the booklet
const mmm_kaersutite = parse(
    Material,
    "0.0205*H₂O+0.0262*Na₂O+0.1267*MgO+0.1191*Al₂O₃+0.4149*SiO₂+0.0101*K₂O+0.1140*CaO+0.0503*TiO₂+0.0016*MnO+0.1129*FeO+0.0014*F",
    name = "MMM1 Kaersutite",
)
const mmm_kyanite = parse(Material, "0.6292*Al₂O₃+0.3708*SiO₂", name = "MMM1 Kyanite")

const mmm_magnetite =
    parse(Material, "0.0003*SiO₂+0.0024*MnO+0.9934*Fe₃O₄", name = "MMM Magnetite")
const mmm_mgalspinel =
    parse(Material, "0.2833*MgO+0.7167*Al₂O₃", name = "MMM Magnesium Aluminum Spinel")
# The SPI book composition for Maganotantalite is clearly missing Fe and Ti
# const mmm_maganotantalite = parse(Material, "0.1410*MnO+0.0560*Nb₂O₅+0.8060*Ta₂O₅", name = "MMM Manganotantalite") # Book value
# This one is based on measurement by NWMR quantified using DTSA-II
const mmm_maganotantalite = parse(
    Material,
    "0.2064*O+0.0035*Ti+0.1007*Mn+0.0185*Fe+0.0986*Nb+0.5744*Ta",
    name = "MMM Maganotantalite",
)

const mmm_molybdenite = parse(Material, "0.3984*S+0.5929*Mo", name = "MMM Molybdenite")
const mmm_monazite = parse(
    Material,
    "0.017*SiO₂+0.2704*P₂O₅+0.0090*CaO+0.0200*Y₂O₃+0.0948*La₂O₃+0.2504*Ce₂O₃+0.0297*Pr₂O₃+0.1100*Nd₂O₃+0.0320*Sm₂O₃+0.0012*Eu₂O₃+0.0256*Gd₂O₃+0.0015*Er₂O₃+0.0030*PbO+0.1180*ThO₂+0.0020*UO₂",
    name = "MMM Monazite (REE,Th,Ca)(P,Si)O₄",
)
const mmm_obsidian = parse(
    Material,
    "0.7629*SiO₂+0.001*TiO₂+0.1303*Al₂O₃+0.0076*FeO+0.0011*MgO+0.0088*CaO+0.0374*Na₂O+0.0005*MnO+0.0431*K₂O",
    name = "MMM Obsidian Na,K,Al,Fe silicate glass",
)
const mmm_olivine = parse(
    Material,
    "0.0817*FeO+0.5044*MgO+0.4101*SiO₂+0.00354*NiO+0.00115*MnO",
    name = "MMM Olivine (Mg,Fe)₂SiO₄",
)
const mmm_orthoclase = parse(
    Material,
    "0.1549*K₂O+0.0091*Na₂O+0.0005*CaO₂+0.00201*Fe₂O₃+0.1674*Al₂O₃+0.6480*SiO₂",
    name = "MMM Orthoclase KAlSi₃O₈",
)
const mmm_pentlandite = parse(
    Material,
    "0.3088*Fe+0.3561*Ni+0.0043*Co+0.3289*S",
    name = "MMM Pentlandite (Fe,Ni)₉S₈",
)
const mmm_phlogpite = parse(
    Material,
    "0.0432*H₂O+0.2898*MgO+0.1222*Al₂O₃+0.4320*SiO₂+0.1129*K₂O",
    name = "MMM Phlogpite KMg₃AlSi₃O₁₀(OH)₂",
)
const mmm_plagioclase = parse(
    Material,
    "0.0436*Na₂O+0.5312*SiO₂+0.0010*MgO+0.2935*Al₂O₃+0.0034*FeO+0.0024*K₂O+0.1193*CaO",
    name = "MMM Plagioclase (Labradorite) (Ca,Na)(Al,Si)₄O₈",
)
const mmm_pollucite = parse(
    Material,
    "0.0011*Li₂O+0.0203*Na₂O+0.1712*Al₂O₃+0.4746*SiO₂+0.0018*K₂+0.0018*Rb₂O+0.3422*Cs₂O",
    name = "MMM Pollucite (CsSi₂AlO₆)",
)
const mmm_pyrite = parse(Material, "0.4630*Fe+0.5341*S", name = "MMM Pyrite")
const mmm_pyrope = parse(
    Material,
    "0.4143*SiO₂+0.0050*TiO₂+0.2158*Al₂O₃+0.0206*Cr₂O₃+0.0876*FeO+0.0031*MnO+0.2035*MgO+0.0435*CaO",
    name = "MMM Pyrope Mg₃Al₂So₃O₁₂",
)
const mmm_quartz = parse(Material, "SiO2", name = "MMM Quartz")
const mmm_rhodonite = parse(
    Material,
    "0.0434*FeO+0.4678*SiO2+0.0195*MgO+0.4230*MnO+0.0463*CaO",
    name = "MMM Rhodonite MnSiO₃",
)
const mmm_scheelite =
    parse(Material, "0.1951*CaO+0.002*MoO₃+0.8032*WO₃", name = "MMM Scheelite CaWO₄")
const mmm_spessartine = parse(
    Material,
    "0.2026*Al₂O₃+0.3629*SiO₂+0.0102*CaO+0.0014*TiO₂+0.4102*MnO+0.012*FeO+0.0006*SnO₂",
    name = "MMM Spessartine Mn₃Al₂Si₃O₁₂",
)
const mmm_sphene = parse(
    Material,
    "0.0136*Al₂O₃+0.3083*SiO₂+0.2882*CaO+0.3781*TiO₂+0.0005*MnO+0.0066*FeO",
    name = "MMM Sphene CaTiSiO₅",
)
const mmm_spodumene = parse(
    Material,
    "0.0798*Li₂O+0.0012*Na₂O+0.2728*Al₂O₃+0.6423*SiO₂",
    name = "MMM Spodumene LiAlSi₂O₆",
)
const mmm_stibnite = parse(Material, "Sb₂S₃", name = "MMM Stibnite Sb₂S₃")
const mmm_strontium_titanate =
    parse(Material, "SrTiO₃", name = "MMM Strontium Titanate SrTiO₃")
const mmm_tugtupite = parse(
    Material,
    "0.5124*SiO₂+0.1089*Al₂O₃+0.2563*Na₂O+0.07574*ClO+0.0522*BeO",
    name = "MMM Tugtupite Na₄BeAlSi₄O₁₂Cl",
)
const mmm_uvarovite = parse(
    Material,
    "0.0005*MgO+0.0668*Al₂O₃+0.3623*SiO₂+0.3489*CaO+0.0158*TiO₂+0.0029*V₂O₅+0.1903*Cr₂O₃+0.0124*Fe₂O₃",
    name = "MMM Uvarovite Ca₃Cr₂Si₃O₁₂",
)
const mmm_willemite = parse(
    Material,
    "0.6612*ZnO+0.2830*SiO₂+0.0595*MnO+0.0002*FeO",
    name = "MMM Willemite (Zn,Mn)₂SiO₄",
)
const mmm_wollastonite =
    parse(Material, "0.5173*SiO₂+0.4828*CaO", name = "MMM Wollastonite CaSiO₈")
const mmm_zircon = parse(
    Material,
    "0.3265*SiO₂+0.0009*P₂O₅+0.0006*Y₂O₃+0.6615*ZrO₂+0.00106*HfO₂",
    name = "MMM Zircon ZrSiO₄",
)
const mmm_carbon = parse(Material, "C", name = "MMM Carbon C")
const mmm_copper = parse(Material, "Cu", name = "MMM Copper Cu")

const mengason_mineral_mount1 = (
    mmm_albite,
    mmm_almadine,
    mmm_andradite,
    mmm_anhydrite,
    mmm_arsenopyrite,
    mmm_barite,
    mmm_benitoite,
    mmm_biotite,
    mmm_bustamite,
    mmm_calcite,
    mmm_celestite,
    mmm_calcopyrite,
    mmm_chlorite,
    mmm_chromite,
    mmm_cinnabar,
    mmm_crocoite,
    mmm_cryolite,
    mmm_calcite,
    mmm_chromiumoxide,
    mmm_diopside,
    mmm_dolomite,
    mmm_fluorapatite,
    mmm_galena,
    mmm_hematite,
    mmm_jadeite,
    mmm_kaersutite,
    mmm_kyanite,
)

const mengason_mineral_mount2 = (
    mmm_magnetite,
    mmm_mgalspinel,
    mmm_maganotantalite,
    mmm_molybdenite,
    mmm_monazite,
    mmm_obsidian,
    mmm_olivine,
    mmm_orthoclase,
    mmm_pentlandite,
    mmm_phlogpite,
    mmm_plagioclase,
    mmm_pollucite,
    mmm_pyrite,
    mmm_pyrope,
    mmm_quartz,
    mmm_rhodonite,
    mmm_scheelite,
    mmm_spessartine,
    mmm_sphene,
    mmm_spodumene,
    mmm_stibnite,
    mmm_strontium_titanate,
    mmm_tugtupite,
    mmm_uvarovite,
    mmm_willemite,
    mmm_wollastonite,
    mmm_zircon,
    mmm_carbon,
    mmm_copper,
)

function createmineralartifact()
    # This is the path to the Artifacts.toml we will manipulate
    artifacts_toml = joinpath(@__DIR__, "Artifacts.toml")

    # Query the `Artifacts.toml` file for the hash bound to the name "iris"
    # (returns `nothing` if no such binding exists)
    mineral_hash = artifact_hash("rplraw", artifacts_toml)

    # If the name was not bound, or the hash it was bound to does not exist, create it!
    if mineral_hash === nothing || !artifact_exists(mineral_hash)
        # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
        mineral_hash = create_artifact() do artifact_dir
            print("Downloading mineral database.")
            tarball = joinpath(artifact_dir, "mineral.tar.gz")
            download(
                "https://drive.google.com/uc?export=download&id=1Ackbz0YtaliQNCdZmPfj9uWTwhVNBzy8",
                tarball,
            )
            Pkg.probe_platform_engines!()
            Pkg.unpack(tarball, artifact_dir, verbose = true)
            rm(tarball)
        end
        # Now bind that hash within our `Artifacts.toml`.  `force = true` means that if it already exists,
        # just overwrite with the new content-hash.  Unless the source files change, we do not expect
        # the content hash to change, so this should not cause unnecessary version control churn.
        bind_artifact!(artifacts_toml, "rplraw", mineral_hash)
    end

    # Get the path of the iris dataset, either newly created or previously generated.
    # this should be something like `~/.julia/artifacts/dbd04e28be047a54fbe9bf67e934be5b5e0d357a`
    return artifact_path(mineral_hash)
end

function loadmineraldata(parseit::Bool = false)::DataFrame
    minpath = createmineralartifact()
    res = CSV.File(joinpath(minpath, "RRUFF_Export_20191025_022204.csv")) |> DataFrame
    if parseit
        function parseelm(str)
            if str == "Ln" # Lanthanide
                return Set{Element}(elements[z(n"La"):z(n"Lu")])
            elseif str == "An" # Actinide
                return Set{Element}(elements[z(n"Ac"):z(n"U")])
            elseif str == "REE" # Rare-earth element
                return Set{Element}([
                    n"Ce",
                    n"Dy",
                    n"Er",
                    n"Eu",
                    n"Gd",
                    n"Ho",
                    n"La",
                    n"Lu",
                    n"Nd",
                    n"Pr",
                    n"Pm",
                    n"Sm",
                    n"Sc",
                    n"Tb",
                    n"Tm",
                    n"Yb",
                    n"Y",
                ])
            else
                return Set{Element}([parse(Element, str)])
            end
        end
        function parseelms(str)
            if !ismissing(str)
                return mapreduce(
                    parseelm,
                    union,
                    filter(s -> length(s) > 0, split(str, c -> isspace(c))),
                    init = Set{Element}(),
                )
            else
                return Set{Element}()
            end
        end
        function matormissing(row)
            str = row["IMA Chemistry (plain)"]
            try
                # The formula with '+' represent valences not sums
                if isnothing(findfirst(c -> c == '+', str))
                    m = parse(Material, str, name = row["Mineral Name"])
                else
                    return missing
                end
            catch
                return missing
            end
        end
        res[:, :Elements] .= parseelms.(res[:, "Chemistry Elements"])
        res[:, :Material] .= matormissing.(eachrow(res))
    end
    return res
end

"""
    loadsmithsoniandata(; clean=false)

Load compositional data associated with the Smithsonian Microbeam Standards data set as a DataFrame. Setting clean=true will replace "<0.XXX" with 0.0,
replace "missing" with 0.0 and parse string values as Float64.
The data source is https://naturalhistory.si.edu/research/mineral-sciences/collections-overview/reference-materials/smithsonian-microbeam-standards
"""
function loadsmithsoniandata(; clean = false)
    tmp =
        CSV.File(
            joinpath(@__DIR__, "..", "data", "smithsonian_microbeam_standards.csv"),
            header = 2,
        ) |> DataFrame
    if clean
        res = copy(tmp[:, 1:4])
        colnames = propertynames(tmp)
        val(i) =
            ismissing(i) ? 0.0 :
            (i isa String ? (startswith(i, "<") ? 0.0 : parse(Float64, i)) : i)
        for i = 5:ncol(tmp)
            res[:, colnames[i]] = val.(tmp[:, i])
        end
        return res
    else
        return tmp
    end
end

"""
    parsedsmithsoniandata()::Dict{String, Material}

Converts the data from loadsmithsoniandata() into Material representation.
"""
function parsedsmithsoniandata()::Dict{String,Material}
    data = loadsmithsoniandata(clean = true)
    pm(cn) =
        cn == "H2O-" ? "H2O" : (cn == "nB2O5" ? "B2O5" : (cn == "REE2O3" ? "La2O3" : cn))
    cols = Dict{String,String}(colname => pm(colname) for colname in names(data)[5:67])
    res = Dict{String,Material}()
    for r in eachrow(data)
        str = ""
        for (cn, mat) in cols
            if r[cn] > 0.0
                str *= "+$(r[cn]/100.0)*$mat"
            end
        end
        mat = parse(Material, str[2:end], name = r[:Name])
        mat[:CatalogNumber] = r["Catalog Number"]
        mat[:EZID] = r[:EZID]
        mat[:ISGN] = r[:IGSN]
        if r[:REE2O3] > 0.0
            mat[:Note] = "Ambiguous REE2O3 replaced with La2O3."
        end
        res[r[:Name]] = mat
    end
    return res
end


text:
# x
const dummy_single_variable   = (x::Vector{MOI.VariableIndex}) -> MOI.SingleVariable(x[1])
# x₁, x₂
const dummy_vectorofvariables = (x::Vector{MOI.VariableIndex}) -> MOI.VectorOfVariables(x)
# 1.0 * x + 0.0
const dummy_scalar_affine     = (x::Vector{MOI.VariableIndex}) -> MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0)
# 1.0 * x + 1.0 * x^2 + 0.0
const dummy_scalar_quadratic  = (x::Vector{MOI.VariableIndex}) -> MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(1.0, x), MOI.ScalarQuadraticTerm.(1.0, x, x), 0.0)
# x₁ +    + 0.0
#    + x₂ + 0.0
const dummy_vector_affine     = (x::Vector{MOI.VariableIndex}) -> MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1:length(x), MOI.ScalarAffineTerm.(1.0, x)), zeros(Float64, length(x)))
# x₁ +    + x₁^2        + 0.0
#    + x₂ +      + x₂^2 + 0.0
const dummy_vector_quadratic  = (x::Vector{MOI.VariableIndex}) -> MOI.VectorQuadraticFunction(
    MOI.VectorAffineTerm.(1:length(x), MOI.ScalarAffineTerm.(1.0, x)),           # affine component
    MOI.VectorQuadraticTerm.(1:length(x), MOI.ScalarQuadraticTerm.(1.0, x, x)),  # affine component
    zeros(Float64, length(x))                                                    # constant term
)

const BasicConstraintTests = Dict(
    (MOI.SingleVariable, MOI.LessThan{Float64})    => ( dummy_single_variable, 1, MOI.LessThan(1.0) ),
    (MOI.SingleVariable, MOI.GreaterThan{Float64}) => ( dummy_single_variable, 1, MOI.GreaterThan(1.0) ),
    (MOI.SingleVariable, MOI.EqualTo{Float64})     => ( dummy_single_variable, 1, MOI.EqualTo(1.0) ),
    (MOI.SingleVariable, MOI.Interval{Float64})    => ( dummy_single_variable, 1, MOI.Interval(1.0, 2.0) ),

    (MOI.SingleVariable, MOI.ZeroOne)                 => ( dummy_single_variable, 1, MOI.ZeroOne() ),
    (MOI.SingleVariable, MOI.Integer)                 => ( dummy_single_variable, 1, MOI.Integer() ),
    (MOI.SingleVariable, MOI.Semicontinuous{Float64}) => ( dummy_single_variable, 1, MOI.Semicontinuous(1.0, 2.0) ),
    (MOI.SingleVariable, MOI.Semiinteger{Float64})    => ( dummy_single_variable, 1, MOI.Semiinteger(1.0, 2.0) ),

    (MOI.VectorOfVariables, MOI.SOS1{Float64}) => ( dummy_vectorofvariables, 2, MOI.SOS1([1.0, 2.0]) ),
    (MOI.VectorOfVariables, MOI.SOS2{Float64}) => ( dummy_vectorofvariables, 2, MOI.SOS2([1.0, 2.0]) ),
    (MOI.VectorOfVariables, MOI.Reals)         => ( dummy_vectorofvariables, 2, MOI.Reals(2) ),
    (MOI.VectorOfVariables, MOI.Zeros)         => ( dummy_vectorofvariables, 2, MOI.Zeros(2) ),
    (MOI.VectorOfVariables, MOI.Nonpositives)  => ( dummy_vectorofvariables, 2, MOI.Nonpositives(2) ),
    (MOI.VectorOfVariables, MOI.Nonnegatives)  => ( dummy_vectorofvariables, 2, MOI.Nonnegatives(2) ),

    (MOI.VectorOfVariables, MOI.SecondOrderCone)        => ( dummy_vectorofvariables, 3, MOI.SecondOrderCone(3) ),
    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => ( dummy_vectorofvariables, 3, MOI.RotatedSecondOrderCone(3) ),
    (MOI.VectorOfVariables, MOI.GeometricMeanCone)      => ( dummy_vectorofvariables, 3, MOI.GeometricMeanCone(3) ),
    (MOI.VectorOfVariables, MOI.ExponentialCone)        => ( dummy_vectorofvariables, 3, MOI.ExponentialCone() ),
    (MOI.VectorOfVariables, MOI.DualExponentialCone)    => ( dummy_vectorofvariables, 3, MOI.DualExponentialCone() ),
    (MOI.VectorOfVariables, MOI.PowerCone)              => ( dummy_vectorofvariables, 3, MOI.PowerCone(2.0) ),
    (MOI.VectorOfVariables, MOI.DualPowerCone)          => ( dummy_vectorofvariables, 3, MOI.DualPowerCone(2.0) ),

    (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle) => ( dummy_vectorofvariables,  7, MOI.PositiveSemidefiniteConeTriangle(3) ),
    (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeSquare)   => ( dummy_vectorofvariables, 10, MOI.PositiveSemidefiniteConeSquare(3) ),
    (MOI.VectorOfVariables, MOI.LogDetConeTriangle)               => ( dummy_vectorofvariables,  8, MOI.LogDetConeTriangle(3) ),
    (MOI.VectorOfVariables, MOI.LogDetConeSquare)                 => ( dummy_vectorofvariables, 11, MOI.LogDetConeSquare(3) ),
    (MOI.VectorOfVariables, MOI.RootDetConeTriangle)              => ( dummy_vectorofvariables,  7, MOI.RootDetConeTriangle(3) ),
    (MOI.VectorOfVariables, MOI.RootDetConeSquare)                => ( dummy_vectorofvariables, 10, MOI.RootDetConeSquare(3) ),

    (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => ( dummy_scalar_affine, 1, MOI.LessThan(1.0) ),
    (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => ( dummy_scalar_affine, 1, MOI.GreaterThan(1.0) ),
    (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})     => ( dummy_scalar_affine, 1, MOI.EqualTo(1.0) ),
    (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})    => ( dummy_scalar_affine, 1, MOI.Interval(1.0, 2.0) ),

    (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})    => ( dummy_scalar_quadratic, 1, MOI.LessThan(1.0) ),
    (MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}) => ( dummy_scalar_quadratic, 1, MOI.GreaterThan(1.0) ),
    (MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64})     => ( dummy_scalar_quadratic, 1, MOI.EqualTo(1.0) ),
    (MOI.ScalarQuadraticFunction{Float64}, MOI.Interval{Float64})    => ( dummy_scalar_quadratic, 1, MOI.Interval(1.0, 2.0) ),

    (MOI.VectorAffineFunction{Float64}, MOI.Reals)        => ( dummy_vector_affine, 2, MOI.Reals(2) ),
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros)        => ( dummy_vector_affine, 2, MOI.Zeros(2) ),
    (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) => ( dummy_vector_affine, 2, MOI.Nonpositives(2) ),
    (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_affine, 2, MOI.Nonnegatives(2) ),

    (MOI.VectorQuadraticFunction{Float64}, MOI.Reals)        => ( dummy_vector_quadratic, 2, MOI.Reals(2) ),
    (MOI.VectorQuadraticFunction{Float64}, MOI.Zeros)        => ( dummy_vector_quadratic, 2, MOI.Zeros(2) ),
    (MOI.VectorQuadraticFunction{Float64}, MOI.Nonpositives) => ( dummy_vector_quadratic, 2, MOI.Nonpositives(2) ),
    (MOI.VectorQuadraticFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_quadratic, 2, MOI.Nonnegatives(2) )
)
"""
    basic_constraint_tests(model::MOI.ModelLike, config::TestConfig;
        delete                  = true,
        get_constraint_function = true,
        get_constraint_set      = true,
        include                 = Any[],
        exclude                 = Any[]
        )

Test basic constraint functionality for all function-in-set combinations that
are supported by `model`.

See also `basic_constraint_test_helper`.

If `delete = true`, it will test the deletion of constraints.
If `get_constraint_function = true`, it will test the getting of `ConstraintFunction`.
If `get_constraint_set = true`, it will test the getting of `ConstraintSet`.

`include` and `exclude` can be used to run a subset of the tests, although only
one can be used in each call.

### Examples

    basic_constraint_tests(model, config; exclude = [
        (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    ])

    basic_constraint_tests(model, config;
        delete = false,
        include = [
            (MOI.VectorOfVariables, MOI.SOS2{Float64})
        ]
    )
"""
function basic_constraint_tests(model::MOI.ModelLike, config::TestConfig;
    delete                  = true,
    get_constraint_function = true,
    get_constraint_set      = true,
    include                 = Any[],
    exclude                 = Any[]
    )
    if length(include) > 0 && length(exclude) > 0
        error("Don't use `include` and `exclude` together.")
    end
    test_keys = length(include) > 0 ? include : Iterators.filter(x->!(x in exclude), keys(BasicConstraintTests))
    for (F,S) in test_keys
        if MOI.supports_constraint(model, F, S)
            @testset "$(F)-$(S)" begin
                (cf, N, set) = BasicConstraintTests[(F,S)]
                basic_constraint_test_helper(model, config, cf, set, N;
                    delete                  = delete,
                    get_constraint_function = get_constraint_function,
                    get_constraint_set      = get_constraint_set
                )
            end
        end
    end
end

"""
    basic_constraint_test_helper(model::MOI.ModelLike, config::TestConfig, func::Function, set::MOI.AbstractSet, N::Int;
        delete::Bool                  = true,
        get_constraint_function::Bool = true,
        get_constraint_set::Bool      = true
    )

A helper function for `basic_constraint_tests`.

This function tests the basic functionality of the constraint type `F`-in-`S`
where `S = typeof(set)`, and `func` is a function that takes a vector of `N`
variables and returns a `MOI.AbstractFunction` of type `F`.

If `delete = true`, it will test the deletion of constraints.
If `get_constraint_function = true`, it will test the getting of `ConstraintFunction`.
If `get_constraint_set = true`, it will test the getting of `ConstraintSet`.

### Example

    basic_constraint_test_helper(model, config,
        (x) -> MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.LessThan(1.0),
        1;
        delete=false
    )

"""
function basic_constraint_test_helper(model::MOI.ModelLike, config::TestConfig, func::Function, set::MOI.AbstractSet, N::Int=1;
    delete::Bool                  = true,
    get_constraint_function::Bool = true,
    get_constraint_set::Bool      = true
    )
    MOI.empty!(model)
    x = MOI.add_variables(model, N)
    constraint_function = func(x)
    F, S = typeof(constraint_function), typeof(set)

    @test MOI.supports_constraint(model, F, S)

    @testset "add_constraint" begin
        @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 0
        c = MOI.add_constraint(model, constraint_function, set)
        @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 1

        @testset "ConstraintName" begin
            @test MOI.get(model, MOI.ConstraintName(), c) == ""
            @test MOI.supports(model, MOI.ConstraintName(), typeof(c))
            MOI.set(model, MOI.ConstraintName(), c, "c")
            @test MOI.get(model, MOI.ConstraintName(), c) == "c"
        end

        if get_constraint_function
            @testset "ConstraintFunction" begin
                @test MOI.get(model, MOI.ConstraintFunction(), c) ≈ constraint_function
            end
        end
        if get_constraint_set
            @testset "ConstraintSet" begin
                @test MOI.get(model, MOI.ConstraintSet(), c) == set
            end
        end
    end

    @testset "ListOfConstraintIndices" begin
        c_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        @test length(c_indices) == MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 1
    end

    if F != MOI.SingleVariable
        # We can't add multiple single variable constraints as these are
        # interpreted as bounds etc.
        @testset "add_constraints" begin
            n = MOI.get(model, MOI.NumberOfConstraints{F,S}())
            cc = MOI.add_constraints(model,
                [constraint_function, constraint_function],
                [set, set]
            )
            @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 3
            @test length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}())) == 3
        end
    end

    @testset "is_valid" begin
        c_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        @test all(MOI.is_valid.(model, c_indices))
    end

    if delete
        @testset "delete" begin
            c_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            MOI.delete(model, c_indices[1])
            @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == length(c_indices)-1
            @test !MOI.is_valid(model, c_indices[1])
        end
    end
end


text:
module EchogramUtils

using Statistics

export db2pow, pow2db, mag2db, db2mag, hysthresh, bwselect, vertically_smooth, vertically_bin, INmask,INfilter, meandb, stddb

"""
    db2pow(ydb)

Convert decibels to power
"""
db2pow(ydb::Real) = 10^(ydb/10)

"""
    pow2db(y)

Convert power to decibels.
"""
pow2db(y::Real) = y >= 0 ? 10log10(y) : NaN

"""
    mag2db(y)

Convert magnitude to decibels.
"""
mag2db(y::Real) = 20log10(y)

"""
    db2mag(ydb)

Convert decibels to magnitude.
"""
db2mag(ydb::Real) = 10^(ydb/20)

"""
    meandb(ydb)

Remember, the Average of the Log is Not the Log of the Average,
so averaging dB requires conversion to and from the linear domain.
"""
function meandb(ydb)
    a = db2pow.(ydb)
    pow2db(mean(a))
end

function stddb(ydb)
    a = db2pow.(ydb)
    pow2db(std(a))
end

"""
    vertically_smooth(A, r, thickness)

Taking an echogram array `A` and a range vector `r`, smooth values
vertically by taking the mean over succesive `thickness` bins.

N.B. if your data is in dB, consider converting to linear (perhaps
using `dB2pow`).

"""
function vertically_smooth(A, r; thickness=1.43)
    b = zeros(size(A))
    m,n = size(A)
    for j =1:n
        for i in 0:thickness:r[end,j]
            mask = (r[:,j] .>= i)  .& (r[:,j] .< i+thickness)
            b[mask,j] .= mean(A[mask,j])
        end
    end
    return b
end

"""
    vertically_bin(A, r, thickness)

Taking an echogram array `A` and a range vector `r`, bin values
vertically by taking the mean over succesive `thickness` sized bins.

N.B. if your data is in dB, consider converting to linear (perhaps
using `dB2pow`).

"""
function vertically_bin(A, r; thickness=1.43)
    b = []
    m,n = size(A)
    for j =1:n
        for i in 0:thickness:r[end,j]
            mask = (r[:,j] .>= i)  .& (r[:,j] .< i+thickness)
            push!(b,mean(A[mask,j]))
        end
    end
    return vcat(b...)
end

"""
    INmask(A, delta)

Impulse noise filter based based on the two-sided comparison method
described by Anderson et al. (2005) and further described in Ryan et
al. (2015).

It is often desirable to linearly average to a coarser vertical
resolution (perhaps using `vertically_smooth`) before calling IN.

`A` is an echogram array and `delta` is the threshold. A sample is
marked as true if its value is more than `delta` greater than samples
on either side.

Returns a `BitArray` with the same dimensions as `A`.

"""
function INmask(A; delta=10)
    m,n = size(A)

    a = A[:,1:end-2]
    b = A[:,2:end-1]
    c = A[:,3:end]
    m1 = (b.-a) .> delta
    m2 = (b.-c) .> delta

    hcat(falses(m), m1 .& m2, falses(m))
end

function mymedian(a,b)
    (a + b) /2
end

function INfilter(A; delta=10)
    mask = INmask(A,delta=delta)

    r = A[1:end,2:end]
    l = A[1:end,1:end-1]
    x = copy(A)

    x[mask] .= mymedian.(l[mask[1:end,2:end]], r[mask[1:end,1:end-1]])
    return x
end

"""
    hysthresh(im, T1, T2)

Hysteresis thresholding

A simple port of Peter Kovesi's MATLAB method to Julia.

See http://www.peterkovesi.com/matlabfns/Spatial/hysthresh.m

Usage: bw = hysthresh(im, T1, T2)

Arguments:
             im  - image to be thresholded (assumed to be non-negative)
             T1  - upper threshold value
             T2  - lower threshold value
                   (T1 and T2 can be entered in any order, the larger of the
                   two values is used as the upper threshold)
Returns:
             bw  - the thresholded image (containing values 0 or 1)

Function performs hysteresis thresholding of an image.
All pixels with values above threshold T1 are marked as edges
All pixels that are connected to points that have been marked as edges
and with values above threshold T2 are also marked as edges.

Copyright (c) 1996-2005 Peter Kovesi
School of Computer Science & Software Engineering
The University of Western Australia
http://www.csse.uwa.edu.au/

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

"""
function hysthresh(im, T1, T2)

    if T1 < T2  # T1 and T2 reversed - swap values
        tmp = T1;
        T1 = T2;
        T2 = tmp;
    end

    # Edge points above lower threshold.
    aboveT2 = im .> T2

    # Row and colum coords of points above upper threshold.
    aboveT1r, aboveT1c = findn(im .> T1);

    # Obtain all connected regions in aboveT2 that include a point that has a
    # value above T1

    bw = bwselect(aboveT2, aboveT1c, aboveT1r);

end

"""
    bwselect(BW, c , r)

Select objects in a binary image.

Similar to the MATLAB function of the same name.
"""
function bwselect(BW, c, r)
    # constants
    north = CartesianIndex(-1,  0)
    south = CartesianIndex( 1,  0)
    east  = CartesianIndex( 0,  1)
    west  = CartesianIndex( 0, -1)

    queue = CartesianIndex.(r,c)

    m,n = size(BW)
    out = falses(m,n)

    while !isempty(queue)
        node = pop!(queue)

        if BW[node]
            wnode = node
            enode = node + east

             # Move west until node is false
            while checkbounds(Bool, BW, wnode) && BW[wnode]
                out[wnode] = true
                if checkbounds(Bool, BW, wnode + north) && !out[wnode + north]
                    push!(queue, wnode + north)
                end
                if checkbounds(Bool, BW, wnode + south) && !out[wnode + south]
                    push!(queue, wnode + south)
                end
                wnode += west
            end

            # Move east until node is false
            while checkbounds(Bool, BW, enode) && BW[enode]
                out[enode] = true
                if checkbounds(Bool, BW, enode + north) && !out[enode + north]
                    push!(queue, enode + north)
                end
                if checkbounds(Bool, BW, enode + south) && !out[enode + south]
                    push!(queue, enode + south)
                end
                enode += east
            end
        end

    end
    return out
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 1d1c0a56-7946-11eb-051f-fb214769d3f0
using PyPlot, BandedMatrices, SparseArrays, LinearAlgebra; PyPlot.svg(true);

# ╔═╡ 20390c4e-7981-11eb-16e4-a5f422480175
md"""We're implementing the method of [Oliver 1968](https://link.springer.com/article/10.1007/BF02166688), a method for solving recurrence relations of the form
```math
P_n(i) y(i) + P_{n-1}(i) y(i+1) + \cdots + P_0(i) y(i+n) = R(i)
```
where the recurrence is unstable in both the forward and backwards solutions. Oliver's method transforms an unstable forward recurrence into a boundary value problem. This can be accomplished by solving a linear system,

```math
\mathbf{P} \vec{y} = \vec{\rho}.
```
We consider the general case in which we specify ``n`` initial values, such that the terms ``Y(a+1), Y(a+2), \ldots, Y(a+n-m)`` are specified, and on the other end, ``Y(b-m+1), \ldots, Y(b)`` are specified for some integer ``m \in [0,n]``, the number of boundary values specified on the end of the recurrence. In this case, the solution is ``\vec{y}``,

```math
\vec{y} = [ Y(a + n - m + 1), \cdots , Y(b - m) ]^T.
```

"""

# ╔═╡ ca5b9a8e-7981-11eb-2dec-dfb58c50ef7d
"""
Construct a boundary value problem from a recurrence relation. The functions for coefficients and the right side dispatch on this abstract type.

`N` and `M` correspond to the total number of specified initial boundary conditions, and the number specified on the end of the recurrence, respectively.
"""
abstract type OliverProblem{T,N,M} end

# ╔═╡ c82988f2-7986-11eb-184c-4dc4685efe09
struct Example61161{T,N,M} <: OliverProblem{T,N,M} end

# ╔═╡ 1a9b22ea-7989-11eb-1123-d923b0f891f1
md"""
Consider an example recurrence relation, which we will call `Example61161`
```math
6 y(s) - 11 y(s+1) + 6 y(s+2) - y(s+3) = 0.
```
with initial conditions ``y(1) = 2``, ``y(2)=5``, and ``y(3) = 15``. This recurrence relation has an explicit solution
```math
y(s) = 1 - 2^{s-1} + 2 \cdot 3^{s-1}.
```
In this case, we actually have constant coefficients so ``P_0(i) = -1``, ``P_1(i) = 6``, ``P_2(i) = -11``, and ``P_3(i) = 6``. Further, ``R(i) = 0``.
"""

# ╔═╡ c4023f50-799d-11eb-18a4-f98202a8f420
y(s) = 1 - 2^(s-1) + 2 * 3^(s-1)

# ╔═╡ c00af4c4-7988-11eb-1dee-4f042a18d53b
"""
	OliverP(OP::OliverProblem, s, i)

Computes ``P_s(i)``, the coefficient ``s`` at solution index ``i``.
"""
function OliverP(OP::Example61161, s, i)
	if s == 0
		return -1
	elseif s == 1
		return 6
	elseif s == 2
		return -11
	elseif s == 3
		return 6
	end
end

# ╔═╡ 67cda80a-798b-11eb-020d-7de13497098c
"""
	OliverR(OP::OliverProblem, i)

Computes ``R(i)``, the right side of the recurrence relation.
"""
OliverR(OP::Example61161, i) = 0

# ╔═╡ c87eb55c-798b-11eb-05db-ab39057c8d20
"""
	assembleP{T}(OP::OliverProblem{M,N}, a, b) where {M, N}

Generic description of the ``\\mathbf{P}`` matrix within the linear system of Oliver's method.
"""
function assembleP(OP::OliverProblem{T,N,M}, a, b) where {T,N,M}
	Psize = b - N - a
	𝐏 = BandedMatrix{T}(undef, (Psize, Psize), (N-M,M))
	for row in 1:Psize
		for col in max(1, row-N+M):min(row + M, Psize)
			Δm = col - row
			𝐏[row, col] = OliverP(OP, M - Δm, a + row)
		end
	end
	return 𝐏
end

# ╔═╡ 20107d58-7981-11eb-1823-b3f867edeb85
assembleP(Example61161{Float64,3,1}(), 0, 10)

# ╔═╡ aaa88a3e-7946-11eb-3301-2720cf0f2301
"""
	assembleρ{T}(OP::OliverProblem{M,N}, a, b) where {M, N}

Generic description of the ``\\vec{\\rho}`` vector within the linear system of Oliver's method.
"""
function assembleρ(OP::OliverProblem{T,N,M}, a, b, YBC) where {T,N,M}
	Psize = b - N - a
	ρ = zeros(T, Psize)

	for row in 1:(N-M)
		Σ = zero(T)
		for s in 1:(N - M - row + 1)
			# every row, increase by 1 the sum UB, YBC arg, P arg
			Σ = Σ + OliverP(OP, N - s + 1, a + row) * YBC[a + s + row - 1]
		end
		ρ[row] = OliverR(OP, a + row) - Σ
	end

	rₛ, rₑ = (b - N - M + 1 - a), (b - N - a)
	for row in rₛ:rₑ
		row_from_end = rₑ - row + 1
		Σ = zero(T)
		for s in 0:(M - row_from_end)
			# every row, decrease by 1 the sum UB, YBC arg, P arg
			Σ = Σ + OliverP(OP, s, b - N - row_from_end + 1) *
				YBC[b - s - row_from_end + 1]
		end
		ρ[row] = OliverR(OP, a+row) - Σ
	end

	return ρ
end

# ╔═╡ f29b3902-7946-11eb-3fcc-0be7bcbf2c6c
md"""Recall the initial conditions ``y(1) = 2``, ``y(2)=5``, and ``y(3) = 15``."""

# ╔═╡ 22514dac-799a-11eb-2940-db93f451d9df
begin
	a, b, n, m = 1, 9, 3, 3
	# YBC = Dict([1,2,3] .=> [2,5,15])
	YBC = Dict( (a+1):(b+1) .=> y.((a+1):(b+1)))
	ex = Example61161{Float64,n,m}()
	P = assembleP(ex, a, b)
	ρ = assembleρ(ex, a, b, YBC)

	# """$(P * y.((a+n-m+1):(b-m))) $(ρ)"""
	sol = P \ ρ
	ref = y.((a+n-m+1):(b-m))
	"""$((sol .- ref) ./ ref)"""
end

# ╔═╡ 3e4fe168-7a07-11eb-0a47-0b31536f8773

# ╔═╡ Cell order:
# ╠═1d1c0a56-7946-11eb-051f-fb214769d3f0
# ╟─20390c4e-7981-11eb-16e4-a5f422480175
# ╠═ca5b9a8e-7981-11eb-2dec-dfb58c50ef7d
# ╠═c82988f2-7986-11eb-184c-4dc4685efe09
# ╟─1a9b22ea-7989-11eb-1123-d923b0f891f1
# ╠═c4023f50-799d-11eb-18a4-f98202a8f420
# ╠═c00af4c4-7988-11eb-1dee-4f042a18d53b
# ╠═67cda80a-798b-11eb-020d-7de13497098c
# ╠═c87eb55c-798b-11eb-05db-ab39057c8d20
# ╠═20107d58-7981-11eb-1823-b3f867edeb85
# ╠═aaa88a3e-7946-11eb-3301-2720cf0f2301
# ╟─f29b3902-7946-11eb-3fcc-0be7bcbf2c6c
# ╠═22514dac-799a-11eb-2940-db93f451d9df
# ╠═3e4fe168-7a07-11eb-0a47-0b31536f8773


text:
type ParentFunctionalCPD{D} <: CPD{D}
    target::NodeName
    parents::NodeNames
    accessor::Function # calling this gives you the distribution from the assignment

    function ParentFunctionalCPD{D}(target::NodeName, accessor::Function) where D
        new(target, NodeName[], accessor)
    end
    function ParentFunctionalCPD{D}(target::NodeName, parents::NodeNames, accessor::Function) where D
        new(target, parents, accessor)
    end
end

name(cpd::ParentFunctionalCPD) = cpd.target
parents(cpd::ParentFunctionalCPD) = cpd.parents
(cpd::ParentFunctionalCPD)(a::Assignment) = cpd.accessor(a, parents(cpd))
(cpd::ParentFunctionalCPD)() = (cpd)(Assignment()) # cpd()
(cpd::ParentFunctionalCPD)(pair::Pair{NodeName}...) = (cpd)(Assignment(pair)) # cpd(:A=>1)


text:
"""
This implements a cylindrical cutoff for the bond environments:
* A central bond (i,j) is within the cutoff if r_ij < rcutbond.
* A neighbour atom at position rkij (relative position to midpoint) is within
the environment if - after transformation to (r, θ, z) coordinates, it satisfies
`r <= rcutenv` and `abs(z) <= zcutenv`.

This struct implements the resulting filter under `env_filter`.
"""
struct CylindricalCutoff{T}
   rcutbond::T
   rcutenv::T
   zcutenv::T
end

env_filter(r, z, cutoff::CylindricalCutoff) =
      (r <= cutoff.rcutenv) && (abs(z) <= cutoff.zcutenv)

struct ACEBondPotential{TM} <: AbstractCalculator
   models::Dict{Tuple{AtomicNumber, AtomicNumber}, TM}
   cutoff::CylindricalCutoff{Float64}
end

struct ACEBondPotentialBasis{TM} <: JuLIP.MLIPs.IPBasis
   models::Dict{Tuple{AtomicNumber, AtomicNumber}, TM}  # model = basis
   inds::Dict{Tuple{AtomicNumber, AtomicNumber}, UnitRange{Int}}
   cutoff::CylindricalCutoff{Float64}
end

ACEBondCalc = Union{ACEBondPotential, ACEBondPotentialBasis}

function _get_basisinds(V::ACEBondPotential)
   inds = Dict{Tuple{AtomicNumber, AtomicNumber}, UnitRange{Int}}()
   zz = sort(collect(keys(V.models)))
   i0 = 0
   for z in zz
      mo = V.models[z]
      len = length(mo.basis)
      inds[z] = (i0+1):(i0+len)   # to generalize for general models
      i0 += len
   end
   return inds
end

_get_basisinds(V::ACEBondPotentialBasis) = V.inds

function basis(V::ACEBondPotential)
   models = Dict( [zz => model.basis for (zz, model) in V.models]... )
   inds = _get_basisinds(V)
   return ACEBondPotentialBasis(models, inds, V.cutoff)
end

# TODO:
#   - nparams
#   - get_params
#   - set_params!

Base.length(basis::ACEBondPotentialBasis) =
      sum(length(inds) for (_, inds) in basis.inds)

# --------------------------------------------------------

import JuLIP: energy, forces, virial
import ACE: evaluate, evaluate_d, grad_config

# overload the initiation of the bonds iterator to correctly extract the
# right cutoffs.
bonds(at::Atoms, calc::ACEBondCalc) =
         bonds( at, calc.cutoff.rcutbond,
                calc.cutoff.rcutbond/2 + calc.cutoff.zcutenv,
                (r, z) -> env_filter(r, z, calc.cutoff) )

_get_model(calc::ACEBondCalc, zi, zj) =
      calc.models[(min(zi, zj), max(zi,zj))]

function energy(calc::ACEBondPotential, at::Atoms)
   E = 0.0
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, calc)
      # find the right ace model
      ace = _get_model(calc, at.Z[i], at.Z[j])
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, at.Z[i], at.Z[j], Rs, Zs)
      # evaluate
      Eij = evaluate(ace, env)
      E += Eij.val
   end
   return E
end

function forces(calc::ACEBondPotential, at::Atoms)
   F = zeros(SVector{3, Float64}, length(at))
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, calc)
      Zi, Zj = at.Z[i], at.Z[j]
      # find the right ace model
      ace = _get_model(calc, Zi, Zj)
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, Zi, Zj, Rs, Zs)
      # evaluate
      dV_cyl = grad_config(ace, env)
      # transform back?
      dV_drrij, dV_dRs = rrule_eucl2cyl(rrij::SVector, Zi, Zj, Rs, Zs, dV_cyl)
      # assemble the forces
      F[i] += dV_drrij
      F[j] -= dV_drrij
      for (k, dv) in zip(Js, dV_dRs)
         F[k] -= dv
         F[i] += 0.5 * dv
         F[j] += 0.5 * dv
      end
   end
   return F
end

# site_virial(dV::AbstractVector{JVec{T1}}, R::AbstractVector{JVec{T2}}
#             ) where {T1, T2} =  (
#       length(R) > 0 ? (- sum( dVi * Ri' for (dVi, Ri) in zip(dV, R) ))
#                     : zero(JMat{fltype_intersect(T1, T2)})
#       )

function virial(calc::ACEBondPotential, at::Atoms{T}) where {T}
   V = zero(SMatrix{3, 3, T})
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, calc)
      Zi, Zj = at.Z[i], at.Z[j]
      # find the right ace model for this bond
      ace = _get_model(calc, Zi, Zj)
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, Zi, Zj, Rs, Zs)
      # evaluate
      dV_cyl = grad_config(ace, env)
      # transform back?
      dV_drrij, dV_dRs = rrule_eucl2cyl(rrij::SVector, Zi, Zj, Rs, Zs, dV_cyl)
      # assemble the virial
      #   dV_dRs contain derivative relative to midpoint
      #   dV_drrij contain derivative w.r.t. rrij
      V -= dV_drrij * rrij'
      for q = 1:length(dV_dRs)
         V -= dV_dRs[q] * Rs[q]'
      end
   end
   return V
end

function energy(basis::ACEBondPotentialBasis, at::Atoms)
   E = zeros(Float64, length(basis))
   Et = zeros(Float64, length(basis))
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, basis)
      # find the right ace model
      ace = _get_model(basis, at.Z[i], at.Z[j])
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, at.Z[i], at.Z[j], Rs, Zs)
      # evaluate
      ACE.evaluate!(Et, ace, ACE.ACEConfig(env))
      E += Et
   end
   return E
end

function forces(basis::ACEBondPotentialBasis, at::Atoms)
   F = zeros(SVector{3, Float64}, length(basis), length(at))
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, basis)
      Zi, Zj = at.Z[i], at.Z[j]
      # find the right ace model
      ace = _get_model(basis, Zi, Zj)
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, Zi, Zj, Rs, Zs)
      # evaluate
      dB_cyl = evaluate_d(ace, env)
      # transform back?
      dB_drrij, dB_dRs = rrule_eucl2cyl(rrij::SVector, Zi, Zj, Rs, Zs, dB_cyl)
      # assemble the forces
      F[:, i] += dB_drrij
      F[:, j] -= dB_drrij
      for n = 1:length(Js)
         k = Js[n]
         dv = dB_dRs[:, n]
         F[:, k] -= dv
         F[:, i] += 0.5 * dv
         F[:, j] += 0.5 * dv
      end
   end
   return F  # TODO: this is probably the wrong format
end

function virial(basis::ACEBondPotentialBasis, at::Atoms{T}) where {T}
   V = zeros(SMatrix{3, 3, T}, length(basis))
   for (i, j, rrij, Js, Rs, Zs) in bonds(at, basis)
      Zi, Zj = at.Z[i], at.Z[j]
      # find the right ace model
      ace = _get_model(basis, Zi, Zj)
      # transform the euclidean to cylindrical coordinates
      env = eucl2cyl(rrij, Zi, Zj, Rs, Zs)
      # evaluate
      dB_cyl = evaluate_d(ace, env)
      # transform back?
      dB_drrij, dB_dRs = rrule_eucl2cyl(rrij::SVector, Zi, Zj, Rs, Zs, dB_cyl)
      # assemble the virials
      for iB = 1:length(basis)
         V[iB] -= dB_drrij[iB] * rrij'
         for q = 1:length(Rs)
            V[iB] -= dB_dRs[iB, q] * Rs[q]'
         end
      end
   end
   return V  # TODO: double-check the format
end


text:
macro wandbsave(logdir)
    quote
        if !("JULIA_NO_WANDB" in keys(ENV)) || ENV["JULIA_NO_WANDB"] == "" || ENV["JULIA_NO_WANDB"] == "false" || ENV["JULIA_NO_WANDB"] == false
            wandb.save($(esc(logdir)))
        end
    end
end


text:
module MonteCarlo

using Reexport
@reexport using MonteCarloObservable
using StableDQMC, LightXML, Parameters

using Printf, SparseArrays, LinearAlgebra, Dates, Random

include("helpers.jl")
include("abstract.jl")

include("flavors/MC/MC.jl")
include("flavors/DQMC/DQMC.jl")
include("flavors/DQMC_bond/DQMC_bond.jl")

include("lattices/square.jl")
include("lattices/chain.jl")
include("lattices/cubic.jl")
include("lattices/ALPS.jl")

include("models/Ising/IsingModel.jl")
include("models/HubbardAttractive/HubbardModelAttractive.jl")
include("models/HoppingBF/HoppingBFModel.jl")

include("../test/testfunctions.jl")

export reset!
export run!
export IsingModel
export HubbardModelAttractive
export HoppingBFModel
export MC
export DQMC
export DQMC_bond
export greens
export observables

export initialize_stack #for test
export build_stack #for test
export calculate_greens#for test
export  propagate #for test
#export add_slice_sequence_left #for test
export multiply_slice_matrix_left!
end # module

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
module BouncingSphere3D_1

using Modia
using ModiaLang.StaticArrays

BouncingSphere = Model3D(
    boxHeigth = 0.1,
    groundMaterial = VisualMaterial(color="DarkGreen", transparency=0.5),
    gravField = UniformGravityField(g=9.81, n=[0, -1, 0]),
    world = Object3D(feature=Scene(gravityField=:gravField,
                                          visualizeFrames=false,
                                          defaultFrameLength=0.2,
                                          gap=0.2,
                                          enableContactDetection=true,
                                          visualizeContactPoints=false)),
    worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
    ground = Object3D(parent=:world,
                      translation=:[0.0,-boxHeigth/2,0.0],
                      feature=Solid(shape=Box(lengthX=4.0, lengthY=:boxHeigth, lengthZ=0.7),
                                    visualMaterial=:groundMaterial,
                                    solidMaterial="Steel",
                                    collision=true)),
    sphere = Object3D(feature=Solid(shape=Sphere(diameter=0.2),
                                    visualMaterial=VisualMaterial(color="Blue"),
                                    solidMaterial="Steel",
                                    massProperties=MassPropertiesFromShapeAndMass(mass=0.001),
                                    collision=true)),
    free = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=SVector{3,Float64}(0.0, 1.0, 0.0)))
)

bouncingSphere = @instantiateModel(BouncingSphere, unitless=true, log=false, logStateSelection=false, logCode=false)

#@show bouncingSphere.parameterExpressions
#@show bouncingSphere.parameters

stopTime = 2.2
dtmax = 0.1
tolerance = 1e-8
requiredFinalStates = [0.0, 0.10092547226369847, 0.0, 0.0, 0.01950941258387679, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
simulate!(bouncingSphere, stopTime=stopTime, tolerance=tolerance, dtmax=dtmax, log=true, logStates=true, logEvents=true,
          requiredFinalStates_atol=0.01, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(bouncingSphere, "free.r", figure=1)

end


text:
using JLD2
using Statistics

mutable struct IntrinsicRewardLogger <: LoggerKeyData
    intrinsic_reward::Array{Float64,1}
    log_interval::Int

    function IntrinsicRewardLogger(logger_init_info)
        num_logged_steps = logger_init_info[LoggerInitKey.TOTAL_STEPS]
        new(zeros(num_logged_steps), 1)
    end
end

function lg_step!(self::IntrinsicRewardLogger, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval)
        self.intrinsic_reward[ind] = get_reward(agent.intrinsic_reward, agent)
    end
end

function lg_episode_end!(self::IntrinsicRewardLogger, cur_step_in_episode, cur_step_total)
end

function save_log(self::IntrinsicRewardLogger, save_dict::Dict)
    save_dict[:intrinsic_reward] = self.intrinsic_reward
end

function extract_demon_weights(learner_weights,gvf_i,demon_learner)
    demon_weights = if demon_learner isa SRLearner
        immediate_rewards = learner_weights[2][gvf_i,:]
        num_SFs = demon_learner.feature_size
        ψ_weights = learner_weights[1][num_SFs*(gvf_i-1)+1:num_SFs*(gvf_i)]
        flattenall([immediate_rewards, ψ_weights])
    elseif demon_learner isa QLearner
        num_actions = demon_learner.num_actions
        gvf_w = learner_weights[(gvf_i-1)*num_actions + 1: gvf_i * num_actions,:]
        flattenall(gvf_w)
    else
        @warn "Not a valid type of demon learner for monitoring weight change contribution per demon"
    end
    return demon_weights
end

mutable struct WC_Demon_Logger <: LoggerKeyData
    wc_per_demon::Array{Float64,2}
    gvf_weights::Dict
    log_interval::Int
    function WC_Demon_Logger(logger_init_info)
        new(zeros(4,logger_init_info[LoggerInitKey.TOTAL_STEPS]), Dict(),1)
    end
end

function lg_start!(self::WC_Demon_Logger, env, agent)
    learner_weights = get_weights(agent.demon_learner)
    for gvf_i in 1:4
        gvf_weights = deepcopy(extract_demon_weights(learner_weights,gvf_i, agent.demon_learner))
        self.gvf_weights[gvf_i] = gvf_weights
    end
end

function lg_step!(self::WC_Demon_Logger, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    learner_weights = deepcopy(get_weights(agent.demon_learner))
    for gvf_i in 1:4
        gvf_weights = deepcopy(extract_demon_weights(learner_weights,gvf_i, agent.demon_learner))
        wc = sum(abs.(self.gvf_weights[gvf_i] - gvf_weights))
        self.wc_per_demon[gvf_i,cur_step_total] = wc
        self.gvf_weights[gvf_i] = gvf_weights
    end
end

function lg_episode_end!(self::WC_Demon_Logger, cur_step_in_episode, cur_step_total)
end

function save_log(self::WC_Demon_Logger, save_dict::Dict)
    save_dict[:wc_per_demon] = self.wc_per_demon
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity