text:
@testset "BufferedConfigRecorder" begin
    isfile("testfile.confs") && rm("testfile.confs")
    r = MonteCarlo.BufferedConfigRecorder{Matrix{Float64}}("testfile.confs", 2, 20)

    # Initialization
    @test r.filename == MonteCarlo.FilePath(false, "testfile.confs", "testfile.confs")
    @test length(r.buffer) == 20
    @test eltype(r.buffer) == Matrix{Float64}
    @test r.rate == 2
    @test r.idx == 1
    @test r.chunk == 1
    @test r.total_length == 0
    @test r.save_idx == -1
    @test !isdefined(r.buffer, 1)

    # Basic utility functions
    @test length(r) == r.total_length
    @test lastindex(r) == r.total_length
    @test isempty(r)

    # push to empty memory buffer
    c = rand(4, 4)
    MonteCarlo._push!(r, c)

    @test r.idx == 2
    @test r.chunk == 1
    @test r.total_length == 1
    @test r.save_idx == -1
    @test r.buffer[1] == c
    @test !isdefined(r.buffer, 2)

    @test length(r) == r.total_length
    @test lastindex(r) == r.total_length
    @test !isempty(r)
    @test r[1] == c

    # fill memory buffer - no save yet
    for _ in 2:20
        MonteCarlo._push!(r, c)
    end

    @test r.idx == 21
    @test r.chunk == 1
    @test r.total_length == 20
    @test r.save_idx == -1
    @test r.buffer[20] == c
    @test r[20] == c

    # overflow push - should save and reset buffer
    c2 = rand(4, 4)
    MonteCarlo._push!(r, c2)

    @test r.idx == 2
    @test r.chunk == 2
    @test r.total_length == 21
    @test r.save_idx == 20
    @test r.buffer[1] == c2
    @test r[21] == c2

    # query last chunk - should load and mark idx for reload
    @test r[20] == c
    @test r[1] == c
    @test r.idx == -1
    @test r.chunk == 1
    @test r.total_length == 21
    @test r.save_idx == 21

    # push to wrong chunk - should load correct chunk
    c3 = rand(4, 4)
    MonteCarlo._push!(r, c3)

    @test r.idx == 3
    @test r.chunk == 2
    @test r.total_length == 22
    @test r.save_idx == 21
    @test r.buffer[1] == c2 # this would fail if chunk not loaded
    @test r.buffer[2] == c3
    @test r.buffer[3] == c

    # test merge!
    c4 = rand(4, 4)
    cs = ConfigRecorder{Matrix{Float64}}([c4 for _ in 1:17], 10)
    merge!(r, cs)

    @test r.idx == 20
    @test r.chunk == 2
    @test r.total_length == 39
    @test r.save_idx == 21
    @test r.buffer[1] == c2   # old
    @test r.buffer[2] == c3   # old
    @test r.buffer[3] == c4   # new
    @test r.buffer[19] == c4  # new
    @test r.buffer[20] == c   # to be overwritten
    rm("testfile.confs")
end

function test_mc(mc, x)
    # Check if loaded/replayed mc matches original
    for f in fieldnames(typeof(mc.p))
        @test getfield(mc.p, f) == getfield(x.p, f)
    end
    @test mc.conf == x.conf
    @test mc.model.L == x.model.L
    @test mc.model.dims == x.model.dims
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    # @test mc.model.neighs == x.model.neighs
    @test mc.model.energy[] == x.model.energy[]
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.thermalization_measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

println("MC")
@time @testset "MC" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(
        model, beta = 0.66, thermalization = 33, sweeps = 123,
        recorder = ConfigRecorder(MC, IsingModel, 1)
    )
    @time run!(mc, verbose=false)
    save("testfile.jld2", mc)
    x = load("testfile.jld2")
    rm("testfile.jld2")
    test_mc(mc, x)

    x.measurements = MonteCarlo.default_measurements(mc, model)
    x.last_sweep = 0
    @time replay!(x, verbose=false)
    test_mc(mc, x)

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(
        model, beta=1.0, thermalization = 0, sweeps=10_000_000,
        measure_rate=10_000, recorder=ConfigRecorder(MC, IsingModel, 10_000)
    )
    @time state = run!(
        mc, verbose = false,
        safe_before = now() + Second(1),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    @info mc.last_sweep
    @time mc, state = resume!(
        "resumable_testfile.jld",
        verbose = false,
        safe_before = now() + Second(10),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld"
    )
    @info mc.last_sweep

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(
        model, beta=1.0,
        thermalization=0, sweeps=10_000length(cs),
        measure_rate=10_000, recorder=ConfigRecorder(MC, IsingModel, 10_000)
    )
    @time state = run!(mc, verbose = false)
    @test mc.configs.configs == cs.configs
    @test mc.configs.rate == cs.rate
    rm("resumable_testfile.jld")
end

function test_dqmc(mc, x)
    for f in fieldnames(typeof(mc.parameters))
        @test getfield(mc.parameters, f) == getfield(x.parameters, f)
    end
    # @test mc.conf == x.conf
    @test mc.model.mu == x.model.mu
    @test mc.model.t == x.model.t
    @test mc.model.U == x.model.U
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    @test mc.model.flv == x.model.flv
    @test mc.scheduler == x.scheduler
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                getfield(v, f).B == getfield(x.measurements[k], f).B
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            v isa MonteCarlo.DQMCMeasurement && f == :temp && continue
            v isa MonteCarlo.DQMCMeasurement && f == :kernel && continue
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                # TODO: implement ≈ for LightObservable, LogBinner, etc
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(getfield(v, f).B.compressors)
                    r = r && (a.B.compressors[i].value ≈ b.B.compressors[i].value)
                    r = r && (a.B.compressors[i].switch ≈ b.B.compressors[i].switch)
                end
                r = r && (a.B.x_sum ≈ b.B.x_sum)
                r = r && (a.B.x2_sum ≈ b.B.x2_sum)
                r = r && (a.B.count ≈ b.B.count)
            elseif getfield(v, f) isa LogBinner
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(a.compressors)
                    r = r && (a.compressors[i].value ≈ b.compressors[i].value)
                    r = r && (a.compressors[i].switch ≈ b.compressors[i].switch)
                end
                r = r && (a.x_sum ≈ b.x_sum)
                r = r && (a.x2_sum ≈ b.x2_sum)
                r = r && (a.count ≈ b.count)
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

for file in readdir()
    if endswith(file, "jld") || endswith(file, "jld2") || endswith(file, ".confs")
        rm(file)
    end
end

println("DQMC")
@time @testset "DQMC" begin
    isfile("testfile.confs") && rm("testfile.confs")
    model = HubbardModelAttractive(4, 2, t = 1.7, U = 5.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 21, sweeps = 117, measure_rate = 1,
        recorder = BufferedConfigRecorder(DQMC, HubbardModelAttractive, "testfile.confs", rate = 1)
    )
    mc[:CDC] = charge_density_correlation(mc, model)
    t = time()
    run!(mc, verbose=false)
    t = time() - t
    @info t

    save("testfile.jld", mc)
    x = load("testfile.jld")
    rm("testfile.jld")
    @test mc.conf == x.conf

    # Repeat these tests once with x being replayed rather than loaded
    test_dqmc(mc, x)

    # Check everything again with x being a replayed simulation
    x[:CDC] = charge_density_correlation(x, model)
    x.last_sweep = 0
    @time replay!(x, verbose=false)
    test_dqmc(mc, x)
    isfile("testfile.confs") && rm("testfile.confs")

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = HubbardModelAttractive(2, 2, t = 1.7, U = 5.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 0, sweeps = 10_000_000, measure_rate = 100,
        recorder = BufferedConfigRecorder(DQMC, HubbardModelAttractive, "testfile.confs", rate = 100)
    )
    mc[:CDC] = charge_density_correlation(mc, model)

    @time state = run!(
        mc, verbose = false,
        safe_before = now() + Second(2),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = [mc.recorder[i] for i in 1:length(mc.recorder)]
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    @time mc, state = resume!(
        "resumable_testfile.jld2",
        verbose = false,
        safe_before = now() + Second(10),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = [mc.recorder[i] for i in 1:length(mc.recorder)]
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld2")
    isfile("testfile.confs") && rm("testfile.confs")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = HubbardModelAttractive(2, 2, t = 1.7, U = 5.5)
    mc = DQMC(
        model, beta = 1.0, thermalization = 0, sweeps = 100length(cs), measure_rate = 100,
        recorder = BufferedConfigRecorder(DQMC, HubbardModelAttractive, "testfile.confs", rate = 100)
    )
    @time state = run!(mc, verbose = false)
    matches = true
    for i in 1:length(mc.recorder)
        matches = matches && (mc.recorder[i] == cs[i])
    end
    @test matches
    rm("resumable_testfile.jld2")
    isfile("testfile.confs") && rm("testfile.confs")
end


text:
module Hamburg

using CSV, DataFrames, Pipe
using Documenter, DocStringExtensions
export dataset

include("covid-19/Covid19.jl")
include("pollen/Pollen.jl")

"""
$(SIGNATURES)
Takes topic and name of the specific dataset and loads it into a DataFrame.

# Examples
```julia-repl
julia> dataset("covid-19", "infected")
```
"""
function dataset(topic::String, dataset::String)::DataFrame
    CSV.read(joinpath(@__DIR__, topic, "$dataset.csv"), DataFrame)
end

end # module

ext: jl


text:
function raindrops_v2(number)
    if number % 3 == 0
        if number % 5 == 0
            if number % 7 == 0
                return "PlingPlangPlong"
            else
                return "PlingPlang"
            end
        elseif number % 7 == 0
            return "PlingPlong"
        else
            return "Pling"
        end
    elseif number % 5 == 0
        if number % 7 == 0
            return "PlangPlong"
        else
            return "Plang"
        end
    elseif number % 7 == 0
        return "Plong"
    else
        return string(number)
    end
end

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
@testset "SAMIN" begin
    prob = MVP.UnconstrainedProblems.examples["Himmelblau"]

    xtrue = prob.solutions
    f = OptimTestProblems.MultivariateProblems.objective(prob)
    x0 = prob.initial_x
    res = optimize(f, x0.-100., x0.+100.0, x0, Optim.SAMIN(), Optim.Options(iterations=1000000))
    @test Optim.minimum(res) < 1e-6
end


text:
######################################################################
# CuboidMesh
######################################################################
struct CuboidMesh{D,T,D´}
    ticks::NTuple{D,Vector{T}}
    simpitr::MarchingSimplices{D,D´}
end

function Base.show(io::IO, mesh::CuboidMesh{D}) where {D}
    i = get(io, :indent, "")
    print(io,
"$(i)CuboidMesh{$D}: a mesh of a $(D)D parameter cuboid
$i  Ranges     : $(extrema.(mesh.ticks))
$i  Axes ticks : $(length.(mesh.ticks))
$i  Simplices  : $(size(mesh.simpitr)) -> $(length(mesh.simpitr))")
end

"""
    cuboid(ticks...; subticks = 13)

Create a `CuboidMesh` of L-dimensional marching-tetrahedra over a cuboid aligned with the
Cartesian axes. The dimension `L` is given by the number of `ticks`, each of the form `(x₁,
x₂,...)`. The interval between `xⱼ` and `xⱼ₊₁` ticks in axis `i` are further subdivided to
have a number of subticks including endpoints. The number is `subticks` if `subticks` is an
`Integer`, `subticks[i]` if `subticks = (s₁, s₂,...)` or `subticks[i][j]` if `subticks =
((s₁₁, s₁₂,...), (s₂₁, s₂₂,...), ...)`.

# Examples

```jldoctest
julia> cuboid((-π, π), (0, 2π); subticks = 25)
CuboidMesh{2}: a mesh of a 2D parameter cuboid
  Ranges     : ((-3.141592653589793, 3.141592653589793), (0.0, 6.283185307179586))
  Axes ticks : (25, 25)
  Simplices  : (24, 24, 2) -> 1152

julia> cuboid((-π, π), (0, 2π); subticks = (10, 10))
CuboidMesh{2}: a mesh of a 2D parameter cuboid
  Ranges     : ((-3.141592653589793, 3.141592653589793), (0.0, 6.283185307179586))
  Axes ticks : (10, 10)
  Simplices  : (9, 9, 2) -> 162
```

# External links
- Marching tetrahedra (https://en.wikipedia.org/wiki/Marching_tetrahedra) in Wikipedia
"""
cuboid(ticks::Vararg{Tuple,L}; subticks = 13) where {L} = _cuboid(sanitize_subticks(subticks, ticks), ticks...)

sanitize_subticks(st::NTuple{L,Any}, t::NTuple{L,Any}) where {L} = _sanitize_subticks.(st, t)
sanitize_subticks(st, t) = _sanitize_subticks.(Ref(st), t)
_sanitize_subticks(st::Number, t::Tuple{Vararg{Number}}) = Base.tail(_sanitize_subticks.(st, t))
_sanitize_subticks(st::Number, ::Number) = Int(st)
_sanitize_subticks(st::Tuple{Vararg{Number,L´}}, ::Tuple{Vararg{Number,L}}) where {L,L´} = L´ == L - 1 ? st :
    throw(ArgumentError("Malformed `subticks`. The number of subticks for each axis should be one less than the number of ticks for that axis"))

function _cuboid(subticks, axesticks::Vararg{Tuple,L}) where {L}
    allticks = ntuple(Val(L)) do i
        allaxisticks = typeof(1.0)[]  # We want the machine's float type, without committing to Float64
        axisticks = axesticks[i]
        nticks = length(axisticks)
        foreach(1:nticks-1) do j
            append!(allaxisticks, range(axisticks[j], axisticks[j+1], length = subticks[i][j]))
            j == nticks-1 || pop!(allaxisticks)
        end
        allaxisticks
    end
    simpitr = marchingsimplices(CartesianIndices(eachindex.(allticks)))
    return CuboidMesh(allticks, simpitr)
end

Base.eachindex(mesh::CuboidMesh) = CartesianIndices(eachindex.(mesh.ticks))

Base.getindex(mesh::CuboidMesh, n::CartesianIndex) = SVector(getindex.(mesh.ticks, Tuple(n)))
Base.getindex(mesh::CuboidMesh, i...) = mesh[eachindex(mesh)[i...]]

Base.size(mesh::CuboidMesh, i...) = size(eachindex(mesh), i...)

Base.length(mesh::CuboidMesh) = prod(length.(mesh.ticks))

vertices(mesh::CuboidMesh) = (SVector(v) for v in Iterators.product(mesh.ticks...))

nvertices(mesh::CuboidMesh) = length(vertices(mesh))

neighbors(mesh::CuboidMesh, n::CartesianIndex) = marchingneighbors(eachindex(mesh), n)

neighbors_forward(mesh::CuboidMesh, n::CartesianIndex) = marchingneighbors_forward(eachindex(mesh), n)

# function neighbors(mesh::CuboidMesh, j::Int)
#     c = eachindex(mesh)
#     l = LinearIndices(c)
#     return (l[i] for i in neighbors(mesh, c[j]))
# end

marchingsimplices(m::CuboidMesh) = m.simpitr

unitperms(m::CuboidMesh) = m.simpitr.unitperms

unitsimplex(m::CuboidMesh, i) = m.simpitr.simps[i]

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using ODEClass
using Test

@testset "ODEClass.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: chapmanbe/ODEClass

max_issues_repo_name: chapmanbe/ODEClass

max_forks_repo_name: chapmanbe/ODEClass


text:
using SpecialFunctions

@testset "special functions" begin
    # The intention here is not to check the accuracy of the libraries,
    # rather the sanity of library calls:
    piq = Float128(pi)
    halfq = Float128(0.5)
    @test gamma(halfq) ≈ sqrt(piq)
    for func in (erf, erfc, besselj0, besselj1, bessely0, bessely1, lgamma)
        @test func(halfq) ≈ func(0.5)
    end
    for func in (bessely, besselj)
        @test func(3,halfq) ≈ func(3,0.5)
    end
end


text:
#this file includes the creation of the various julian meshes, returning true if all of
#them are creates succesfully with no errors

using FEniCS

test_pytriangle = pyUnitTriangleMesh()
#test_pytetrahedron = pyUnitTetrahedronMesh() # deprecated? https://fenicsproject.org/docs/dolfin/2019.1.0/python/_autogenerated/dolfin.cpp.generation.html
test_pyquadmesh =pyUnitQuadMesh(10,10) # is one of the experimental meshes, but
#works directly as a PyObject. Investigate this at a later stage?
test_pyinterval = pyUnitIntervalMesh(10)
test_pysquare = pyUnitSquareMesh(10,10,"crossed")
test_pycube = pyUnitCubeMesh(10,10,10)
test_pybox = pyBoxMesh(Point([0.0,0.0,0.0]),Point([1.0,1.0,1.0]),10,10,10)
test_pyrectangle = pyRectangleMesh(Point([0.0, 0.0]), Point([10.0, 4.0]), 10, 10)
test_meshcopy = pyMesh("./dolphin.xml")
true


text:
abstract type Bound{T} end

eltype(::Bound{T}) where T = T

struct OpenBound{T<:Real} <: Bound{T}
    value::T
    function OpenBound{T}(x::Real) where {T<:Real}
        !(T <: Integer) || error("Bounds must be closed for integers")
        if T <: AbstractFloat
            !isnan(x) || error("Bound value must not be NaN")
            !isinf(x) || error("Bound value must not be Inf/-Inf")
        end
        new{T}(x)
    end
end
OpenBound(x::T) where T<:Real = OpenBound{T}(x)

convert(::Type{OpenBound{T}}, b::OpenBound{T}) where T<:Real = b
convert(::Type{OpenBound{T}}, b::OpenBound) where T<:Real = OpenBound{T}(b.value)
string(b::OpenBound) = string("OpenBound(", b.value, ")")

struct ClosedBound{T<:Real} <: Bound{T}
    value::T
    function ClosedBound{T}(x::Real) where {T<:Real}
        if T <: AbstractFloat
            !isnan(x) || error("Bound value must not be NaN")
            !isinf(x) || error("Bound value must not be Inf/-Inf")
        end
        new{T}(x)
    end
end
ClosedBound(x::T) where T<:Real = ClosedBound{T}(x)

convert(::Type{ClosedBound{T}}, b::ClosedBound{T}) where T<:Real = b
convert(::Type{ClosedBound{T}}, b::ClosedBound) where T<:Real = ClosedBound{T}(b.value)
string(b::ClosedBound) = string("ClosedBound(", b.value, ")")

struct NullBound{T<:Real} <: Bound{T} end
NullBound(::Type{T}) where {T<:Real} = NullBound{T}()

convert(::Type{NullBound{T}}, b::NullBound{T}) where T<:Real= b
convert(::Type{NullBound{T}}, b::NullBound) where T<:Real = NullBound{T}()
string(b::NullBound{T}) where T = string("NullBound(", T, ")")

checkvalue(a::NullBound,   x::Real) = true
checkvalue(a::OpenBound,   x::Real) = a.value <  x
checkvalue(a::ClosedBound, x::Real) = a.value <= x

checkvalue(x::Real, a::NullBound)   = true
checkvalue(x::Real, a::OpenBound)   = x <  a.value
checkvalue(x::Real, a::ClosedBound) = x <= a.value

promote_rule(::Type{Bound{T1}},::Type{Bound{T2}}) where {T1,T2} = Bound{promote_rule(T1,T2)}

function show(io::IO, b::T) where T<:Bound
    print(io, string(b))
end


text:
##==============================================================================
## FileDataEntryBlob Types
##==============================================================================
export FileDataEntry
"""
    $(TYPEDEF)
Data Entry in a file.
"""
struct FileDataEntry <: AbstractDataEntry
    label::Symbol
    id::UUID
    folder::String
    hash::String #using bytes2hex or perhaps Vector{Uint8}?
    createdTimestamp::ZonedDateTime

    function FileDataEntry(label, id, folder, hash, timestamp)
        if !isdir(folder)
            @info "Folder '$folder' doesn't exist - creating."
            # create new folder
            mkpath(folder)
        end
        return new(label, id, folder, hash, timestamp)
    end
end

# @generated function ==(x::FileDataEntry, y::FileDataEntry)
#     mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x))
# end

#
# getHash(entry::AbstractDataEntry) = hex2bytes(entry.hash)

##==============================================================================
## FileDataEntry Common
##==============================================================================
blobfilename(entry::FileDataEntry) = joinpath(entry.folder,"$(entry.id).dat")
entryfilename(entry::FileDataEntry) = joinpath(entry.folder,"$(entry.id).json")

##==============================================================================
## FileDataEntry Blob CRUD
##==============================================================================

function getDataBlob(dfg::AbstractDFG, entry::FileDataEntry)
    if isfile(blobfilename(entry))
        open(blobfilename(entry)) do f
            return read(f)
        end
    else
        error("Could not find file '$(blobfilename(entry))'.")
        # return nothing
    end
end

function addDataBlob!(dfg::AbstractDFG, entry::FileDataEntry, data::Vector{UInt8})
    if isfile(blobfilename(entry))
        error("Key '$(entry.id)' blob already exists.")
    elseif isfile(entryfilename(entry))
        error("Key '$(entry.id)' entry already exists, but no blob.")
    else
        open(blobfilename(entry), "w") do f
            write(f, data)
        end
        open(entryfilename(entry), "w") do f
            JSON.print(f, entry)
        end
        return getDataBlob(dfg, entry)::Vector{UInt8}
    end
end

function updateDataBlob!(dfg::AbstractDFG, entry::FileDataEntry, data::Vector{UInt8})
    if !isfile(blobfilename(entry))
        @warn "Entry '$(entry.id)' does not exist, adding."
        return addDataBlob!(dfg, entry, data)
    else
        # perhaps add an explicit force update flag and error otherwise
        @warn "Key '$(entry.id)' already exists, data will be overwritten."
        deleteDataBlob!(dfg, entry)
        return addDataBlob!(dfg, entry, data)
    end
end

function deleteDataBlob!(dfg::AbstractDFG, entry::FileDataEntry)
    data = getDataBlob(dfg, entry)
    rm(blobfilename(entry))
    rm(entryfilename(entry))
    return data
end

##==============================================================================
## FileDataEntry CRUD Helpers
##==============================================================================

function addData!(::Type{FileDataEntry}, dfg::AbstractDFG, label::Symbol, key::Symbol, folder::String, blob::Vector{UInt8}, timestamp=now(localzone());
                  id::UUID = uuid4(), hashfunction = sha256)
    fde = FileDataEntry(key, id, folder, bytes2hex(hashfunction(blob)), timestamp)
    de = addDataEntry!(dfg, label, fde)
    db = addDataBlob!(dfg, fde, blob)
    return de=>db
end


text:
module ESDLPlots
using ImageMagick
export plotTS, plotMAP, plotXY, plotScatter, plotMAPRGB
export plotlyjs, gadfly, gr, pyplot
using ESDL.Cubes
import ESDL.DAT
import ESDL.DAT: getFrontPerm
import ESDL.Cubes.Axes.axname
import ESDL.Cubes: findAxis, AbstractCubeData
import Reactive: Signal
import Interact: slider, dropdown, Observable, observe
import Colors: RGB, @colorant_str, colormap,  distinguishable_colors
import FixedPointNumbers: Normed
import Measures
import Compose
import Images
import DataStructures: OrderedDict
import PlotUtils: optimize_ticks, cgrad
import Compose: rectangle, text, line, compose, context, stroke, svgattribute, bitmap, HCenter, VBottom, HRight, VCenter

import Plots
import Plots: plotlyjs, gr, pyplot, plot, bar, scatter
import StatsPlots: groupedbar

const U8=Normed{UInt8,8}

abstract type ESDLPlot end
"Expression to evaluate after the data is loaded"
getafterEx(::ESDLPlot)=Expr(:block)

"Setting fixed variables"
getFixedVars(::ESDLPlot,cube)=Expr(:block)

mutable struct FixedAx
  axis
  widgetlabel::String
  musthave::Bool
  isimmu::Bool
  position::Int
end

mutable struct FixedVar
  depAxis
  varVal
  widgetlabel::String
  musthave::Bool
end

include("maps.jl")
include("other.jl")

toYr(tx::TimeAxis)=((tx.values.startyear+(tx.values.startst-1)/tx.values.NPY):(1.0/tx.values.NPY):(tx.values.stopyear+(tx.values.stopst-1)/tx.values.NPY))-(tx.values.startyear+(tx.values.startst-1)/tx.values.NPY)

r1(x)=reshape(x,length(x))
prepAx(x)=x.values
prepAx(x::TimeAxis)=toYr(x)
function repAx(x,idim,ax)
  l=length(x)
  inrep=prod(size(x)[1:idim-1])
  outrep=div(l,(inrep*size(x,idim)))
  repeat(collect(ax),inner=[inrep],outer=[outrep])
end
function count_to(f,c,i)
  ni=0
  for ind=1:i
    f(c[ind]) && (ni+=1)
  end
  return ni
end

getWidget(x::CategoricalAxis;label=axname(x))       = dropdown(x.values,label=label)
getWidget(x::RangeAxis{T};label=axname(x)) where {T<:Real} = (last(x.values)-first(x.values)) > 0 ? slider(x.values,label=label) : slider(reverse(x.values),label=label)
getWidget(x::RangeAxis;label=axname(x))             = slider(x.values,label=label)
getWidget(x::SpatialPointAxis;label="Spatial Point")= slider(1:length(x),label=label)

plotTS(x;kwargs...)=plotXY(x;xaxis="time",kwargs...)

function setPlotAxis(a::FixedAx,axlist,fixedAxes,customobs,positionobs)
  ix = a.axis===nothing ? nothing : findAxis(a.axis,axlist)
  if ix !== nothing
    push!(fixedAxes,axlist[ix])
    push!(customobs,ix)
    a.axis=ix
    positionobs[ix]=0
  else
    push!(customobs,a)
    a.axis = a.isimmu ? error("Axis $a must be selected.") : 0
  end
end
function setPlotAxis(a::FixedVar,axlist,fixedAxes,customobs,positionobs)
  a.depAxis=findAxis(a.depAxis,axlist)
  if a.varVal!=nothing
    push!(customobs,a.varVal)
    positionobs[a.depAxis]=0
  else
    push!(customobs,a)
  end
end

import Interact: observe, Widget
import InteractBase: throttle
cart(i::Integer) = CartesianIndex((i,))

function createWidgets(axlist,availableAxis,availableIndices,axlabels,widgets,axtuples,customobs,positionobs)

  if !isempty(availableAxis)
    for (icust,at) in enumerate(customobs)
      if isa(at,FixedAx)
        options = collect(at.musthave ? zip(axlabels[availableIndices],availableIndices) : zip(["None";axlabels[availableIndices]],[0;availableIndices]))
        axmenu  = dropdown(OrderedDict(options),label=at.widgetlabel,value=options[1][2],value_label=options[1][1])
        sax = observe(axmenu)
        widgets[Symbol(at.widgetlabel)]=axmenu
        customobs[icust] = sax
      elseif isa(at,FixedVar)
        w=getWidget(axlist[at.depAxis],label=at.widgetlabel)
        widgets[Symbol(at.widgetlabel)]=w
        customobs[icust] = throttle(1.0,observe(w))
      else
        #println("Skipping selected")
      end
    end
    for i in availableIndices
      w=getWidget(axlist[i])
      widgets[Symbol(axname(axlist[i]))]=w
      positionobs[i] = throttle(1.0,observe(w))
    end
  else
    for (icust,at) in enumerate(customobs)
      if isa(at,FixedAx)
        at.musthave && error("No axis left to put on $(at.widgetlabel)")
        customobs[icust] = -1
      end
    end
  end
end

const namedcolms=Dict(
:viridis=>[cgrad(:viridis)[ix] for ix in range(0,stop=1,length=100)],
:magma=>[cgrad(:magma)[ix] for ix in range(0,stop=1,length=100)],
:inferno=>[cgrad(:inferno)[ix] for ix in range(0,stop=1,length=100)],
:plasma=>[cgrad(:plasma)[ix] for ix in range(0,stop=1,length=100)])
typed_dminmax(::Type{<:Union{T,Missing}},dmin,dmax) where T = typed_dminmax(T,dmin,dmax)
typed_dminmax2(::Type{<:Union{T,Missing}},dmin,dmax) where T = typed_dminmax2(T,dmin,dmax)
typed_dminmax(::Type{T},dmin,dmax) where {T<:Integer}=(Int(dmin),Int(dmax))
typed_dminmax(::Type{T},dmin,dmax) where {T<:AbstractFloat}=(Float64(dmin),Float64(dmax))
typed_dminmax2(::Type{T},dmin,dmax) where {T<:Integer}=(isa(dmin,Tuple) ? (Int(dmin[1]),Int(dmin[2]),Int(dmin[3])) : (Int(dmin),Int(dmin),Int(dmin)), isa(dmax,Tuple) ? (Int(dmax[1]),Int(dmax[2]),Int(dmax[3])) : (Int(dmax),Int(dmax),Int(dmax)))
typed_dminmax2(::Type{T},dmin,dmax) where {T<:AbstractFloat}=(isa(dmin,Tuple) ? (Float64(dmin[1]),Float64(dmin[2]),Float64(dmin[3])) : (Float64(dmin),Float64(dmin),Float64(dmin)), isa(dmax,Tuple) ? (Float64(dmax[1]),Float64(dmax[2]),Float64(dmax[3])) : (Float64(dmax),Float64(dmax),Float64(dmax)))

mygetval(i)=i
mygetval(i::Observable)=i[]
mygetval(i::FixedAx)=-1

import Interact
import InteractBase
import Observables
import Widgets
function plotGeneric(plotObj::ESDLPlot, cube::AbstractCubeData{T};kwargs...) where T

  axlist=caxes(cube)

  axlabels=map(axname,axlist)
  fixedAxes=CubeAxis[]
  customobs=[]
  widgets=OrderedDict{Symbol,Any}()
  positionobs=Array{Any}(undef,ndims(cube))

  pAxVars=plotAxVars(plotObj)

  foreach(t->setPlotAxis(t,axlist,fixedAxes,customobs,positionobs),pAxVars)

  for (sy,val) in kwargs
    ix = findAxis(string(sy),axlist)
    if ix !== nothing
      push!(fixedAxes,axlist[ix])
      positionobs[ix] = val
    end
  end

  availableIndices=findall(ax->!in(ax,fixedAxes),axlist)
  availableAxis=axlist[availableIndices]

  createWidgets(axlist,availableAxis,availableIndices,axlabels,widgets,pAxVars,customobs,positionobs)

  allobs = [customobs;positionobs]

  if any(i->isassigned(allobs,i) && isa(allobs[i],Observable),1:length(allobs))

    s=Observable{Any}(" ")

    map!((i...)->plotCall(plotObj,cube,i...),s,customobs...,positionobs...);

    ifirstwig = findfirst(i->isassigned(allobs,i) && isa(allobs[i],Observable),1:length(allobs))
    firstwig = allobs[ifirstwig]
    firstwig[] = firstwig[]

    layout = (Widgets.manipulatelayout)((Widgets.get_backend)())
    return Widget{:ESDLPlot}(widgets,output=s,layout = layout)

  else

    return plotCall(plotObj,cube,customobs...,positionobs...)

  end
end

end # module


text:
"""
	meca(cmd0::String="", arg1=nothing; kwargs...)

Plot focal mechanisms.

Full option list at [`meca`]($(GMTdoc)meca.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)

- $(GMT.opt_B)
- **C** | **offset** :: [Type => Bool | Str | GMTcpt]

    Offsets focal mechanisms to the longitude, latitude specified in the last two columns of the input
    ($(GMTdoc)meca.html#c)
- **D** | **depth_limits** :: [Type => Str | Tuple]

    Plots events between depmin and depmax.
    ($(GMTdoc)meca.html#d)
- **E** | **fill_extensive** :: [Type => Str | Number]

    Selects filling of extensive quadrants. [Default is white].
    ($(GMTdoc)meca.html#e)
- **Fa** | **Fe** | **Fg** | **Fo** | **Fp** | **Fr** | **Ft** | **Fz** :: [Type => ]

    Sets one or more attributes.
    ($(GMTdoc)meca.html#f)
- **G** | **fill** :: [Type => Str | Number]

    Selects shade, color or pattern for filling the sectors [Default is no fill].
    ($(GMTdoc)meca.html#g)
- $(GMT.opt_P)
- **L** | **outline_pen** | **pen_outline** :: [Type => Str | Number | Tuple]

    Draws the “beach ball” outline with pen attributes instead of with the default pen set by **pen**
    ($(GMTdoc)meca.html#l)
- **M** | **same_size** | **samesize** :: [Type => Bool]

    Use the same size for any magnitude. Size is given with **S**
    ($(GMTdoc)meca.html#m)
- **N** | **no_clip** | **noclip** :: [Type => Str | []]

    Do NOT skip symbols that fall outside frame boundary.
    ($(GMTdoc)meca.html#n)
- **Sc|aki** | **Sc|CMT|gcmt** | **Sm|mt|moment_tensor** | ... :: [Type => Str]

    Selects the meaning of the columns in the input data.
    ($(GMTdoc)meca.html#s)
- **T** | **nodal** :: [Type => Number | Str]

    Plots the nodal planes and outlines the bubble which is transparent.
    ($(GMTdoc)meca.html#t)
- **W** | **pen** :: [Type => Str | Tuple]

    Set pen attributes for all lines and the outline of symbols.
    ($(GMTdoc)meca.html#w)
- **Z** | **color** | **cmap** :: [Type => Number | Str | GMTcpt]

    Give a CPT and let compressive part color be determined by the z-value in the third column.
    ($(GMTdoc)meca.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)

Example: Plot a focal mechanism using the Aki & Richards convention

```julia
    psmeca([0.0 3.0 0.0 0 45 90 5 0 0], aki=true, fill=:black, region=(-1,4,0,6), proj=:Merc, show=1)
```
The same but add a Label
```julia
    psmeca(mat2ds([0.0 3.0 0.0 0 45 90 5 0 0], "Thrust"), aki=true, fill=:black, region=(-1,4,0,6), proj=:Merc, show=1)
```
"""
function meca(cmd0::String="", arg1=nothing; first=true, kwargs...)

    proggy = (IamModern[1]) ? "meca "  : "psmeca "
	length(kwargs) == 0 && return monolitic(proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	common_mecas(cmd0, arg1, d, proggy, first, K, O)
end

# ---------------------------------------------------------------------------------------------------
meca!(cmd0::String="", arg1=nothing; kw...) = meca(cmd0, arg1; first=false, kw...)
meca(arg1,  cmd0::String=""; kw...) = meca(cmd0, arg1; first=true, kw...)
meca!(arg1, cmd0::String=""; kw...) = meca(cmd0, arg1; first=false, kw...)

const psmeca  = meca 			# Alias
const psmeca! = meca!			# Alias

"""
	coupe(cmd0::String="", arg1=nothing; kwargs...)

Plot cross-sections of focal mechanisms.

Full option list at [`pscoupe`]($(GMTdoc)coupe.html).
Essentially the same as **meca** plus **A**. Run `gmthelp(coupe)` to see the list.
"""
function coupe(cmd0::String="", arg1=nothing; first=true, kwargs...)

    proggy = (IamModern[1]) ? "coupe "  : "pscoupe "
	length(kwargs) == 0 && return monolitic(proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	common_mecas(cmd0, arg1, d, proggy, first, K, O)
end

coupe!(cmd0::String="", arg1=nothing; kw...) = coupe(cmd0, arg1; first=false, kw...)
coupe(arg1,  cmd0::String=""; kw...) = coupe(cmd0, arg1; first=true, kw...)
coupe!(arg1, cmd0::String=""; kw...) = coupe(cmd0, arg1; first=false, kw...)

const pscoupe  = coupe 			# Alias
const pscoupe! = coupe!			# Alias

# ---------------------------------------------------------------------------------------------------
function common_mecas(cmd0, arg1, d, proggy, first, K, O)
	# Helper function to both psmeca & pscoupe

	arg2 = nothing		# May be needed if GMTcpt type is sent in via C
	N_args = (arg1 === nothing) ? 0 : 1

	if (occursin("meca", proggy))
		cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12cd/0d")
		cmd, = parse_common_opts(d, cmd, [:UVXY :c :di :e :p :t :params], first)
		cmd  = parse_these_opts(cmd, d, [[:C :offset], [:D :depth_limits]])
	else
		cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX14c/10c")
		cmd, = parse_common_opts(d, cmd, [:UVXY :c :di :e :p :t :params], first)
		cmd_ = add_opt(d, "", "Aa", [:Aa :cross_ll_pts], (lon1="", lat1="", lon2="", lat2="", dip="", width="", dmin="", dmax="", frame="_+f"))
		if (cmd_ == "")
			cmd_ = add_opt(d, "", "Ab", [:Ab :cross_ll_azim], (lon1="", lat1="", strike="", length="", dip="", width="", dmin="", dmax="", frame="_+f"))
		end
		if (cmd_ == "")
			cmd_ = add_opt(d, "", "Ac", [:Ac :cross_xy_pts], (x1="", y1="", x2="", y2="", dip="", width="", dmin="", dmax="", frame="_+f"))
		end
		if (cmd_ == "")
			cmd_ = add_opt(d, "", "Ad", [:Ad :cross_xy_azim], (x1="", y1="", strike="", length="", dip="", width="", dmin="", dmax="", frame="_+f"))
		end
		(cmd_ == "" && !show_kwargs[1]) && error("Specifying cross-section type is mandatory")
		cmd *= cmd_
	end

	cmd = add_opt_fill(cmd, d, [:E :fill_extensive], 'E')
	cmd  = parse_these_opts(cmd, d, [[:L :outline_pen :pen_outline], [:M :same_size :samesize],
	                                 [:N :no_clip :noclip], [:T :nodal]])
	cmd  = parse_these_opts(cmd, d, [[:Fa :PT_axes], [:Fe :T_axis_color], [:Fg :P_axis_color], [:Fo :psvelo],
	                                 [:Fp :P_axis_pen], [:Fr :label_box], [:Ft :T_axis_pen], [:Fz :zero_trace]])
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	#(occursin("coupe", proggy)) && (cmd = add_opt(d, cmd, "Q", [:Q]))

	# If file name sent in, read it and compute a tight -R if it was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)
	n_cols = (isa(arg1, Matrix) || isa(arg1, GMTdataset)) ? size(arg1,2) : 0

	if     (haskey(d, :Sa) || haskey(d, :aki))  symbs = [:Sa :aki]
	elseif (haskey(d, :Sc) || haskey(d, :CMT) || haskey(d, :gcmt))  symbs = [:Sc :CMT :gcmt]
	elseif (haskey(d, :Sm) || haskey(d, :mt) || haskey(d, :moment_tensor))  symbs = [:Sm :mt :moment_tensor]
	elseif (haskey(d, :Sd) || haskey(d, :mt_closest) || haskey(d, :moment_closest))  symbs = [:Sd :mt_closest :moment_closest]
	elseif (haskey(d, :Sz) || haskey(d, :mt_deviatoric) || haskey(d, :moment_deviatoric))  symbs = [:Sz :mt_deviatoric :moment_deviatoric]
	elseif (haskey(d, :Sp) || haskey(d, :partial))  symbs = [:Sp :partial]
	elseif (haskey(d, :Sx) || haskey(d, :principal) || haskey(d, :principal_axis))  symbs = [:Sx :principal :principal_axis]
	elseif (haskey(d, :Sy) || haskey(d, :principal_closest))     symbs = [:Sy :principal_closest]
	elseif (haskey(d, :St) || haskey(d, :principal_deviatoric))  symbs = [:St :principal_deviatoric]
	elseif (n_cols >= 7 && n_cols < 10)  d[:Sa] = true; symbs = [:Sa]		# Implicit Aki
	elseif (n_cols >= 11 && n_cols < 14) d[:Sc] = true; symbs = [:Sc]		# Implicit CMT
	elseif (show_kwargs[1])  symbs = [:Sa :aki :Sc :CMT :gcmt :Sm :mt :Sd :mt_closest :moment_closest :Sz :mt_deviatoric :moment_deviatoric :Sp :partial :Sx :principal :principal_axis :Sy :principal_closest :St :principal_deviatoric]
	else  error("Must select one convention")
	end
	cmd_ = add_opt(d, "", string(symbs[1]), symbs, (scale="", angle="+a", font=("+f", font), justify="+j", offset="+o"))
	if (length(cmd_) != 0 && (length(cmd_) == 4 || cmd_[5] == '+'))		# If scale not given search for the 'scale' kwarg
		cmd = ((val = find_in_dict(d, [:scale])[1]) !== nothing) ? cmd * string(cmd_,val) : cmd * cmd_ * "2.5c"
	end

	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:Z :C :color :cmap], 'Z', N_args, arg1, arg2)
	cmd *= opt_pen(d, 'W', [:W :pen])

	finish_PS_module(d, proggy * cmd, "", K, O, true, arg1, arg2)
end


text:
module MultiModalMuSig

using DataFrames
using LinearAlgebra
using Statistics
using SpecialFunctions
using NLopt

export IMMCTM, MMCTM, ILDA, LDA, fit!, format_counts_lda, format_counts_ctm, format_counts_mmctm

include("IMMCTM.jl")
include("MMCTM.jl")
include("ILDA.jl")
include("LDA.jl")
include("common.jl")
include("utils.jl")

end # module


text:
using Turing
using Turing.RandomMeasures
using Random
using AutoGibbs
using BenchmarkTools
using Plots
using StatsPlots

include("models.jl")
include("models_comparison.jl")

function benchmark_models(combinations...; N=100, bparams...)
    old_prog = Turing.PROGRESS[]
    Turing.turnprogress(false)
    suite = BenchmarkGroup()

    benchmarks = map(combinations) do (name, (m, alg))
        suite[name] = @benchmarkable sample($m, $alg, $N)
    end

    tune!(suite)
    results = run(suite;bparams...)

    Turing.turnprogress(old_prog)
    return results
end

function compare_gmm()
    data = [randn(10) .- 0.5; randn(5) .+ 1]
    m1 = gmm_example(data)
    m2 = gmm_tarray_example(data)

    cond_gmm_z = StaticConditional(m1, :z)
    alg1 = Gibbs(cond_gmm_z, MH(:w, :μ))
    alg2 = Gibbs(cond_gmm_z, HMC(0.01, 10, :w, :μ))
    alg3 = Gibbs(PG(20, :z), MH(:w, :μ))
    alg4 = Gibbs(PG(20, :z), HMC(0.01, 10, :w, :μ))

    b = benchmark_models("AG+MH" => (m1, alg1),
                         "AG+HMC" => (m1, alg2),
                         "PG+MH" => (m2, alg3),
                         "PG+HMC" => (m2, alg4);
                         N=1000, verbose=true, samples=10)
    display(b)
end

compare_gmm()


text:
extractionfn = PartialEF("partial")
@testset "DefaultDS" begin
    @test_throws ErrorException("Invalid outputType") DefaultDS("dimension", outputName="name", outputType="type")
    @test_throws TypeError DefaultDS("dimension", outputName="name", outputType=1)
    @test_throws TypeError DefaultDS("dimension", outputName=0, outputType="string")
end

@testset "ExtractionDS" begin
    @test_throws ErrorException("Invalid outputType") ExtractionDS("dimension", PartialEF("partial"), outputName="name", outputType="type")
    @test_throws TypeError ExtractionDS("dimension", extractionfn, outputName="name", outputType=1)
    @test_throws TypeError ExtractionDS("dimension", extractionfn, outputName=0, outputType="string")
end

dimspec = DefaultDS("dimension")
@testset "ListFiltered" begin
    @test_throws TypeError ListFiltered(dimspec, ["a", "b"], isWhitelist=1)
end

@testset "LookupDS" begin
    @test_throws TypeError LookupDS("dimension", "name", outputName=1)
end

lookupmap = Map(true, Dict())
@testset "MapLookupDS" begin
    @test_throws TypeError MapLookupDS("dimension", lookupmap, outputName=1)
    @test_throws TypeError MapLookupDS("dimension", lookupmap, retainMissingValue=1)
    @test_throws TypeError MapLookupDS("dimension", lookupmap, replaceMissingValueWith=1)
    @test_throws TypeError MapLookupDS("dimension", lookupmap, optimize=1)
    @test_throws ErrorException("Cannon specify replaceMissingValueWith when retainMissingValue == true") MapLookupDS("dimension", lookupmap, retainMissingValue=true, replaceMissingValueWith="a")
end

@testset "RegexEF" begin
    @test_throws ErrorException("index must be a non-negative integer") RegexEF("a", index=-1)
    @test_throws ErrorException("index must be a non-negative integer") RegexEF("a", index="a")
    @test_throws TypeError RegexEF("a", replaceMissingValue=1)
end

@testset "SubstringEF" begin
    @test_throws ErrorException("length must be a non-negative integer") SubstringEF(1, length=-1)
    @test_throws ErrorException("length must be a non-negative integer") SubstringEF(1, length="a")
end

@testset "TimeFormatEF" begin
    @test_throws TypeError TimeFormatEF(format=1)
    @test_throws TypeError TimeFormatEF(timeZone=1)
    @test_throws TypeError TimeFormatEF(locale=1)
    @test_throws TypeError TimeFormatEF(granularity=1)
    @test_throws TypeError TimeFormatEF(asMillis=1)
end

@testset "TimeParseEF" begin
    @test_throws TypeError TimeParseEF("format", "format", joda=1)
end

@testset "JavaScriptEF" begin
    @test_throws TypeError(:typeassert, "", Bool, 1) JavaScriptEF("function", injective=1)
end

@testset "RegisteredLookupEF" begin
    @test_throws TypeError(:typeassert, "", Bool, 1) RegisteredLookupEF("lookup", injective=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) RegisteredLookupEF("lookup", retainMissingValue=1)
    @test_throws TypeError(:typeassert, "", String, 1) RegisteredLookupEF("lookup", replaceMissingValueWith=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) RegisteredLookupEF("lookup", optimize=1)
    @test_throws ErrorException("Cannon specify replaceMissingValueWith when retainMissingValue == true") RegisteredLookupEF("lookup", retainMissingValue=true, replaceMissingValueWith="a")
end

@testset "InlineLookupEF" begin
    @test_throws TypeError(:typeassert, "", Bool, 1) InlineLookupEF(lookupmap, injective=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) InlineLookupEF(lookupmap, retainMissingValue=1)
    @test_throws TypeError(:typeassert, "", String, 1) InlineLookupEF(lookupmap, replaceMissingValueWith=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) InlineLookupEF(lookupmap, optimize=1)
    @test_throws ErrorException("Cannon specify replaceMissingValueWith when retainMissingValue == true") InlineLookupEF(lookupmap, retainMissingValue=true, replaceMissingValueWith="a")
end

@testset "StringFormatEF" begin
    @test_throws ErrorException("Invalid nullHandling") StringFormatEF("format", nullHandling=1)
    @test_throws ErrorException("Invalid nullHandling") StringFormatEF("format", nullHandling="a")
end

@testset "UpperEF" begin
    @test_throws TypeError(:typeassert, "", String, 1) UpperEF(locale=1)
end

@testset "LowerEF" begin
    @test_throws TypeError(:typeassert, "", String, 1) LowerEF(locale=1)
end

@testset "BucketEF" begin
    @test_throws ErrorException("size must be a non-negative integer") BucketEF(size="a")
    @test_throws ErrorException("size must be a non-negative integer") BucketEF(size=-1)
    @test_throws ErrorException("offset must be a non-negative integer") BucketEF(offset="a")
    @test_throws ErrorException("offset must be a non-negative integer") BucketEF(offset=-1)
end


text:
# This is model 13.4 from Statistical Rethinking Edition 1.

# \toc

# ## Data

import CSV

using DataFrames
using Turing
using TuringModels

data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv")
df = CSV.read(data_path, DataFrame; delim=';')

dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept)))
df.male = [g == "male" ? 1 : 0 for g in df.gender]
df.dept_id = [dept_map[de] for de in df.dept]
df

# ## Model

@model function m13_4(applications, dept_id, male, admit)
    sigma_dept ~ truncated(Cauchy(0, 2), 0, Inf)
    a ~ Normal(0, 10)
    a_dept ~ filldist(Normal(a, sigma_dept), 6)

    logit_p = a_dept[dept_id]

    admit .~ BinomialLogit.(applications, logit_p)
end;

# ## Output

chains = sample(
    m13_4(df.applications, df.dept_id, df.male, df.admit),
    Turing.NUTS(0.65),
    5000
)

# \defaultoutput{}

# ## Original output

m_13_4_rethinking = """
Inference for Stan model: 0c5c36512c20c41995d1bd25be525138.
    3 chains, each with iter=4500; warmup=500; thin=1;
    post-warmup draws per chain=4000, total post-warmup draws=12000.

                   mean se_mean   sd     2.5%      25%      50%      75%    97.5%
    a_dept[1]      0.59    0.00 0.07     0.46     0.54     0.59     0.64     0.73
    a_dept[2]      0.54    0.00 0.09     0.37     0.48     0.54     0.60     0.71
    a_dept[3]     -0.62    0.00 0.07    -0.75    -0.66    -0.62    -0.57    -0.49
    a_dept[4]     -0.67    0.00 0.08    -0.81    -0.72    -0.67    -0.62    -0.52
    a_dept[5]     -1.09    0.00 0.10    -1.28    -1.15    -1.09    -1.03    -0.90
    a_dept[6]     -2.65    0.00 0.15    -2.96    -2.75    -2.65    -2.55    -2.37
    a             -0.63    0.01 0.63    -1.89    -1.01    -0.64    -0.26     0.66
    sigma_dept     1.45    0.01 0.57     0.77     1.08     1.33     1.68     2.90
    lp__       -2602.37    0.03 2.08 -2607.45 -2603.49 -2602.01 -2600.86 -2599.35
               n_eff Rhat
    a_dept[1]  11285    1
    a_dept[2]  11943    1
    a_dept[3]  12707    1
    a_dept[4]  11592    1
    a_dept[5]  14598    1
    a_dept[6]  13363    1
    a           6965    1
    sigma_dept  5643    1
    lp__        5322    1
""";


text:
module ROSTransforms

    # tf functions
    include("tf.jl")
    # tf2_ros functions
    include("tf2_ros.jl")

end


text:
# Unit tests for analysis/plotting.jl

@testset "analysis/plotting.jl" begin
    # TODO
end


text:
# -------------------------------------------------------
#                    Abstract types
# -------------------------------------------------------
abstract type Checkerboard end
abstract type CBTrue <: Checkerboard end
abstract type CBFalse <: Checkerboard end
abstract type CBGeneric <: CBTrue end
abstract type CBAssaad <: CBTrue end

abstract type AbstractDQMC{C<:Checkerboard, GreensEltype<:Number, HoppingEltype<:Number} end

(@isdefined DQMC_CBTrue) || (global const DQMC_CBTrue = AbstractDQMC{C} where C<:CBTrue)
(@isdefined DQMC_CBFalse) || (global const DQMC_CBFalse = AbstractDQMC{C} where C<:CBFalse)

# -------------------------------------------------------
#                    Includes
# -------------------------------------------------------
using MonteCarloObservable
using TimerOutputs
using FFTW
using HDF5
using JLD
using LightXML

# using Iterators
using Dates
using LinearAlgebra
using SparseArrays
using Statistics
using Printf
using Random

(@isdefined TIMING) || (global const TIMING = false)
macro mytimeit(exprs...)
    if TIMING
        return :(@timeit($(esc.(exprs)...)))
    else
        return esc(exprs[end])
    end
end

"""
Macro to define substitute macros.
Concretely, `@def mymacro body` will define a macro
such that `@mymacro` will be replaced by `body`.
"""
macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end
export @def

include("helpers.jl")
include("parameters.jl")
include("lattice.jl")
include("stack.jl")
include("slice_matrices.jl")
include("linalg.jl")
include("statistics.jl") # here for convenience
include("hoppings.jl")
include("hoppings_checkerboard.jl")
include("hoppings_checkerboard_generic.jl")
include("interactions.jl")
include("action.jl")
include("local_updates.jl")
include("global_updates.jl")
# include("observable.jl")
include("boson_measurements.jl")
include("fermion_measurements.jl")

# tools for apps
include("tools/xml_parameters.jl")
include("tools/hdf5_parameters.jl")

# -------------------------------------------------------
#                    Concrete types
# -------------------------------------------------------
mutable struct Analysis
  acc_rate::Float64
  acc_rate_global::Float64
  prop_global::Int
  acc_global::Int
  to::TimerOutputs.TimerOutput

  Analysis() = new(0.0,0.0,0,0,TimerOutput())
end

mutable struct DQMC{C<:Checkerboard, GreensEltype<:Number, HoppingEltype<:Number} <: AbstractDQMC{C, GreensEltype, HoppingEltype}
  p::Params
  l::Lattice{HoppingEltype}
  s::Stack{GreensEltype}
  a::Analysis
end

DQMC(p::Params; force_generic_chkr=false) = begin
  CB = CBFalse

  # choose checkerboard variant
  if p.chkr
    if p.Nhoppings == "none" && p.NNhoppings == "none"
      CB = iseven(p.L) && !force_generic_chkr ? CBAssaad : CBGeneric
    else
      CB = CBGeneric
    end
  end
  println()
  @show CB
  @show p.hoppings
  @show p.Nhoppings
  @show p.NNhoppings
  @show p.lattice_file
  println()

  ### SET DATATYPES
  G = ComplexF64
  H = ComplexF64
  if !p.Bfield
    H = Float64;
    G = p.opdim > 1 ? ComplexF64 : Float64; # O(1) -> real GF
  end

  mc = DQMC{CB,G,H}(p, Lattice{H}(), Stack{G}(), Analysis())
  load_lattice(mc)
  mc
end

# type helpers
@inline geltype(mc::DQMC{CB, G, H}) where {CB, G, H} = G
@inline heltype(mc::DQMC{CB, G, H}) where {CB, G, H} = H
@inline cbtype(mc::DQMC{CB, G, H}) where {CB, G, H} = CB

"""
Memory usage of the `AbstractDQMC` in MB.
"""
memory_usage(mc::AbstractDQMC) = round(Base.summarysize(mc) / 1024 / 1024, digits=1)

# -------------------------------------------------------
#                  Monte Carlo
# -------------------------------------------------------
function init!(mc::DQMC)
  Random.seed!(mc.p.seed); # init RNG
  init!(mc, rand(mc.p.opdim,mc.l.sites,mc.p.slices), false)
end
function init!(mc::DQMC, start_conf, init_seed=true)
  a = mc.a
  @mytimeit a.to "init mc" begin
  init_seed && Random.seed!(mc.p.seed); # init RNG

  # Init hsfield
  println("\nInitializing HS field")
  mc.p.hsfield = start_conf
  println("Initializing boson action\n")
  mc.p.boson_action = calc_boson_action(mc)

  # stack init and test
  initialize_stack(mc)
  println("Building stack")
  build_stack(mc)
  println("Initial propagate: ", mc.s.current_slice, " ", mc.s.direction)
  propagate(mc)
  end

  TIMING && show(TimerOutputs.flatten(a.to); allocations = true)
  nothing
end

function run!(mc::DQMC)
  println("\n\nMC Thermalize - ", mc.p.thermalization*2)
  flush(stdout)
  thermalize!(mc)

  h5open(mc.p.output_file, "r+") do f
    HDF5.has(f, "resume/box") && o_delete(f, "resume/box")
    HDF5.has(f, "resume/box_global") && o_delete(f, "resume/box_global")
    write(f, "resume/box", mc.p.box)
    write(f, "resume/box_global", mc.p.box_global)
  end

  println("\n\nMC Measure - ", mc.p.measurements*2)
  flush(stdout)
  measure!(mc)
  nothing
end

function resume!(mc::DQMC, lastconf, prevmeasurements::Int)
  p = mc.p

  # Init hsfield
  println("\nLoading last HS field")
  p.hsfield = copy(lastconf)
  println("Initializing boson action\n")
  p.boson_action = calc_boson_action(mc)

  jldopen(p.output_file, "r") do fjld
    f = fjld.plain
    box = read(f, "resume/box")
    box_global = read(f, "resume/box_global")
    p.box = box
    p.box_global = box_global

    if HDF5.has(f, "thermal_init/udsweep")
      try
        mc.a.to = read(fjld, "thermal_init/udsweep")
      catch err
        @warn "Resuming measurements but couldn't load \"thermal_init/udsweep\" timer."
      end
    else
      @warn "Resuming measurements but couldn't find \"thermal_init/udsweep\" timer."
    end
  end

  println("\n\nMC Measure (resuming) - ", p.measurements*2, " (total $((p.measurements + prevmeasurements)*2))")
  flush(stdout)
  measure!(mc; prevmeasurements = prevmeasurements)

  nothing
end

function thermalize!(mc::DQMC)
  a = mc.a
  p = mc.p

  a.acc_rate = 0.0
  a.acc_rate_global = 0.0
  a.prop_global = 0
  a.acc_global = 0

  reset_timer!(a.to)
  if p.prethermalized > 0
    # load old a.to["udsweep"] timer
    jldopen(p.output_file, "r") do f
      a.to = read(f["thermal_init/udsweep"])
    end
  end

  for i in (p.prethermalized+1):p.thermalization
    udswdur = @elapsed @timeit a.to "udsweep" for u in 1:2 * p.slices
      update(mc, i)
    end
    @printf("\tsweep duration: %.4fs\n", udswdur/2)
    flush(stdout)

    if mod(i, 10) == 0
      a.acc_rate = a.acc_rate / (10 * 2 * p.slices)
      a.acc_rate_global = a.acc_rate_global / (10 / p.global_rate)
      println("\n\t", i*2)
      @printf("\t\tsweep dur (total mean): %.4fs\n", TimerOutputs.time(a.to["udsweep"])/2 *10.0^(-9)/TimerOutputs.ncalls(a.to["udsweep"]))
      @printf("\t\tacc rate (local) : %.1f%%\n", a.acc_rate*100)
      if p.global_updates
        @printf("\t\tacc rate (global): %.1f%%\n", a.acc_rate_global*100)
        @printf("\t\tacc rate (global, overall): %.1f%%\n", a.acc_global/a.prop_global*100)
      end

      # adaption (only during thermalization)
      if a.acc_rate < 0.5
        @printf("\t\tshrinking box: %.2f\n", 0.9*p.box)
        p.box = 0.9 * p.box
      else
        @printf("\t\tenlarging box: %.2f\n", 1.1*p.box)
        p.box = 1.1 * p.box
      end

      if p.global_updates
        if a.acc_global/a.prop_global < 0.5
          @printf("\t\tshrinking box_global: %.2f\n", 0.9*p.box_global)
          p.box_global = 0.9 * p.box_global
        else
          @printf("\t\tenlarging box_global: %.2f\n", 1.1*p.box_global)
          p.box_global = 1.1 * p.box_global
        end
      end
      a.acc_rate = 0.0
      a.acc_rate_global = 0.0
      println()
      flush(stdout)
    end

    # Save thermal configuration for "resume"
    if (shutdown = now() >= p.walltimelimit) || i%p.write_every_nth == 0
      _save_thermal(mc, i)
      if shutdown
        println("Approaching wall-time limit. Safely exiting.")
        exit(42)
      end
    end

  end

  if TIMING
    # display(a.to);
    show(TimerOutputs.flatten(a.to); allocations = true)
    # save(p.output_file[1:end-10]*"timings.jld", "to", a.to); # TODO
    rm(p.output_file)
    exit();
  end
  nothing
end

_save_thermal(mc, i) = begin
    jldopen(mc.p.output_file, "r+") do fjld
      f = fjld.plain
      HDF5.has(f, "thermal_init/conf") && HDF5.o_delete(f, "thermal_init/conf")
      HDF5.has(f, "thermal_init/prethermalized") && HDF5.o_delete(f, "thermal_init/prethermalized")
      HDF5.has(f, "thermal_init/udsweep") && HDF5.o_delete(f, "thermal_init/udsweep")
      write(fjld, "thermal_init/conf", mc.p.hsfield)
      write(fjld, "thermal_init/prethermalized", i)
      write(fjld, "thermal_init/udsweep", mc.a.to["udsweep"])
      (i == mc.p.thermalization) && saverng(f; group="thermal_init/rng") # for future features
    end
end

"""
Create Observable objects and store them in a NamedTuple
"""
function _create_observables(mc::AbstractDQMC, cs)
  observables = mc.p.obs
  outfile = mc.p.output_file
  obs = NamedTuple()

  if :greens in observables
    g = Observable(typeof(mc.s.greens), "greens"; alloc=cs,
                                                  inmemory=false,
                                                  outfile=outfile,
                                                  group="obs/greens")
    obs = add(obs, greens = g)
  end

  if :boson_action in observables
    boson_action = Observable(Float64, "boson_action"; alloc=cs,
                                                       inmemory=false,
                                                       outfile=outfile,
                                                       group="obs/boson_action")
    obs = add(obs, boson_action = boson_action)
  end

  return obs
end

"""
Measure observables like the equal-time Green's function if requested
"""
function _measure_observables!(mc, obs)
        # bosonic
        if haskey(obs, :boson_action)
          push!(obs[:boson_action], mc.p.boson_action)
        end

        # fermionic
        if haskey(obs, :greens)
          g = wrap_greens(mc,mc.s.greens,mc.s.current_slice,1)
          effective_greens2greens!(mc, g)
          push!(obs[:greens], g)
        end
        nothing
end

function measure!(mc::DQMC; prevmeasurements = 0)
  a = mc.a
  l = mc.l
  p = mc.p
  s = mc.s

  initialize_stack(mc)
  println("Renewing stack")
  build_stack(mc)
  println("Initial propagate: ", s.current_slice, " ", s.direction)
  propagate(mc)

  cs = choose_chunk_size(mc)

  configurations = Observable(typeof(p.hsfield), "configurations"; alloc=cs, inmemory=false, outfile=p.output_file, group="obs/configurations")

  obs = _create_observables(mc, cs)

  i_start = 1
  i_end = p.measurements

  if p.resume
    restorerng(p.output_file; group="resume/rng")
    togo = mod1(prevmeasurements, p.write_every_nth)-1
    i_start = prevmeasurements-togo+1
    i_end = p.measurements + prevmeasurements
    cs = configurations.alloc
    println("Overriding cs to match cs of previous run (cs = $(cs)).")
  end

  @show cs

  acc_rate = 0.0
  acc_rate_global = 0.0
  for i in i_start:i_end
    @timeit a.to "udsweep" for u in 1:2 * p.slices
      update(mc, i)

      # if s.current_slice == 1 && s.direction == 1 && (i-1)%p.write_every_nth == 0 # measure criterium
      if s.current_slice == p.slices && s.direction == -1 && (i-1)%p.write_every_nth == 0 # measure criterium
        dumping = (length(configurations)+1)%cs == 0
        dumping && println("Dumping... (2*i=$(2*i))")
        # push!(boson_action, p.boson_action)

        push!(configurations, p.hsfield)

        _measure_observables!(mc, obs)

        dumping && saverng(p.output_file; group="resume/rng")
        dumping && println("Dumping block of $cs datapoints was a success")
        flush(stdout)

        if !dumping # measuring but we aren't dumping
          # Estimate whether we'll make it to another measurement before hitting WTL. If not, flush and restart.
          to_udsweep = mc.a.to["udsweep"]
          udsd = TimerOutputs.time(to_udsweep) *10.0^(-9)/TimerOutputs.ncalls(to_udsweep)
          secs_to_meas = p.write_every_nth * udsd
          secs_to_meas *= 1.1 # add 10 percent because we might be slower
          next_meas_date = now() + Millisecond(ceil(Int, secs_to_meas*1000))

          if next_meas_date >= p.walltimelimit
            println("Approaching wall-time limit. Won't make it to next measurement. Safely exiting. (i = $(i)). Current date: $(Dates.format(now(), "d.u yyyy HH:MM")).")
            println("Flushing configurations which haven't been dumped yet.")
            flush(configurations)
            for o in obs
              flush(o)
            end
            # flush(boson_action)
            saverng(p.output_file; group="resume/rng")
            exit(42)
          end
        end
      end
    end

    if mod(i, 100) == 0
      a.acc_rate = a.acc_rate / (100 * 2 * p.slices)
      a.acc_rate_global = a.acc_rate_global / (100 / p.global_rate)
      println("\t", i*2)
      @printf("\t\tsweep dur (total mean): %.4fs\n", TimerOutputs.time(a.to["udsweep"])/2 *10.0^(-9)/TimerOutputs.ncalls(a.to["udsweep"]))
      @printf("\t\tacc rate (local) : %.1f%%\n", a.acc_rate*100)
      if p.global_updates
        @printf("\t\tacc rate (global): %.1f%%\n", a.acc_rate_global*100)
        @printf("\t\tacc rate (global, overall): %.1f%%\n", a.acc_global/a.prop_global*100)
      end
      a.acc_rate = 0.0
      a.acc_rate_global = 0.0
      flush(stdout)
    end

    if (approaching_wtl = now() >= p.walltimelimit)
      println("Approaching wall-time limit. Safely exiting immediately. (i = $(i)). Current date: $(Dates.format(now(), "d.u yyyy HH:MM")).")
      println("Won't flush because state wouldn't be resumable.")
      exit(42)
    end
  end

  println("Final flush.")
  flush(configurations)
  for o in obs
    flush(o)
  end

  nothing
end

function update(mc::DQMC, i::Int)
  p = mc.p
  s = mc.s
  l = mc.l
  a = mc.a

  propagate(mc)

  if p.global_updates && (s.current_slice == p.slices && s.direction == -1 && mod(i, p.global_rate) == 0)
    a.prop_global += 1
    @mytimeit a.to "global updates" b = global_update(mc)
    a.acc_rate_global += b
    a.acc_global += b
  end

  a.acc_rate += local_updates(mc)
  nothing
end

# -------------------------------------------------------
#                     Other stuff
# -------------------------------------------------------
# cosmetics
import Base.summary
import Base.show
Base.summary(mc::DQMC) = "DQMC"
function Base.show(io::IO, mc::DQMC{C}) where C<:Checkerboard
  print(io, "DQMC of O($(mc.p.opdim)) model\n")
  print(io, "r = ", mc.p.r, ", λ = ", mc.p.lambda, ", c = ", mc.p.c, ", u = ", mc.p.u, "\n")
  print(io, "Beta: ", mc.p.beta, " (T ≈ $(round(1/mc.p.beta, digits=3)))", "\n")
  print(io, "L: ", mc.p.L, " (N = ", mc.l.sites , ")\n")
  print(io, "Checkerboard: ", C, "\n")
  print(io, "B-field: ", mc.p.Bfield)
end
Base.show(io::IO, m::MIME"text/plain", mc::DQMC) = print(io, mc)

"""
Calculate DateTime where wall-time limit will be reached.

Example call: wtl2DateTime("3-12:42:05", now())
"""
function wtl2DateTime(wts::AbstractString, start_time::DateTime)
  @assert occursin("-", wts)
  @assert occursin(":", wts)
  @assert length(wts) >= 10

  tmp = split(wts, "-")

  d = parse(Int, tmp[1])
  h, m, s = parse.(Int, split(tmp[2], ":"))

  start_time + Dates.Day(d) + Dates.Hour(h) + Dates.Minute(m) + Dates.Second(s)
end

function set_walltimelimit!(p, start_time)
  if "WALLTIMELIMIT" in keys(ENV)
    p.walltimelimit = wtl2DateTime(ENV["WALLTIMELIMIT"], start_time)
    @show ENV["WALLTIMELIMIT"]
    @show p.walltimelimit
  elseif occursin("cheops", gethostname())
    p.walltimelimit = wtl2DateTime("9-23:30:00", start_time) # CHEOPS
    println("Set CHEOPS walltime limit, i.e. 9-23:30:00.")
  elseif occursin("jw", gethostname())
    p.walltimelimit = wtl2DateTime("0-23:30:00", start_time) # JUWELS, was "0-23:30:00"
    println("Set JUWELS walltime limit, i.e. 0-23:30:00.")
  end

  nothing
end

"""

    csheuristics(wctl, udsd, writeeverynth; maxcs=100)

Find a reasonable value for the chunk size `cs` based on wallclocktime limit (DateTime),
ud-sweep duration `udsd` (in seconds) and `writeeverynth` to avoid

* non-progressing loop of resubmitting the simulation because we never make it to a dump
* memory overflow (we use a maximum chunk size `maxcs`)
* slowdown of simulation by IO (for too small chunk size)

A value of 0 for `wctl` is interpreted as no wallclocktime limit and chunk size `maxcs` will be returned.
"""
function csheuristics(wctl::DateTime, udsd::Real, writeeverynth::Int; maxcs::Int=100)
    wcts = (wctl - now()).value / 1000. # seconds till we hit the wctl
    # wcts == 0.0 && (return maxcs) # no wallclocktime limit

    num_uds = wcts/(udsd * writeeverynth) # Float64: number of push!s we will (in theory) perform in the given time
    num_uds *= 0.80 # 20% buffer, i.e. things might take longer than expected.
    cs = max(floor(Int, min(num_uds, 100)), 1)
end

"""
Choose a apropriate chunk size for the given Monte Carlo simulation.
"""
function choose_chunk_size(mc::AbstractDQMC)
    p = mc.p
    if !("udsweep" in keys(mc.a.to.inner_timers))
      @warn "No \"udsweep\" timer found. Choosing default `cs=100`!! Job might never dump to disk!"
      return 100
    end
    to_udsweep = mc.a.to["udsweep"]

    udsd = TimerOutputs.time(to_udsweep) *10.0^(-9)/TimerOutputs.ncalls(to_udsweep)
    cs = csheuristics(p.walltimelimit, udsd, p.write_every_nth)
    cs = min(cs, floor(Int, p.measurements/p.write_every_nth)) # cs musn't be larger than # measurements
    p.edrun && (cs = 1000) # this should probably better set maxcs in csheuristics call

    secs_to_dump = cs * p.write_every_nth * udsd
    dump_date = now() + Millisecond(ceil(Int, secs_to_dump*1000))
    println("Chose a chunk size of $cs. Should dump to file around $(formatdate(dump_date)). Walltime limit is $(formatdate(p.walltimelimit)).")

    return cs
end

formatdate(d) = Dates.format(d, "d.u yyyy HH:MM")

"""
Draw a random number from a uniform distribution over an interval [-b, b].
"""
@inline randuniform(b::Float64) = -b + 2 * b * rand() # taken from Distributions.jl: https://tinyurl.com/ycr9jnt4
randuniform(b::Float64, d::Int) = begin
    x = Vector{Float64}(undef, d)
    @inbounds for k in Base.OneTo(d)
      x[k] = randuniform(b)
    end
    x
end

"""
Manual and incomplete(!) serializer for JLD.

It only works for a single most inner timer for which `isempty(to.inner_timers) == true`.
"""
struct MyTimerOutputSerializer
    td::TimerOutputs.TimeData
    name::String
end

JLD.writeas(to::TimerOutput) = MyTimerOutputSerializer(to.accumulated_data, to.name)

function JLD.readas(tos::MyTimerOutputSerializer)
    to = TimerOutput()
    @timeit to tos.name 3+3 # something, will overwrite it
    to[tos.name].accumulated_data = tos.td
    to
end

"""
Allows one to use keyword syntax to add new entries to a `NamedTuple` `obs`.
"""
add(obs::NamedTuple; kw...) = combine(obs, values(kw))
"""
Combine two `NamedTuple`s into a new one.
"""
combine(x::NamedTuple, y::NamedTuple) = NamedTuple{(keys(x)..., keys(y)...)}((x...,y...))

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
#!/usr/bin/env julia

cmd = `th train.lua -input_h5 test.h5 -input_json test.json`
println(cmd)
run(cmd)


text:
println(PROGRAM_FILE);
for x in ARGS;
    println(x);
end


text:
mutable struct Eggbox_Ensemble <: GMC_NS_Ensemble
    path::String

    model_initλ::Function
    models::Vector{Eggbox_Record}

    contour::Float64
    log_Li::Vector{Float64}
    log_Xi::Vector{Float64}
    log_wi::Vector{Float64}
    log_Liwi::Vector{Float64}
    log_Zi::Vector{Float64}
    Hi::Vector{Float64}

    obs::Nothing
    priors::Vector{<:Distribution}
    constants::Vector{<:Real}
    box::Matrix{Float64}

    sample_posterior::Bool
    posterior_samples::Vector{Eggbox_Record}

    GMC_Nmin::Int64

    GMC_τ_death::Float64
    GMC_init_τ::Float64
    GMC_tune_μ::Int64
    GMC_tune_α::Float64
    GMC_tune_PID::NTuple{3,Float64}

    GMC_timestep_η::Float64
    GMC_reflect_η::Float64
    GMC_exhaust_σ::Float64
    GMC_chain_κ::Int64

    t_counter::Int64
end

Eggbox_Ensemble(path::String, no_models::Integer, prior, box, GS_settings...; sample_posterior::Bool=true) =
Eggbox_Ensemble(
    path,
    construct_eggbox_model,
    assemble_EMs(path, no_models, prior, box, Vector{Real}())...,
    [-Inf], #L0 = 0
	[0.], #X0 = 1
	[-Inf], #w0 = 0
	[-Inf], #Liwi0 = 0
	[-Inf], #Z0 = 0
	[0.], #H0 = 0,
    nothing,
    length(prior)==1 ? ([GMC_NS.marginals(prior)...]) : (prior),
    Vector{Real}(),
    length(prior)==1 ? (to_unit_ball.(box,[GMC_NS.marginals(prior)...])) : (to_unit_ball.(box,prior)),
    sample_posterior,
    Vector{Eggbox_Record}(),
    GS_settings...,
	no_models+1)

function Base.show(io::IO, e::Eggbox_Ensemble; progress=true)
    x1o=[0.0, 0.4, 0.8, 0.2, 0.6, 1.0, 0.0, 0.4, 0.8, 0.2, 0.6, 1.0, 0.0,0.4,0.8,0.2,0.6,1.0]
    x2o=[0.0, 0.0, 0.0, 0.2, 0.2, 0.2, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.8, 0.8, 0.8,1.0, 1.0, 1.0]

    x1m=[m.pos[1] for m in e.models]
    x2m=[m.pos[2] for m in e.models]
    x1m.+=1.; x2m.+=1.
    x1m./=2.; x2m./=2

    plt=scatterplot(x1o,x2o,title="Eggbox Ensemble: Contour $(round(e.contour,digits=2))",color=:green,name="Optima")
    scatterplot!(plt, x1m, x2m, color=:magenta, name="Model Particles")
    show(io, plt)
    println()

    (progress && return 19)
end

function assemble_EMs(path::String, no_trajectories::Integer, prior, box, constants)
	ensemble_records = Vector{Eggbox_Record}()
    !isdir(path) && mkpath(path)
    length(prior)==1 ? (marginals=[GMC_NS.marginals(prior)...]) : (marginals=prior)
    box=to_unit_ball.(box,marginals)

    @showprogress 1 "Assembling Normal Model ensemble..." for trajectory_no in 1:no_trajectories
		model_path = string(path,'/',trajectory_no,'.',1)
        if !isfile(model_path)
            proposal=[rand.(prior)...]
            pos=to_unit_ball.(proposal,marginals)
            box_bound!(pos,box)
            θvec=to_prior.(pos,marginals)

            model = construct_eggbox_model(trajectory_no, 1, θvec, pos, [0.], nothing, constants...; v_init=true)

			serialize(model_path, model) #save the model to the ensemble directory
			push!(ensemble_records, Eggbox_Record(trajectory_no,1,pos,model_path,model.log_Li))
		else #interrupted assembly pick up from where we left off
			model = deserialize(model_path)
			push!(ensemble_records, Eggbox_Record(trajectory_no,1,model.pos,model_path,model.log_Li))
		end
	end

	return ensemble_records, minimum([record.log_Li for record in ensemble_records])
end


text:
# Check for memory leaks
using DReal
reset_ctx!()

ctx = Context(qf_nra)
X = Var(ctx, Float64,"Xleaks4", 0.0,1.0)
for i = 1:1000
  push_ctx!(ctx)
  a = rand()
  b = rand()
  add!(ctx,(<=)(ctx,X,b))
  add!(ctx,(>=)(ctx,X,a))
  is_satisfiable(ctx)
  pop_ctx!(ctx)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiStorageV1alpha1VolumeError <: SwaggerModel
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    time::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: time

    function IoK8sApiStorageV1alpha1VolumeError(;message=nothing, time=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1alpha1VolumeError, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiStorageV1alpha1VolumeError, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type IoK8sApiStorageV1alpha1VolumeError

const _property_map_IoK8sApiStorageV1alpha1VolumeError = Dict{Symbol,Symbol}(Symbol("message")=>Symbol("message"), Symbol("time")=>Symbol("time"))
const _property_types_IoK8sApiStorageV1alpha1VolumeError = Dict{Symbol,String}(Symbol("message")=>"String", Symbol("time")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiStorageV1alpha1VolumeError }) = collect(keys(_property_map_IoK8sApiStorageV1alpha1VolumeError))
Swagger.property_type(::Type{ IoK8sApiStorageV1alpha1VolumeError }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiStorageV1alpha1VolumeError[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1alpha1VolumeError }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1alpha1VolumeError[property_name]

function check_required(o::IoK8sApiStorageV1alpha1VolumeError)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1alpha1VolumeError }, name::Symbol, val)
end


text:
ValidMethods = {
  :random_search => BlackBoxOptim.random_search,
  :de_rand_1_bin => BlackBoxOptim.de_rand_1_bin,
  :de_rand_2_bin => BlackBoxOptim.de_rand_2_bin,
  :de_rand_1_bin_radiuslimited => BlackBoxOptim.de_rand_1_bin_radiuslimited,
  :de_rand_2_bin_radiuslimited => BlackBoxOptim.de_rand_2_bin_radiuslimited,
  :adaptive_de_rand_1_bin => BlackBoxOptim.adaptive_de_rand_1_bin,
  :adaptive_de_rand_1_bin_radiuslimited => BlackBoxOptim.adaptive_de_rand_1_bin_radiuslimited,
  :separable_nes => BlackBoxOptim.separable_nes,
  :xnes => BlackBoxOptim.xnes,
  :resampling_memetic_search => BlackBoxOptim.resampling_memetic_searcher,
  :resampling_inheritance_memetic_search => BlackBoxOptim.resampling_inheritance_memetic_searcher,
  :simultaneous_perturbation_stochastic_approximation => BlackBoxOptim.SimultaneousPerturbationSA2,
  :generating_set_search => BlackBoxOptim.GeneratingSetSearcher,
  :probabilistic_descent => BlackBoxOptim.direct_search_probabilistic_descent,
}

MethodNames = collect(keys(ValidMethods))

# Default parameters for all convenience methods that are exported to the end user.
DefaultParameters = {
  :NumDimensions  => :NotSpecified, # Dimension of problem to be optimized
  :SearchRange    => (-10.0, 10.0), # Default search range to use per dimension unless specified
  :SearchSpace    => false, # Search space can be directly specified and takes precedence over Dimension and SearchRange if specified.

  :MaxTime        => false,   # Max time in seconds (takes precedence over the other budget-related params if specified)
  :MaxFuncEvals   => false,   # Max func evals (takes precedence over max iterations, but not max time)
  :MaxSteps       => 10000,   # Max iterations gives the least control since different optimizers have different "size" of their "iterations"
  :MinDeltaFitnessTolerance => 1e-11, # Minimum delta fitness (difference between two consecutive best fitness improvements) we can accept before terminating
  :FitnessTolerance => 1e-8,  # Stop optimization when the best fitness found is within this distance of the actual optimum (if known)

  :NumRepetitions => 1,     # Number of repetitions to run for each optimizer for each problem

  :ShowTrace      => true,  # Print tracing information during the optimization
  :TraceInterval  => 0.50,  # Minimum number of seconds between consecutive trace messages printed to STDOUT
  :SaveTrace      => false,
  :SaveFitnessTraceToCsv => false, # Save a csv file with information about the major fitness improvement events (only the first event in each fitness magnitude class is saved)
  :SaveParameters => false, # Save parameters to a json file for later scrutiny

  :RandomizeRngSeed => true, # Randomize the RngSeed value before using any random numbers.
  :RngSeed        => 1234,   # The specific random seed to set before any random numbers are generated. The seed is randomly selected if RandomizeRngSeed is true, and this parameter is updated with its actual value.

  :PopulationSize => 50
}

# Create a problem given
#   a problem or
#   a function and a search range or
#   a function and a
# while possibly updating the params with the specific dimension and search
# space to be used.
function setup_problem(functionOrProblem; parameters = Dict())

  params = Parameters(parameters, DefaultParameters)

  # If an OptimizationProblem was given it takes precedence over the search range param setting.
  if issubtype(typeof(functionOrProblem), OptimizationProblem)

    # If a fixed dim problem was given it takes precedence over the dimension param setting.
    if is_fixed_dimensional(functionOrProblem)
      problem = functionOrProblem
    else
      # If an anydim problem was given the dimension param must have been specified.
      if params[:NumDimensions] == :NotSpecified
        throw(ArgumentError("You MUST specify the number of dimensions in a solution when an any-dimensional problem is given"))
      else
        problem = as_fixed_dim_problem(functionOrProblem, parameters[:NumDimensions])
      end
    end

  elseif typeof(functionOrProblem) == Function

    # Check that a valid search space has been stated and create the search_space
    # based on it, or bail out.
    if typeof(params[:SearchRange]) == typeof( (0.0, 1.0) )
      if params[:NumDimensions] == :NotSpecified
        throw(ArgumentError("You MUST specify the number of dimensions in a solution when giving a search range $(searchRange)"))
      end
      params[:SearchSpace] = symmetric_search_space(params[:NumDimensions], params[:SearchRange])
    elseif typeof(params[:SearchRange]) == typeof( [(0.0, 1.0)] )
      params[:SearchSpace] = RangePerDimSearchSpace(params[:SearchRange])
      params[:NumDimensions] = length(params[:SearchRange])
    else
      throw(ArgumentError("Invalid search range specification."))
    end

    # Now create an optimization problem with the given information. We currently reuse the type
    # from our pre-defined problems so some of the data for the constructor is dummy.
    problem = fixeddim_problem(functionOrProblem;
      search_space = params[:SearchSpace], range = params[:SearchRange],
      dims = params[:NumDimensions]
    )

  end

  params[:SearchSpace] = search_space(problem)

  return problem, params

end

function compare_optimizers(functionOrProblem::Union(Function, OptimizationProblem);
  max_time = false, search_space = false, search_range = (0.0, 1.0), dimensions = 2,
  methods = MethodNames, parameters = Dict())

  evaluator = 1.0

  results = Any[]
  for(m in methods)
    tic()
    best, fitness, reason, etime, parameters = bboptimize(functionOrProblem; method = m, parameters = parameters,
      max_time = max_time, search_space = search_space, dimensions = dimensions,
      search_range = search_range)
    push!( results,  (m, best, fitness, toq()) )
    evaluator = parameters[:Evaluator]
    delete!(parameters, :Evaluator)
  end

  sorted = sort( results, by = (t) -> t[3] )

  if parameters[:ShowTrace]
    println("\n********************************************************************************")
    println(describe(evaluator))
    for(i in 1:length(sorted))
      println("$(i). $(sorted[i][1]), fitness = $(sorted[i][3]), time = $(sorted[i][4])")
    end
    println("********************************************************************************\n")
  end

  return sorted

end

function compare_optimizers(problems::Dict{Any, FixedDimProblem}; max_time = false,
  methods = MethodNames, parameters = Dict())

  # Lets create an array where we will save how the methods ranks per problem.
  ranks = zeros(length(methods), length(problems))
  fitnesses = zeros(Float64, length(methods), length(problems))
  times = zeros(Float64, length(methods), length(problems))

  problems = collect(problems)

  for i in 1:length(problems)
    name, p = problems[i]
    res = compare_optimizers(p; max_time = max_time, methods = methods, parameters = parameters)
    for(j in 1:length(res))
      method, best, fitness, elapsedtime = res[j]
      index = findfirst(methods, method)
      ranks[index, i] = j
      fitnesses[index, i] = fitness
      times[index, i] = elapsedtime
    end
  end

  avg_ranks = round(mean(ranks, 2), 2)
  avg_fitness = round(mean(fitnesses, 2), 3)
  avg_times = round(mean(times, 2), 2)

  perm = sortperm(avg_ranks[:])
  println("\nBy avg rank:")
  for(i in 1:length(methods))
    j = perm[i]
    print("\n$(i). $(methods[j]), avg rank = $(avg_ranks[j]), avg fitness = $(avg_fitness[j]), avg time = $(avg_times[j]), ranks = ")
    showcompact(ranks[j,:][:])
  end

  perm = sortperm(avg_fitness[:])
  println("\n\nBy avg fitness:")
  for(i in 1:length(methods))
    j = perm[i]
    print("\n$(i). $(methods[j]), avg rank = $(avg_ranks[j]), avg fitness = $(avg_fitness[j]), avg time = $(avg_times[j]), ranks = ")
    showcompact(ranks[j,:][:])
  end

  return ranks, fitnesses
end

function bboptimize(functionOrProblem; max_time = false,
  search_space = false, search_range = (0.0, 1.0), dimensions = 2,
  method = :adaptive_de_rand_1_bin_radiuslimited,
  parameters = Dict())

  params = Parameters(parameters, DefaultParameters)
  params[:MaxTime] = max_time
  params[:SearchSpace] = search_space
  params[:SearchRange] = search_range
  params[:NumDimensions] = dimensions

  problem, params = setup_problem(functionOrProblem; parameters = params)

  # Create a random solution from the search space and ensure that the given function returns a Number.
  ind = rand_individual(params[:SearchSpace])
  res = eval1(ind, problem)
  if !(typeof(res) <: Number)
    throw(ArgumentError("The supplied function does NOT return a number when called with a potential solution (when called with $(ind) it returned $(res)) so we cannot optimize it!"))
  end

  # Check that max_time is larger than zero if it has been specified.
  if params[:MaxTime] != false
    if params[:MaxTime] <= 0.0
      throw(ArgumentError("The max_time must be a positive number"))
    else
      params[:MaxTime] = convert(Float64, params[:MaxTime])
    end
  end

  # Check that a valid number of iterations has been specified. Print warning if higher than 1e8.
  if params[:MaxFuncEvals] != false
    if params[:MaxFuncEvals] < 1
      throw(ArgumentError("The number of function evals MUST be a positive number"))
    elseif params[:MaxFuncEvals] >= 1e8
      println("Number of allowed function evals is $(params[:MaxFuncEvals]); this can take a LONG time")
    end
  end

  # Check that a valid number of iterations has been specified. Print warning if higher than 1e8.
  if params[:MaxSteps] < 1
    throw(ArgumentError("The number of iterations MUST be a positive number"))
  elseif params[:MaxSteps] >= 1e7
    println("Number of allowed iterations is $(params[:MaxSteps]); this can take a LONG time")
  end

  # Check that a valid population size has been given.
  if params[:PopulationSize] < 2
    throw(ArgumentError("The population size MUST be at least 2"))
  end

  # Check that a valid method has been specified and then set up the optimizer
  if (typeof(method) != Symbol) || !any([(method == vm) for vm in MethodNames])
    throw(ArgumentError("The method specified, $(method), is NOT among the valid methods: $(MethodNames)"))
  end
  pop = BlackBoxOptim.rand_individuals_lhs(params[:SearchSpace], params[:PopulationSize])

  params = Parameters(params, {
    :Evaluator    => ProblemEvaluator(problem),
    :Population   => pop,
    :SearchSpace  => search_space
  })
  optimizer_func = ValidMethods[method]
  optimizer = optimizer_func(params)

  run_optimizer_on_problem(optimizer, problem; parameters = params)
end

function tr(msg, parameters, obj = None)
  if parameters[:ShowTrace]
    print(msg)
    if obj != None
      showcompact(transform_horizontal(obj))
    end
  end
  if parameters[:SaveTrace]
    # No saving for now
  end
end

function find_best_individual(e::Evaluator, opt::PopulationOptimizer)
  (best_candidate(e.archive), 1, best_fitness(e.archive))
end

function find_best_individual(e::Evaluator, opt::Optimizer)
  (best_candidate(e.archive), 1, best_fitness(e.archive))
end

# Transform a vector to ensure it is largest in the horizontal direction.
function transform_horizontal(x)
  if typeof(x) <: Array
    nrows, ncols = size(x)
    if ncols < nrows
      x = x'
    end
  end
  return x
end

function run_optimizer_on_problem(opt::Optimizer, problem::OptimizationProblem;
  parameters = Dict())

  if parameters[:RandomizeRngSeed]
    parameters[:RngSeed] = rand(1:int(1e6))
    srand(parameters[:RngSeed])
  end

  # No max time if unspecified. If max time specified it takes precedence over
  # max_steps and MaxFuncEvals. If no max time MaxFuncEvals takes precedence over
  # MaxSteps.
  if parameters[:MaxTime] == false
    max_time = Inf
    if parameters[:MaxFuncEvals] != false
      max_fevals = parameters[:MaxFuncEvals]
      max_steps = Inf
    else
      max_steps = parameters[:MaxSteps]
      max_fevals = Inf
    end
  else
    max_steps = Inf
    max_fevals = Inf
    max_time = parameters[:MaxTime]
  end

  # Now set up an evaluator for this problem. This will handle fitness
  # comparisons, keep track of the number of function evals as well as
  # keep an archive and top list.
  evaluator = get(parameters, :Evaluator, ProblemEvaluator(problem))

  num_better = 0
  num_better_since_last = 0
  tr("Starting optimization with optimizer $(name(opt))\n", parameters)

  step = 1
  t = last_report_time = start_time = time()
  elapsed_time = 0.0

  termination_reason = "" # Will be set in loop below...

  while( true )

    if elapsed_time > max_time
      termination_reason = "Max time reached"
      break
    end

    if num_evals(evaluator) > max_fevals
      termination_reason = "Max number of function evaluations reached"
      break
    end

    if step > max_steps
      termination_reason = "Max number of steps reached"
      break
    end

    if delta_fitness(evaluator.archive) < parameters[:MinDeltaFitnessTolerance]
      termination_reason = "Delta fitness below tolerance"
      break
    end

    if fitness_is_within_ftol(evaluator, parameters[:FitnessTolerance])
      termination_reason = "Within fitness tolerance of optimum"
      break
    end

    # Report on progress every now and then...
    if (t - last_report_time) > parameters[:TraceInterval]
      last_report_time = t
      num_better += num_better_since_last

      # Always print step number, num fevals and elapsed time
      tr(@sprintf("%.2f secs, %d evals , %d steps",
        elapsed_time, num_evals(evaluator), step), parameters)

      # Only print if this optimizer reports on number of better. They return 0
      # if they do not.
      if num_better_since_last > 0
        tr(@sprintf(", improv/step: %.3f (last = %.4f)",
          num_better/step, num_better_since_last/step), parameters)
        num_better_since_last = 0
      end

      # Always print fitness if num_evals > 0
      if num_evals(evaluator) > 0
        tr(@sprintf(", %.9f", best_fitness(evaluator.archive)), parameters)
      end

      tr("\n", parameters)
    end

    if has_ask_tell_interface(opt)

      # They ask and tell interface is more general since you can mix and max
      # elements from several optimizers using it. However, in this top-level
      # execution function we do not make use of this flexibility...
      candidates = ask(opt)
      ranked_candidates = BlackBoxOptim.rank_by_fitness(evaluator, candidates)
      num_better_since_last += tell!(opt, ranked_candidates)

    else

      BlackBoxOptim.step(opt)
      num_better_since_last = 0

    end

    step += 1
    t = time()
    elapsed_time = t - start_time
  end

  step -= 1 # Since it is one too high after while loop above

  tr("\nOptimization stopped after $(step) steps and $(elapsed_time) seconds\n", parameters)
  tr("Termination reason: $(termination_reason)\n", parameters)
  tr("Steps per second = $(step/elapsed_time)\n", parameters)
  tr("Function evals per second = $(num_evals(evaluator)/elapsed_time)\n", parameters)
  tr("Improvements/step = $((num_better+num_better_since_last)/max_steps)\n", parameters)
  tr("Total function evaluations = $(num_evals(evaluator))\n", parameters)

  if typeof(opt) <: PopulationOptimizer
    tr("\nMean value (in population) per position:", parameters, mean(population(opt),1))
    tr("\n\nStd dev (in population) per position:", parameters, std(population(opt),1))
  end

  best, index, fitness = find_best_individual(evaluator, opt)
  tr("\n\nBest candidate found: ", parameters, best)
  tr("\n\nFitness: ", parameters, fitness)
  tr("\n\n", parameters)

  if parameters[:SaveFitnessTraceToCsv]
    timestamp = strftime("%y%m%d_%H%M%S", ifloor(start_time))
    filename = "$(timestamp)_$(problem_summary(evaluator))_$(name(opt)).csv"
    filename = replace(replace(filename, r"\s+", "_"), r"/", "_")
    header_prefix = "Problem,Dimension,Optimizer"
    line_prefix = "$(name(problem)),$(numdims(problem)),$(name(opt))"
    save_fitness_history_to_csv_file(evaluator.archive, filename;
      header_prefix = header_prefix, line_prefix = line_prefix,
      bestfitness = fmin(problem))
  end

  return transform_horizontal(best), fitness, termination_reason, elapsed_time, parameters, num_evals(evaluator)
end

# Summarize a vector of float values by stating its mean, std dev and median.
function report_on_values(desc, v, lpad = "", rpad = "", digits = 3)
  println("$(lpad)$(desc): $(signif(mean(v), digits)) (std. dev = $(signif(std(v), digits)), median = $(signif(median(v), digits)))")
end

# Report on the number of times each key in a count dict was encountered.
# Returns a percentage dict calculated while iterating over the counted items.
function count_dict_report(dict, desc, lpad = "", rpad = "")
  println(desc, ":")
  total = sum(collect(values(dict)))
  pdict = Dict()
  for (r, c) in dict
    pdict[r] = round(100.0*c/total, 2)
    println(lpad, r, ": ", c, " (", pdict[r], "%)", rpad)
  end
  pdict
end

# Print a report based on a result dict from one set of repeated runs of
# an optimization method. Returns the success rate, i.e. number of times the
# termination reason was "Within fitness tolerance...".
function report_from_result_dict(statsdict)
  println("Method: $(statsdict[:method])")
  pdict = count_dict_report(statsdict[:reasoncounts], "  Termination reasons", "    ")
  report_on_values("Fitness", statsdict[:fitnesses], "  ")
  report_on_values("Time", statsdict[:times], "  ")
  report_on_values("Num function evals", statsdict[:numevals], "  ")
  println("  Success rate: ", round(pdict["Within fitness tolerance of optimum"], 3), "%\n")
  pdict["Within fitness tolerance of optimum"]
end

function rank_result_dicts_by(result_dicts, byfunc, desc; rev = false,
  descsummary = "mean", digits = 3, rpad = "")

  ranked = BlackBoxOptim.Utils.assign_ranks_within_tolerance(result_dicts; by = byfunc, tolerance = 1e-3, rev = rev)
  println("Ranked by $(descsummary) $(desc):")
  for (rank, rd, value) in ranked
    println("  $(rank). $(rd[:method]), $(signif(value, digits))$(rpad)")
  end
  println("")

end

function report_on_methods_results_on_one_problem(problem, result_dicts, numrepeats, max_time, ftol)

  println("********************************************************************************\n")

  println("Problem: $(name(problem)), num dims = $(numdims(problem))")
  println("  Num repeats per method = ", numrepeats)
  println("  Fitness tolerance = ", ftol, " (a run is a success if it reaches to within this value of true optimum)")
  println("  Max time budget per run = ", max_time, " secs\n")

  rank_result_dicts_by(result_dicts, (d) -> d[:success_rate], "success rate (to reach within $(ftol) of optimum)";
    descsummary = "median", rev = true, rpad = "%")
  rank_result_dicts_by(result_dicts, (d) -> median(d[:fitnesses]), "fitness"; descsummary = "median")
  rank_result_dicts_by(result_dicts, (d) -> median(d[:times]), "time (in seconds)";
    descsummary = "median", rpad = " secs")
  rank_result_dicts_by(result_dicts, (d) -> int(median(d[:numevals])), "num function evals"; descsummary = "median")

  for rd in result_dicts
    report_from_result_dict(rd)
  end

end

function repeated_bboptimize(numrepeats, problem, dim, methods, max_time, ftol = 1e-5, parameters = Dict{Any, Any}())

  fp = BlackBoxOptim.as_fixed_dim_problem(problem, dim)
  result_dicts = Any[]

  # Just so they are declared
  ps = best_so_far = nothing

  params = Parameters(parameters, {:FitnessTolerance => ftol})

  for m in methods

    ts, fs, nes = zeros(numrepeats), zeros(numrepeats), zeros(Int64, numrepeats)
    rcounts = {"Within fitness tolerance of optimum" => 0}

    for i in 1:numrepeats
      p = fp # BlackBoxOptim.ShiftedAndBiasedProblem(fp)
      best, fs[i], reason, ts[i], ps, nes[i] = bboptimize(p; max_time = max_time,
        method = m, parameters = params)
      rcounts[reason] = 1 + get(rcounts, reason, 0)
    end

    if best_so_far == nothing
      best_so_far = worst_fitness(ps[:Evaluator])
    end

    best_so_far =

    rdict = {:method => m, :fitnesses => fs, :times => ts, :numevals => nes, :reasoncounts => rcounts}
    rdict[:success_rate] = report_from_result_dict(rdict)
    push!(result_dicts, rdict)

  end

  report_on_methods_results_on_one_problem(fp, result_dicts, numrepeats, max_time, ftol)

end


text:
function solve()
    a, b, c, d = [parse(Int, x) for x in split(readline())]

    takahashi = (c+b-1) ÷ b
    aoki = (a+d-1) ÷ d
    takahashi <= aoki ? "Yes" : "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
@testset "subpixel nonmaxima supression" begin
    @info "Test: SubpixelNonmaximaSuppression"

    @testset "Numerical" begin
        img = Gray{N0f8}.(load("algorithms/References/circle.png"))
        nms_ref = Gray{N0f8}.(load("algorithms/References/circle_nms.png"))

        f = SubpixelNonmaximaSuppression()

        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        mag = hypot.(g₁, g₂)

        nms = zeros(eltype(mag), axes(mag))
        offsets = zeros(SVector{2,Float64}, axes(mag))

        f(nms, offsets, mag, g₁, g₂)

        for i in CartesianIndices(nms_ref)
            nms_ref[i] == 0.0 ? (@test nms[i] == 0) : (@test nms[i] == mag[i])
            nms[i] == 0.0 ? (@test offsets[i] == SVector(0.0, 0.0)) : (@test offsets[i] != SVector(0.0, 0.0))
        end
    end
end


text:
module PmapProgressMeter

using ProgressMeter
using Compat

# a global to hold progress meter references
globalProgressMeters = Dict()
globalProgressValues = Dict()
globalPrintLock = Dict()

"Wraps pmap with a progress meter. If the keyword argument passcallback is true, then passes an additional first argument to the function which is a callback taking one argument, n, to add n to the progress meter. If it is false, then updates the progress meter by 1 every time a parallel call
finishes."
function Base.pmap(f::Function, p::Progress, values...; kwargs...)
    global globalProgressMeters
    global globalProgressValues
    global globalPrintLock

    id = randstring(50)
    globalProgressMeters[id] = p
    globalProgressValues[id] = p.counter
    globalPrintLock[id] = ReentrantLock()

    passcallback = false
    kwa = Dict(kwargs)
    if haskey(kwa,:passcallback)
      passcallback = true
      delete!(kwa,:passcallback)
    end

    out = pmap(values...; kwa...) do x...
        if passcallback
          v = f(n -> remotecall(updateProgressMeter, 1, id, n), x...)
        else
          v = f(x...)
          wait(remotecall(updateProgressMeter, 1, id, 1))
        end
        v
    end

    delete!(globalProgressMeters, id)
    out
end

"This is remote-called by all the workers to update the progress."
function updateProgressMeter(id,n)
    global globalProgressMeters
    global globalProgressValues
    global globalPrintLock

    lock(globalPrintLock[id])
    globalProgressValues[id] += n
    update!(globalProgressMeters[id] , globalProgressValues[id])
    unlock(globalPrintLock[id])
end

end # module

max_stars_repo_name: slundberg/PmapProgressMeter

max_issues_repo_name: slundberg/PmapProgressMeter

max_forks_repo_name: slundberg/PmapProgressMeter


text:
Base.call{T <: SymbolicObject}(ex::T) = ex
function Base.call{T <: SymbolicObject}(ex::T, args...)
    xs = free_symbols(ex)
    if length(xs) >= 1
        subs(ex, collect(zip(xs, args))...)
    else
        if classname(ex) == "Symbol"
            ex
        else
            pycall(PyObject(ex), PyAny, args...)
        end
    end
end
Base.call(ex::SymbolicObject, x::Dict) = subs(ex, x)
Base.call(ex::SymbolicObject, x::Pair...) = subs(ex, x...)

## for symbolic functinos (dsolve)
Base.call(u::SymFunction, x::Base.Dict) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
Base.call(u::SymFunction, x::Base.Pair) = throw(ArgumentError("IVPsolutions can only be called with symbolic objects"))
function Base.call(u::SymFunction, x::Sym)
    if u.n == 0
        u.x(PyObject(x))
    else
        __x = Sym("__x")
        diff(u.x(PyObject(__x)), __x, u.n)(__x => x)
    end
end

function Base.call(u::SymFunction, x::Sym, y...)
    if u.n== 0
        u.x(vcat(x, y...)...)
    else
        error("Need to implement derivatives of symbolic functions of two or more variables")
    end
end


text:
using LinearAlgebra
using SparseArrays

function expmv(t::StepRangeLen, A, b::AbstractVecOrMat;
                M = nothing, precision = "double", shift = false)

    t0 = Float64(t.ref)
    tmax = Float64(t.ref + (t.len - 1.) * t.step)
    q = t.len - 1
    n = size(A, 1)

    if shift == true && !hasmethod(tr, typeof(A))
        shift = false
        @warn "Shift set to false as $(typeof(A)) doesn't support tr"
    end

    force_estm = !hasmethod(opnorm, Tuple{typeof(A), typeof(1)})
    temp = (tmax - t0) * normAm(A, 1);

    if (precision == "single" || precision == "half" && temp > 85.496) || ( precision == "double" && temp > 63.152)
       force_estm = true;
    end

    if M == nothing
        (M, alpha, unA) = select_taylor_degree(A, b, precision=precision, shift=shift, force_estm=force_estm)
    end

    tol =
      if precision == "double"
          2.0^(-53)
      elseif precision == "single"
          2.0^(-24)
      elseif precision == "half"
          2.0^(-10)
      end

    X = zeros(eltype(A), n, q+1);
    (m_max, p) = size(M);
    U = diagm(0 => 1:m_max);

    temp, s = degree_selector(tmax - t0, M, U, p)
    h = (tmax - t0)/q;

    X[:,1] = expmv(t0,A,b,M=M,precision=precision,shift=shift);

    mu = 0.
    if shift
        mu = tr(A)/n
        A = A - mu * I
    end

    d = max(1, Integer(floor(q/s)))
    j = Integer(floor(q/d))
    r = q - d * j
    z = X[:,1]
    m_opt, = degree_selector(d, M, U, p)
    dr = d
    K = zeros(eltype(A), n, m_opt+1)

    temp = similar(b)
    f = similar(b)

    for i = 1:j+1
        if i > j
            dr = r
        end
        K[:,1] = z
        m = 0
        for k = 1:dr
            f = copy(z)
            c1 = norm(z, Inf)

            p = 1
            for outer p = 1:m_opt
                if p > m
                    @views K[:,p+1] .= (h/p).*(A*K[:,p])
                end

                @views temp .= (Float64(k)^p).*K[:,p+1]
                f .+= temp
                c2 = norm(temp, Inf)
                if c1 + c2 <= tol * norm(f, Inf)
                    break
                end
                c1 = c2
            end

            m = max(m,p)
            @views X[:, k + (i - 1) * d + 1] .= exp(k * h * mu) * f
        end

        if i <= j
            z = X[:,i*d+1]
        end
    end
    return X
end