text:
#                       The display module for DomTerm.jl
#
# lend idea and code from the `IJulia.jl` package
#
#       J.Wang, 2021

import Base: show, display, redisplay

using Base64: Base64EncodePipe

struct DomTermInline <: AbstractDisplay end

# supported MIME types for DomTerm inline display, in descending order
# of preference (descending "richness")
const DT_mime = [
    #"application/vnd.dataresource+json",
    #["application/vnd.vegalite.v$n+json" for n in 4:-1:2]...,
    #["application/vnd.vega.v$n+json" for n in 5:-1:3]...,
    #"application/vnd.plotly.v1+json",
    "text/html",
    "text/mathml+xml",
    "text/latex",
    "image/svg+xml",
    "image/png",
    "image/jpeg",
    #"application/javascript",
    #"text/markdown",
    "text/plain"
]

DTInlineIOContext(io, KVs::Pair...) = IOContext(
    io,
    :limit=>true, :color=>true, :domterm=>true,
    KVs...
)

# need special handling for showing a string as a textmime
# type, since in that case the string is assumed to be
# raw data unless it is text/plain (JW: handled by the show() function)
# JW: For any type signature that you want to implement in show()
# let the `israwtext()` function with the same signature return false
israwtext(::MIME"text/plain", x::AbstractString) = false
israwtext(::MIME"text/html", x::AbstractString) = false
israwtext(::MIME"text/mathml+xml", x::AbstractString) = false
israwtext(::MIME, x::AbstractString) = true
israwtext(::MIME, x) = false

# convert x to a string of type mime, making sure to use an
# IOContext that tells the underlying show function to limit output
function limitstringmime(mime::MIME, x)
    buf = IOBuffer()
    if istextmime(mime)
        if israwtext(mime, x)
            return String(x)
        else
            # If the raw data `x` is really large and `x` is shown as
            # "image/svg+xml", The `show()` function defined in `ImageShow()`
            # emits tons of SVG commands, sometimes hanges `DomTerm`. We
            # temporarily disable showing large inline SVG images. Use an
            # external viewer instead or show it as PNG.
            if applicable(length, x) && length(x) > 10^5
                print("data too large to show inline. If this is an image, try show it as PNG, or use an external viewer.")
                return ""
            else
                show(DTInlineIOContext(buf), mime, x)
            end
        end
    else
        b64 = Base64EncodePipe(buf)
        if isa(x, Vector{UInt8})
            write(b64, x) # x assumed to be raw binary data
        else
            # if the input array `x` is really large (length > 4M, which means
            # even down-sampled array is still more than 1M long), The `show()`
            # function defined in `ImageShow` emits tons of base64 bytes and
            # hanges `DomTerm`. We temporarily disable showing large inline
            # images. I'm not sure the following 1M bar is proper or not...
            if applicable(length, x) && length(x) > 1*10^6
                print("data too large to show inline. Please use an external viewer.")
                close(b64)
                return ""
            else
                show(DTInlineIOContext(b64), mime, x)
            end
        end
        close(b64)
    end
    return String(take!(buf))
end

# deal with annoying application/x-latex == text/latex synonyms
display(d::DomTermInline, m::MIME"application/x-latex", x) =
    display(d, MIME("text/latex"), limitstringmime(m, x))

# deal with annoying text/javascript == application/javascript synonyms
display(d::DomTermInline, m::MIME"text/javascript", x) =
    display(d, MIME("application/javascript"), limitstringmime(m, x))

function display(d::DomTermInline, M::MIME, x)
    sx = limitstringmime(M, x)
    d = Dict(string(M) => sx)
    if istextmime(M)
        d["text/plain"] = sx # directly show text data, e.g. text/csv
    end

    # JW: I should print it more properly...
    print(sx)
end
displayable(d::DomTermInline, M::MIME) = istextmime(M)

# define `display` function for all supported MIMEs
for mime in DT_mime
    @eval begin
        function display(d::DomTermInline, ::MIME{Symbol($mime)}, x)
            print("\e]72;") # wire byte code for inserting a HTML fragment
            publish_mime_string(MIME($mime), limitstringmime(MIME($mime), x))
            print("\a") # HTML fragment end
        end
        displayable(d::DomTermInline, ::MIME{Symbol($mime)}) = true
    end
end
# register single-argument function `display(x)` for html-type objects
display(d::DomTermInline, x::HTML) = display(d, MIME("text/html"), x)

function publish_mime_string(::MIME"image/svg+xml", x::String)
    # Remove svg file headers. Package `Colors` generates these headers in
    # running `display(MIME("image/svg+xml"), RGB{N0f8}(...))`. DomTerm
    # can not understand them.
    x = replace(x, r"<\?xml.*?>"s=>"")
    x = replace(x, r"<!DOCTYPE.*?>"s=>"")
    x = replace(x, r"\n"=>"") # strip extra spaces off the output
    print(x)
end

function publish_mime_string(::MIME"image/png", x::String)
    if !isempty(x)
        # if the input data used too much memory, x could be empty
        print("<img src='data:image/png;base64,")
        print(x)
        print("' />")
    end
end

publish_mime_string(::MIME, s::String) =  print(s)


text:
@enum InputCategory begin
    BRANCH
    BUS
    DC_BRANCH
    GENERATOR
    LOAD
    RESERVE
    SIMULATION_OBJECTS
    STORAGE
end

ENUM_MAPPINGS = Dict()
for enum in (AngleUnits.AngleUnit, BusTypes.BusType, InputCategory)
    ENUM_MAPPINGS[enum] = Dict()
    for value in instances(enum)
        ENUM_MAPPINGS[enum][lowercase(string(value))] = value
    end
end

"""Get the enum value for the string. Case insensitive."""
function get_enum_value(enum, value::String)
    if !haskey(ENUM_MAPPINGS, enum)
        throw(ArgumentError("enum=$enum is not valid"))
    end

    val = lowercase(value)
    if !haskey(ENUM_MAPPINGS[enum], val)
        throw(ArgumentError("enum=$enum does not have value=$val"))
    end

    return ENUM_MAPPINGS[enum][val]
end

ext: jl


text:
"""
Read Stan samples from a CSV file. Columns that belong to the same variable are grouped into
arrays.

The exported functions are [`read_samples`](@ref) and [`read_sample_matrix`](@ref).
"""
module StanSamples

export read_samples, read_sample_matrix

using ArgCheck: @argcheck
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using ElasticArrays: ElasticArray

####
#### utilities
####

"""
$(SIGNATURES)

Test if the argument is a comment line in Stan sample output.
"""
iscommentline(s::String) = occursin(r"^ *#", s)

"""
$(SIGNATURES)

Return the fields of a line in Stan sample output. The format is CSV,
but never quoted or escaped, so splitting on `,` is sufficient.
"""
fields(line::String) = split(chomp(line), ',')

"""
$(SIGNATURES)

Specification of a variable in the column of the posterior sample.

# Fields

$(FIELDS)
"""
struct ColVar{N}
    "variable name"
    name::Symbol
    "index (may be empty)"
    index::CartesianIndex{N}
end

ColVar(name::Symbol, index::Int...) = ColVar(name, CartesianIndex(tuple(index...)))

"""
$(SIGNATURES)

Test if two `ColVar`s can be merged (same `name` and number of indices).
"""
≅(::ColVar, ::ColVar) = false
≅(v1::ColVar{N}, v2::ColVar{N}) where {N} = v1.name == v2.name

"""
$(SIGNATURES)

Parse a string as a column variable.
"""
function ColVar(s::AbstractString)
    s = split(s, ".")
    name = Symbol(s[1])
    indexes = parse.(Int, s[2:end])
    @argcheck all(indexes .≥ 1) "Non-positive index in $(s)."
    ColVar(name, indexes...)
end

"""
    $(SIGNATURES)

For a vector of indexes, calculate the size (the largest one) and
check that they are contiguous and column-major. Return a tuple of
`Int`s (empty for scalars.)
"""
function combined_size(indexes)
    siz = reduce(max, indexes)
    ran = CartesianIndices(siz)
    # FIXME inelegant collect below
    @argcheck collect(indexes) == vec(collect(ran)) "Non-contiguous indexes."
    siz.I
end

"""
A variable denoting a Stan value, combined from adjacent columns of
with the same variable name. Always has a `name::Symbol`
field. Determines the type of the resulting values.
"""
abstract type StanVar end

"""
$(TYPEDEF)

A scalar (always Float64).
"""
struct StanScalar <: StanVar
    name::Symbol
end

"""
$(TYPEDEF)

An array (always of Float64 elements).
"""
struct StanArray{N} <: StanVar
    name::Symbol
    size::NTuple{N, Int}
end

# this is a shorthand, mainly useful for unit tests
StanArray(name::Symbol, size::Int...) = StanArray(name, size)

"""
$(SIGNATURES)

Empty values for a `var`, can be appended to using `append!`.
"""
empty_values(var::StanScalar) = Vector{Float64}()
empty_values(var::StanArray) = ElasticArray{Float64}(undef, var.size..., 0)

"""
$(SIGNATURES)

Number of columns that correspond to a [`StanVar`](@ref).
"""
ncols(::StanScalar) = 1
ncols(sa::StanArray) = prod(sa.size)

"""
$(SIGNATURES)

Combine column variables into a Stan variable.
"""
function _combine_colvars(colvars)
    var = first(colvars)
    len = findfirst(v -> !(v ≅ var), colvars)
    len = len ≡ nothing ? length(colvars) : len - 1
    siz = combined_size(v.index for v in colvars[1:len])
    if isempty(siz)
        StanScalar(var.name)
    else
        StanArray(var.name, siz)
    end
end

"""
$(SIGNATURES)

Combine column variables, returning a `Tuple` of `StanVar`s.
"""
function combine_colvars(colvars)
    header = StanVar[]
    position = 1
    while position ≤ length(colvars)
        v = _combine_colvars(@view colvars[position:end])
        @argcheck v.name ∉ (h.name for h in header) "Duplicate variable $(v.name)."
        position += ncols(v)
        push!(header, v)
    end
    tuple(header...)
end

"""
$(SIGNATURES)

Read values for `var` from `buffer`, starting at index `1`.
"""
_read_values(var::StanScalar, buffer) = buffer[1]

function _read_values(var::StanArray, buffer)
    a = Array{Float64}(undef, var.size...)
    a[:] .= buffer[1:length(a)]
    a
end

"""
$(SIGNATURES)

Create an empty container for variable values, accessed by variable names.
"""
function empty_vars_values(vars::Tuple)
    NamedTuple{map(var -> var.name, vars)}(map(empty_values, vars))
end

"""
$(SIGNATURES)

Read values from a single line of `io` using the variable
specification `vars`.

The fields are combined into variables and appended into the
corresponding vectors in `vars_values`.

Return `false` for comment lines, `true` lines with data. All other
cases (ie incomplete lines) throw an error. Note that in this case the
vectors in `vars_values` may have an inconsistent length.

`buffer` is pre-allocated for reading in a line at a time.
"""
function read_values(io::IO, vars::Tuple, vars_values::NamedTuple,
                     buffer::Vector{Float64} = Vector{Float64}(undef, length(vars)))
    line = readline(io)
    iscommentline(line) && return false
    buffer .= parse.(Float64, fields(line))
    position = 1
    for var in vars
        a = _read_values(var, @view buffer[position:end])
        append!(vars_values[var.name], a)
        position += ncols(var)
    end
    @assert position == length(buffer) + 1 "Fields remaining after parsing."
    true
end

"""
$(SIGNATURES)

Helper function to read data from a Stan samples CSV file.
"""
function _read_samples(io, vars, vars_values, buffer)
    while !eof(io)
        read_values(io, vars, vars_values, buffer)
    end
    vars_values
end

function read_samples(io::IO)
    while !eof(io)
        line = readline(io)
        if !iscommentline(line)
            colvars = ColVar.(fields(line))
            vars = combine_colvars(colvars)
            vars_values = empty_vars_values(vars)
            buffer = Vector{Float64}(undef, sum(ncols, vars))
            return _read_samples(io, vars, vars_values, buffer)
        end
    end
    error("Could not find non-empty lines.")
end

"""
$(SIGNATURES)

Read Stan samples from a CSV file or a `IO` stream.

Return a container of arrays, accessed by variables names as properties (eg like a
`NamedTuple`, which is in fact the current implementation, but can possibly change). Each
array has samples for a variable, with the last index varying for each draw.

```jldoctest
julia> io = IOBuffer("a,b.1,b.2,c.1.1,c.2.1,c.1.2,c.2.2\n" *
                     "1.0,2.0,3.0,4.0,5.0,6.0,7.0\n" *
                     "8.0,9.0,10.0,11.0,12.0,13.0,14.0");

julia> samples = read_samples(io);

julia> samples.a
2-element Array{Float64,1}:
 1.0
 8.0

julia> samples.b
2×2 ElasticArrays.ElasticArray{Float64,2,1}:
 2.0   9.0
 3.0  10.0

julia> samples.c
2×2×2 ElasticArrays.ElasticArray{Float64,3,2}:
[:, :, 1] =
 4.0  6.0
 5.0  7.0

[:, :, 2] =
 11.0  13.0
 12.0  14.0
```
"""
read_samples(filename::AbstractString) = open(read_samples, filename, "r")

function _read_sample_matrix!(samples, io)
    n = size(samples, 1)
    while !eof(io)
        line = readline(io)
        iscommentline(line) && continue
        sample = parse.(Float64, fields(line))
        @argcheck length(sample) == n "Expected $(n) fields in $(line)"
        append!(samples, sample)
    end
    samples
end

function read_sample_matrix(io::IO)
    while !eof(io)
        line = readline(io)
        iscommentline(line) && continue
        header = fields(line)
        samples = ElasticArray{Float64}(undef, length(header), 0)
        return header, collect(PermutedDimsArray(_read_sample_matrix!(samples, io), (2, 1)))
    end
end

"""
$(SIGNATURES)

Read Stan samples from a CSV file or a `IO` stream. Return the header (a vector of strings,
containing variables names) and the sample as an `AbstractMatrix{Float64}`

```jldoctest
julia> io = IOBuffer("a,b.1,b.2,c.1.1,c.2.1,c.1.2,c.2.2\n" *
       "1.0,2.0,3.0,4.0,5.0,6.0,7.0\n" *
       "8.0,9.0,10.0,11.0,12.0,13.0,14.0");

julia> header, matrix = read_sample_matrix(io);

julia> header
7-element Vector{SubString{String}}:
 "a"
 "b.1"
 "b.2"
 "c.1.1"
 "c.2.1"
 "c.1.2"
 "c.2.2"

julia> matrix
2×7 Matrix{Float64}:
 1.0  2.0   3.0   4.0   5.0   6.0   7.0
 8.0  9.0  10.0  11.0  12.0  13.0  14.0
"""
read_sample_matrix(filename::AbstractString) = open(read_sample_matrix, filename, "r")

end # module


text:
_codegen(s::T) where T <: Number = s

_codegen(s::Symbol) = s
_codegen(ex::Expr) = ex

function _codegen(s::Symbolic{T}) where {T}
    f = operation(s)
    args = _codegen.(arguments(s))
    _codegen(T, f, args...)
end

function _codegen(s::Sym{T}) where {T}
    s.name
end

function _codegen(::Type{T}, f::Function, args...) where {T}
    fsym = gensym(Symbol(f))
    argsyms = Vector{Any}(gensym.(Symbol.(:arg_, 1:length(args))))
    ex = @q begin end

    for (j, (arg, argsym)) in enumerate(zip(args, argsyms))
        t = _codegen(arg)

        # if t isa Expr
        #     push!(ex.args, :($argsym = $t))
        # else
            argsyms[j] = t
        # end

    end

    # push!(ex.args, :($fsym = $f($(argsyms...))))
    push!(ex.args, :($f($(argsyms...))))
    ex
end


text:
module SymbolicDiff

export AbstractSymbolic, AbstractVectorSymbolic, AbstractMatrixSymbolic, SymbolicVariable, SymbolicValue, SymbolicExpression
export SymbolicEnv, SymbolicCache
export symbolic, seval
export @vars, @bind, @expr
export globalenv

import Base
import LinearAlgebra: dot
import SparseMatrix: SparseCSR, SparseCSC, SparseCOO
import SparseArrays: SparseMatrixCSC

include("_symbolic.jl")

include("_env.jl")
include("_eval.jl")
include("_deriv.jl")
include("_deriv2.jl")

include("_vector.jl")
include("_matrix.jl")

include("_operations.jl")
include("_macro.jl")

end


text:
#EDAData
#-------------------------------------------------------------------------------
module EDAData

const defaultPSFReader = :LibPSF #Overwrite with LibPSFC to use alternative algorithm

#=
Wraps lower-level EDA file readers in a more appropriate API.

 - Reads results into DataF1 structures.
 - TODO: Support results of parametric sweeps into MDDatasets.DataRS structures.
=#

using MDDatasets
using CMDimCircuits.NetwAnalysis
import SpiceData #SPICE data reader

#Import PSF reader implementation as "PSFReaderLib" (Hacky, but works):
eval(:(using $defaultPSFReader))
eval(:(const PSFReaderLib = $defaultPSFReader))

import CMDimCircuits.NetwAnalysis: NetworkParameters, NetworkParametersRef, NetworkParametersNoRef
import CMDimCircuits.DSVFiles: DSVReader
import Dates: now

include("base.jl")
include("tr0.jl")
include("psf.jl")
include("snp_reader.jl")
include("snp_writer.jl")
include("show.jl")

#==Exported symbols
===============================================================================#

#==Un-exported interface
================================================================================
	open_tr0(path::String) #Opens Tr0Reader
	open_psf(path::String) #Opens PSFReader

	read_sNp(path::String, numports=X) #Returns NetworkParameters matrix
	read_s2p(path::String) #Returns NetworkParameters matrix
==#

#==Extensions to other modules
================================================================================
	Base.read(::{Tr0Reader/PSFReader}, signame::String)
	Base.names(::{Tr0Reader/PSFReader}) #Returns list of signal names
	Base.close(::{Tr0Reader/PSFReader})
==#

end #EDAData


text:
using CUDA.CUTENSOR
using CUDA
using LinearAlgebra

# using host memory with CUTENSOR doesn't work on Windows
can_pin = !Sys.iswindows()

eltypes = (#(Float16, Float16), #(Float16, Float32),
                        (Float32, Float32), #(Float32, Float64),
                        (Float64, Float64),
                        #(ComplexF16, ComplexF16), (ComplexF16, ComplexF32),
                        (ComplexF32, ComplexF32), #(ComplexF32, ComplexF64),
                        (ComplexF64, ComplexF64))
@testset for NA=2:5, NC = 1:NA-1
    @testset for (eltyA, eltyC) in eltypes
        # setup
        eltyD = eltyC
        dmax = 2^div(18,NA)
        dims = rand(2:dmax, NA)
        p = randperm(NA)
        indsA = collect(('a':'z')[1:NA])
        indsC = indsA[p][1:NC]
        dimsA = dims
        dimsC = dims[p][1:NC]
        A = rand(eltyA, (dimsA...,))
        dA = CuArray(A)
        C = rand(eltyC, (dimsC...,))
        dC = CuArray(C)
        # setup
        can_pin && Mem.pin(A)

        opA = CUTENSOR.CUTENSOR_OP_IDENTITY
        opC = CUTENSOR.CUTENSOR_OP_IDENTITY
        opReduce = CUTENSOR.CUTENSOR_OP_ADD
        # simple case
        dC = CUTENSOR.reduction!(1, dA, indsA, opA, 0, dC, indsC, opC, opReduce)
        C = collect(dC)
        @test reshape(C, (dimsC..., ones(Int,NA-NC)...)) ≈
            sum(permutedims(A, p); dims = ((NC+1:NA)...,))
        if can_pin
            Csimple = zeros(eltyC, dimsC...)
            Mem.pin(Csimple)
            Csimple = CUDA.@sync CUTENSOR.reduction!(1, A, indsA, opA, 0, Csimple, indsC, opC, opReduce)
            @test reshape(Csimple, (dimsC..., ones(Int,NA-NC)...)) ≈
                sum(permutedims(A, p); dims = ((NC+1:NA)...,))
        end

        # using integers as indices
        dC = CUTENSOR.reduction!(1, dA, collect(1:NA), opA, 0, dC, p[1:NC], opC, opReduce)
        C = collect(dC)
        @test reshape(C, (dimsC..., ones(Int,NA-NC)...)) ≈
            sum(permutedims(A, p); dims = ((NC+1:NA)...,))
        if can_pin
            Cinteger = zeros(eltyC, dimsC...)
            Mem.pin(Cinteger)
            Cinteger = CUDA.@sync CUTENSOR.reduction!(1, A, collect(1:NA), opA, 0, Cinteger, p[1:NC], opC, opReduce)
            @test reshape(Cinteger, (dimsC..., ones(Int,NA-NC)...)) ≈
                sum(permutedims(A, p); dims = ((NC+1:NA)...,))
        end

        # multiplication as reduction operator
        opReduce = CUTENSOR.CUTENSOR_OP_MUL
        dC = CUTENSOR.reduction!(1, dA, indsA, opA, 0, dC, indsC, opC, opReduce)
        C = collect(dC)
        @test reshape(C, (dimsC..., ones(Int,NA-NC)...)) ≈
            prod(permutedims(A, p); dims = ((NC+1:NA)...,)) atol=eps(Float16) rtol=Base.rtoldefault(Float16)
        if can_pin
            Cmult = zeros(eltyC, dimsC...)
            Mem.pin(Cmult)
            Cmult = CUDA.@sync CUTENSOR.reduction!(1, A, indsA, opA, 0, Cmult, indsC, opC, opReduce)
            @test reshape(Cmult, (dimsC..., ones(Int,NA-NC)...)) ≈
                prod(permutedims(A, p); dims = ((NC+1:NA)...,)) atol=eps(Float16) rtol=Base.rtoldefault(Float16)
            # NOTE: this test often yields values close to 0 that do not compare approximately
        end

        # with non-trivial coefficients and conjugation
        opA = eltyA <: Complex ? CUTENSOR.CUTENSOR_OP_CONJ :
                                CUTENSOR.CUTENSOR_OP_IDENTITY
        opC = CUTENSOR.CUTENSOR_OP_IDENTITY
        opReduce = CUTENSOR.CUTENSOR_OP_ADD
        C = rand(eltyC, (dimsC...,))
        dC = CuArray(C)
        α = rand(eltyC)
        γ = rand(eltyC)
        dC = CUTENSOR.reduction!(α, dA, indsA, opA, γ, dC, indsC, opC, opReduce)
        @test reshape(collect(dC), (dimsC..., ones(Int,NA-NC)...)) ≈
            α .* conj.(sum(permutedims(A, p); dims = ((NC+1:NA)...,))) .+ γ .* C
    end
end


text:
function fieldnames_types(in_type::Type{<:Any})
    return [FieldName_Type(x,fieldtype(in_type, x)) for x in fieldnames(in_type)]
end

function tuple_field_id(in_type::Type{<:Any})
    tmp = filter(x->x.type == DbId, fieldnames_types(in_type))
    result = length(tmp) == 1 ? tmp[1] : nothing
    result === nothing && throw("Please ensure that one field in your struct $in_type has filedtype DbId")
    return result
end

function _get_struct_iid(inp_struct)
    iid_field_tuple = tuple_field_id(typeof(inp_struct))
    iid_field_tuple === nothing && throw("Your type $(typeof(inp_struct)) has no id field or too many.
                                            This field has to be of field type DbId")
    return getproperty(inp_struct, iid_field_tuple.fieldname)
end

max_issues_repo_name: Humans-of-Julia/TypeDBORM

max_forks_repo_name: Humans-of-Julia/TypeDBORM


text:
export DrillReward,
    flow, dflow!, flowdψ

# -------------------------------------------
# abstract types
# -------------------------------------------

# generic functions to access coefs
idx_drill(d) = OneTo(_nparm(d))
theta_drill(d, theta) = view(theta, idx_drill(d))

# @deprecate length(f::AbstractPayoffFunction) _nparm(f)

_nparm(m::AbstractDrillModel) = _nparm(reward(m))
coefnames(x::AbstractDrillModel) = coefnames(reward(x))

# -------------------------------------------
# Generic functions
# -------------------------------------------

NotDefinedError(m) = throw(error("Not defined for $m"))

# navigate through state space
num_choices(m::AbstractStateSpace, i...) = NotDefinedError(m)
actionspace(m::AbstractStateSpace, i...) = NotDefinedError(m)
_Dgt0(      m::AbstractStateSpace, i...) = NotDefinedError(m)

# navigate through model primitives
reward(    m::AbstractDrillModel) = NotDefinedError(m)
statespace(m::AbstractDrillModel) = NotDefinedError(m)

# From model to things about states
@inline actionspace(m::AbstractDrillModel, i...) = actionspace(statespace(m), i...)
@inline num_choices(m::AbstractDrillModel, i...) = num_choices(statespace(m), i...)
@inline _Dgt0(      m::AbstractDrillModel, i...) = _Dgt0(statespace(m), i...)

@inline _ψ(    m::AbstractDrillModel, state, s::SimulationDraw) = _Dgt0(m, state) ? _ψ2(s) : _ψ1(s)
@inline _dψdθρ(m::AbstractDrillModel, state, s::SimulationDraw) = _Dgt0(m, state) ? azero(s) : _dψ1dθρ(s)

# -------------------------------------------
# some functions to look at stuff
# -------------------------------------------

@inline actionspace(obs::ObservationDrill) = actionspace(_model(obs), _x(obs))
@inline statespace( obs::ObservationDrill) = statespace(_model(obs))
@inline reward(     obs::ObservationDrill) = reward(_model(obs))

@inline _sgnext(            obs::ObservationDrill) = _sgnext(statespace(obs), _x(obs), _y(obs))
@inline _sgnext(d::Integer, obs::ObservationDrill) = _sgnext(statespace(obs), _x(obs), d)
@inline _Dgt0(obs::ObservationDrill) = _Dgt0(statespace(obs), _x(obs))
@inline _d(   obs::ObservationDrill) = _y(obs)

@inline _ψ(    obs::ObservationDrill, ψ::Float64) = ψ
@inline _ψ(    obs::ObservationDrill, s::SimulationDraw) = _ψ(_model(obs), _x(obs), s)
@inline _dψdθρ(obs::ObservationDrill, s::SimulationDraw) = _dψdθρ(_model(obs), _x(obs), s)

@inline expires_today(            obs::ObservationDrill) = expires_today(statespace(obs), _x(obs), _y(obs))
@inline expires_today(d::Integer, obs::ObservationDrill) = expires_today(statespace(obs), _x(obs), d)

# -------------------------------------------
# Payoffs...
# -------------------------------------------

# @inline flow(        d, obs, theta, s) = flow(  reward(_model(obs)),           d,   obs, theta, s)
# @inline dflow!(grad, d, obs, theta, s) = dflow!(reward(_model(obs)), grad,     d,   obs, theta, s)
# @inline flowdψ(grad, d, obs, theta, s) = flowdψ(reward(_model(obs)), grad,     d,   obs, theta, s)
# @inline flow(           obs, theta, s) = flow(  reward(_model(obs)),       _y(obs), obs, theta, s)
# @inline dflow!(grad,    obs, theta, s) = dflow!(reward(_model(obs)), grad, _y(obs), obs, theta, s)
# @inline flowdψ(grad,    obs, theta, s) = flowdψ(reward(_model(obs)), grad, _y(obs), obs, theta, s)
# @deprecate dflowdψ(args...) flowdψ(args...)

@inline flow!(grad, obs, θ, s, dograd) = flow!(grad, reward(_model(obs)), _y(obs), obs, θ, s, dograd)

function dflow(x::AbstractPayoffFunction, d, obs, theta, s)
    length(theta) == _nparm(x) || throw(DimensionMismatch())
    grad = zero(theta)
    dflow!(grad, x, d, obs, theta, s, true)
    return grad
end

function flow(x::AbstractPayoffFunction, d, obs, theta, s)
    grad = Vector{eltype(theta)}(undef, 0)
    return flow!(grad, x, d, obs, theta, s, false)
end

@inline function full_payoff!(grad, d, obs, theta, sim, dograd)
    return flow!(grad, reward(_model(obs)), d, obs, theta, sim, dograd)
end

@inline function full_payoff(d, obs, theta, sim)
    grad = Vector{eltype(theta)}(undef, 0)
    full_payoff!(grad, d, obs, theta, sim, false)
end

# Check Derivatives
#-------------------------------------------------------------

"check gradient of reward function"
function check_flow_grad(m, d, obs, theta, sim)

    # compute analytic
    grad = ones(eltype(theta), length(theta))
    flow!(grad, m, d, obs, theta, sim, true)

    # check finite difference for dθ
    f(x) = flow!(grad, m, d, obs, x, sim, false)
    @test grad ≈ Calculus.gradient(f, theta)

    # check finite difference for dψ
    fdpsi(x) = flow!(grad, m, d, obs, theta, SimulationDrawFD(sim, x), false)
    @test flowdψ(m, d, obs, theta, sim) ≈ Calculus.derivative(fdpsi, 0.0)
end

function check_flow_grad(m, d, obs, theta)
    u, v = randn(2)
    simdraw(x) = SimulationDraw(u, v, theta_drill_ρ(m, x))
    sim = simdraw(theta)
    check_flow_grad(m,d,obs,theta,sim)
end

# @deprecate flow(x::AbstractDrillModel) reward(x)


text:
module DerivableFunctionsBase

using ForwardDiff
using OffsetArrays # To avoid error thrown in _array_for() when using Symbolics.jacobian()
using Symbolics

# Add graceful errors by implementing _GetGrad(::Val) methods

include("Utils.jl")
export GetArgLength

include("DFunctions.jl")
export DFunction, DerivableFunction
export EvalF, EvaldF, EvalddF, In, Out, InOut

include("DifferentiationOperators.jl")
export diff_backends
export GetDeriv, GetGrad, GetJac, GetHess, GetMatrixJac, GetDoubleJac
export GetGrad!, GetJac!, GetHess!, GetMatrixJac!

end


text:
# NaturalSelection.jl
# ===================
#
# A julia package that provides methods for detecting the presence, strength
# and effects of natural selection, in biological data.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

__precompile__()

module NaturalSelection

export
    dNdS_NG86,
    S_N_NG86,
    DS_DN_NG86,
    tajimad,
    mkt,
    make_edge_reference

import BioSequences:
    BioSequences,
    DNA,
    BioSequence,
    DNASequence,
    NucleicAcid,
    Kmer,
    NucAlphs,
    GeneticCode,
    ispurine,
    @kmer_str,
    ncbi_trans_table,
    standard_genetic_code

import GeneticVariation:
    NL79,
    Segregating,
    avg_mut

using PairwiseListMatrices

const DEFAULT_TRANS = BioSequences.ncbi_trans_table[1]
const Codon{T} = BioSequences.Kmer{T, 3}

@inline bitindex(x::Kmer{T,K}, i::Integer) where {T,K} = 2 * (K - i)

include("codons/codon_set.jl")

include("dnds/dnds.jl")
include("tajima.jl")
end


text:
function get_nv(clique)
    length(unique(vcat(clique...)))
end

function get_ne(clique)
    edges = Set([])
    for clique in clique
        for edge in collect(combinations(clique, 2))
            push!(edges, Set(edge))
        end
    end
    return length(edges)
end

function get_nb_lc(clique, cliquetree)
    dstmx = zeros(Int, length(clique), length(clique))
    for i in 1:length(clique) - 1
        for j in i + 1:length(clique)
            value = length(intersect(clique[i], clique[j]))
            dstmx[i, j] = value
            dstmx[j, i] = value
        end
    end
    nb_lc = sum([dstmx[src, dst] * (2 * dstmx[src, dst] + 1) for (src, dst) in cliquetree])
    return nb_lc
end

save_clique(clique::Vector{Vector{Int}}, path::AbstractString) = writedlm(path, clique)

save_cliquetree(cliquetree::Vector{Vector{Int}}, path::AbstractString) = writedlm(path, cliquetree)

function read_cliquetree(path)
    cliquetree = readdlm(path, '\t', Int)
    cliquetree = [cliquetree[i, :] for i in 1:size(cliquetree, 1)]
    return cliquetree
end

function read_clique(path)
    lines = split(read(open(path, "r"), String), '\n')[begin:end-1]
    clique = [parse.(Int, split(line, "\t")) for line in lines]
    return clique
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
# helper function for constrain
function _colconstraints!(col::Vector{Union{T,Null}}, bfunc::Function, keep::Vector{Bool}) where T
    for i ∈ 1:length(keep)
        if !isnull(col[i])
            keep[i] &= bfunc(get(col[i]))
        else
            keep[i] = false
        end
    end
    keep
end

# this is for fixing slowness due to bad dispatching
# performance is better, but it's still slow
function _dispatchConstrainFunc!(f::Function, mask::Union{Vector{Bool},BitArray},
                                 keep::BitArray, cols::Vector{Union{T,Null}}...) where T
    ncols = length(cols)
    # for some reason this completely fixes the performance issues
    # still don't completely understand why
    get_idx(i, j) = get(cols[j][i])
    get_args(i) = (get_idx(i, j) for j ∈ 1:ncols)
    # all arg cols are same length
    for i ∈ 1:length(keep)
        if !mask[i]
            keep[i] = false
        else
            # this is inexplicably slow
            # note that this is even true for the generator alone; not a dispatch issue
            # keep[i] = f((get(col[i]) for col ∈ cols)...)
            keep[i] = f(get_args(i)...)
        end
    end
end

"""
    constrain(df, dict)
    constrain(df, kwargs...)
    constrain(df, cols, func)

Returns a subset of the dataframe `df` for which the column `key` satisfies
`value(df[i, key]) == true`.  Where `(key, value)` are the pairs in `dict`.
Alternatively one can use keyword arguments instead of a `Dict`.

Also, one can pass a function the arguments of which are elements of columns specified
by `cols`.
"""
function constrain(df::AbstractDataFrame,
                   constraints::Dict{K, V})::DataFrame where {K<:Symbol,V<:Function}
    keep = ones(Bool, size(df, 1))
    for (col, bfunc) ∈ constraints
        _colconstraints!(df[col], bfunc, keep)
    end
    df[keep, :]
end

function constrain{K, V<:Array}(df::AbstractDataFrame, constraints::Dict{K, V})::DataFrame
    newdict = Dict(k=>(x -> x ∈ v) for (k, v) ∈ constraints)
    constrain(df, newdict)
end

constrain(df::AbstractDataFrame; kwargs...) = constrain(df, Dict(kwargs))

function constrain(df::AbstractDataFrame, cols::Vector{Symbol}, f::Function)
    keep = BitArray(size(df, 1))
    _dispatchConstrainFunc!(f, completecases(df[cols]), keep, (df[col] for col ∈ cols)...)
    df[keep, :]
end
export constrain

# this is a helper function to @constrain
# TODO fix this so it works for variables used multiple times
function _checkConstraintExpr!(expr::Expr, dict::Dict)
    # the dictionary keys are the column names (as :(:col)) and the values are the symbols
    for (idx, arg) ∈ enumerate(expr.args)
        if isa(arg, QuoteNode)
            newsym = gensym()
            dict[Meta.quot(arg.value)] = newsym
            expr.args[idx] = newsym
        elseif  isa(arg, Expr) && arg.head == :quote
            newsym = gensym()
            dict[Meta.quot(arg.args[1])] = newsym
            expr.args[idx] = newsym
        elseif isa(arg, Expr)
            _checkConstraintExpr!(expr.args[idx], dict)
        end
    end
end

# TODO this is still having performance issues
"""
    @constrain(df, expr)

Constrains the dataframe to rows for which `expr` evaluates to `true`.  `expr` should
specify columns with column names written as symbols.  For example, to do `(a ∈ A) > M`
one should write `:A .> M`.
"""
macro constrain(df, expr)
    dict = Dict()
    _checkConstraintExpr!(expr, dict)
    cols = collect(keys(dict))
    vars = [dict[k] for k ∈ cols]
    cols_expr = Expr(:vect, cols...)
    fun_name = gensym()
    o = quote
        function $fun_name($(vars...))
            $expr
        end
        constrain($df, $cols_expr, $fun_name)
    end
    esc(o)
end
export @constrain


text:
import Term: load_code_and_highlight, highlight_syntax, highlight, theme

@testset "\e[34mHIGHLIGHT" begin

    # @test highlight("test 1 123 33.4 44,5 +1 -2 12 0.5, ,, ...") == "test[#90CAF9][#90CAF9] 1 [/#90CAF9][/#90CAF9][#90CAF9]123[/#90CAF9] [#90CAF9]33.4[/#90CAF9] [#90CAF9]44,5[/#90CAF9] [#90CAF9]+1 [/#90CAF9][#90CAF9]-2 [/#90CAF9][#90CAF9]12[/#90CAF9] [#90CAF9]0.5[/#90CAF9][#90CAF9], [/#90CAF9],[#90CAF9], [/#90CAF9]..."

    @test highlight("this is ::Int64") == "this is [#CE93D8]::Int64[/#CE93D8]"

    @test highlight("print", :func) == "[#FFEE58]print[/#FFEE58]"

    @test highlight("1 + 2", :code) == "[#90CAF9]1 + 2[/#90CAF9]"

    @test highlight("this 1 + `test`") == "this 1 + [#90CAF9]`test`[/#90CAF9]"

    @test highlight(1) == "[#42A5F5]1[/#42A5F5]"

    @test highlight([1, 2, 3]) == "[#42A5F5][1, 2, 3][/#42A5F5]"

    @test highlight(Int32) == "[#CE93D8]Int32[/#CE93D8]"

    @test highlight(print) == "[#FFEE58]print[/#FFEE58]"

    if !Sys.iswindows()

        @test highlight_syntax("""
        This is ::Int64 my style
        """) == "This is \e[1m\e[38;2;252;98;98m::\e[22m\e[39mInt64 my style\n"

        @test highlight_syntax("""
        This is ::Int64 my style
        print(x + 2)
        """) == "This is \e[1m\e[38;2;252;98;98m::\e[22m\e[39mInt64 my style\n\e[1m\e[38;2;255;245;157mprint\e[22m\e[39m\e[38;2;252;116;116m(\e[39mx \e[1m\e[38;2;252;98;98m+\e[22m\e[39m \e[38;2;144;202;249m2\e[39m\e[38;2;252;116;116m)\e[39m\n"

        @test load_code_and_highlight("14_test_highlight.jl", 7) == "[red bold]❯[/red bold] [white]7[/white] \n  [grey39]8[/grey39]  @test [(255, 245, 157)  bold  ]highlight[/(255, 245, 157)  bold  ][(252, 116, 116)    ]([/(252, 116, 116)    ][(165, 214, 167)    ]\"this is ::Int64\"[/(165, 214, 167)    ][(252, 116, 116)    ])[/(252, 116, 116)    ] [(252, 98, 98)  bold  ]==[/(252, 98, 98)  bold  ] [(165, 214, 167)    ]\"this is [[#CE93D8]]::Int64[[/#CE93D8]]\"[/(165, 214, 167)    ]\n  [grey39]9[/grey39] \n  [grey39]10[/grey39]  @test [(255, 245, 157)  bold  ]highlight[/(255, 245, 157)  bold  ][(252, 116, 116)    ]([/(252, 116, 116)    ][(165, 214, 167)    ]\"print\"[/(165, 214, 167)    ][(252, 116, 116)    ],[/(252, 116, 116)    ] [(165, 214, 167)    ]:func[/(165, 214, 167)    ][(252, 116, 116)    ])[/(252, 116, 116)    ] [(252, 98, 98)  bold  ]==[/(252, 98, 98)  bold  ] [(165, 214, 167)    ]\"[[#FFEE58]]print[[/#FFEE58]]\"[/(165, 214, 167)    ]\n  [grey39]11[/grey39] \n  [grey39]12[/grey39]  @test [(255, 245, 157)  bold  ]highlight[/(255, 245, 157)  bold  ][(252, 116, 116)    ]([/(252, 116, 116)    ][(165, 214, 167)    ]\"1 + 2\"[/(165, 214, 167)    ][(252, 116, 116)    ],[/(252, 116, 116)    ] [(165, 214, 167)    ]:code[/(165, 214, 167)    ][(252, 116, 116)    ])[/(252, 116, 116)    ] [(252, 98, 98)  bold  ]==[/(252, 98, 98)  bold  ] [(165, 214, 167)    ]\"[[#90CAF9]]1 + 2[[/#90CAF9]]\"[/(165, 214, 167)    ]\n  [grey39]13[/grey39] "

    end
end


text:
abstract type BatchDim end

struct BatchDimFirst <: BatchDim end
struct BatchDimLast <: BatchDim end

"""
    BatchViewCollated(data, size; droplast = false)

A batch view of container `data` with collated batches of
size `size`.
"""
@with_kw_noshow struct BatchViewCollated{TData}
    data::TData
    size::Int
    count::Int
    partial::Bool
    batchdim::BatchDim = BatchDimLast()
end

Base.show(io::IO, bvcollated::BatchViewCollated) = print(io, "batchviewcollated() with $(bvcollated.count) batches of size $(bvcollated.size)")

const batchviewcollated = BatchViewCollated

function BatchViewCollated(data, size; partial = true, batchdim = BatchDimLast())
    count = nobs(data) ÷ size
    if partial && (nobs(data) % size > 0)
        count += 1
    end
    return BatchViewCollated(data, size, count, partial, batchdim)
end

Base.length(bv::BatchViewCollated) = bv.count
LearnBase.nobs(bv::BatchViewCollated) = length(bv)

function LearnBase.getobs(bv::BatchViewCollated, idx::Int)
    idxs = batchindices(nobs(bv.data), bv.size, idx)
    collate([getobs(bv.data, idx) for idx in idxs])
end

function LearnBase.getobs!(buf, bv::BatchViewCollated, idx::Int)
    indices = batchindices(nobs(bv.data), bv.size, idx)
    for (idx, obs) in zip(indices, obsslices(buf, bv.batchdim))
        obs_ = getobs!(obs, bv.data, idx)
        # if data container does not implement, getobs!, this is needed
        if obs_ !== obs
            copyrec!(obs, obs_)
        end
    end

    if bv.partial && (idx == nobs(bv)) && (nobs(bv.data) % bv.size > 0)
        # This will mess up the buffer in the `RingBuffer`
        # by not reinserting the correct buffer.
        # It is okay only because it will that buffer will
        # not be used since the partial batch is the last
        # batch.
        return obsslice(buf, 1:(nobs(bv.data) % bv.size), bv.batchdim)
    end

    return buf
end

# batch view helpers

"""
    obsslices(batch, batchdim = BatchDimLast())

Iterate over views of all observations in a `batch`.
`batch` can be a batched array, a tuple of batches, or a
dict of batches.

```julia
batch = rand(10, 10, 4)  # batch size is 4
iter = obsslices(batch, BatchDimLast())
@assert size(first(iter)) == (10, 10)

iter2 = obsslices(batch, BatchDimFirst())
@assert size(first(iter)) == (10, 4)
```
"""
obsslices(batch, batchdim = BatchDimLast()) =
    (obsslice(batch, i, batchdim) for i in 1:_batchsize(batch, batchdim))

function obsslice(batch::AbstractArray{T, N}, i, ::BatchDimLast) where {T, N}
    return view(batch, [(:) for _ in 1:N-1]..., i)
end

function obsslice(batch::AbstractArray{T, N}, i, ::BatchDimFirst) where {T, N}
    return view(batch, i, [(:) for _ in 2:N]...)
end

function obsslice(batch::Tuple, i, batchdim)
    return Tuple(obsslice(batch[j], i, batchdim) for j in 1:length(batch))
end

function obsslice(batch::NamedTuple, i, batchdim)
    return (; zip(keys(batch), obsslice(values(batch), i, batchdim))...)
end

function obsslice(batch::Dict, i, batchdim)
    return Dict(k => obsslice(v, i, batchdim) for (k, v) in batch)
end

# Utils

_batchsize(batch::Tuple, batchdim) = _batchsize(batch[1], batchdim)
_batchsize(batch::Dict, batchdim) = _batchsize(batch[first(keys(batch))], batchdim)
_batchsize(batch::AbstractArray{T, N}, ::BatchDimLast) where {T, N} = size(batch, N)
_batchsize(batch::AbstractArray, ::BatchDimFirst) = size(batch, 1)

copyrec!(dst::AbstractArray, src::AbstractArray) = copy!(dst, src)
copyrec!(dst::Tuple, src::Tuple) = foreach((a, b) -> copyrec!(a, b), dst, src)
copyrec!(dst::NamedTuple, src::NamedTuple) = copyrec!(values(dst), values(src))
copyrec!(dst::Dict, src::Dict) = foreach((a, b) -> copyrec!(dst, src), values(dst), values(src))

"""
    batchindices(n, size, i)

Get the indices of batch `i` with batch size `size` of
a collection with `n` elements.

Might be a partial batch if `i` is the last batch and
`n` is not divisible by `size`.
"""
function batchindices(n, size::Int, i::Int)
    from = (i - 1) * size + 1
    to = min(n, from + size - 1)
    return from:to
end


text:
"""
Helper functions for working with frontend assets (including JS, CSS, etc files).
"""
module Assets

import Genie, Genie.Configuration

export include_asset, css_asset, js_asset

### PUBLIC ###

"""
    include_asset(asset_type::Union{String,Symbol}, file_name::Union{String,Symbol};
                  fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String

Returns the path to an asset. `asset_type` can be one of `:js`, `:css`. The `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicating if a fingerprint (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function include_asset(asset_type::Union{String,Symbol}, file_name::Union{String,Symbol};
                        fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  asset_type = string(asset_type)
  file_name = string(file_name)

  suffix = fingerprinted ? "-" * Genie.ASSET_FINGERPRINT * ".$asset_type" : ".$asset_type"
  "/$asset_type/$(file_name)$(suffix)"
end

"""
    css_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String

Path to a css asset. The `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicating if a fingerprint (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function css_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  include_asset(:css, file_name, fingerprinted = fingerprinted)
end

"""
    js_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String

Path to a js asset. `file_name` should not include the extension.
`fingerprinted` is a `Bool` indicating if a fingerprint (unique hash) should be added to the asset's filename (used in production to invalidate caches).
"""
function js_asset(file_name::String; fingerprinted::Bool = Genie.config.assets_fingerprinted) :: String
  include_asset(:js, file_name, fingerprinted = fingerprinted)
end

end


text:
@testset "Running Source Tests" begin

    # ------------------------- #
    # --- Run source tests  --- #
    # ------------------------- #

    # --- Move immediate --- #
    core = test_core(TestPrograms.move_immediate())
    # Clock core.
    runfor(core, 10)
    # Test that the immediate was moved into the final dmem location.
    # Note that we must convert to index 1 when accessing the dmem array in Core.
    @test core.dmem[1] == 100

    # --- Move Return Address --- #
    core = test_core(TestPrograms.move_return_address())
    runfor(core, 10)
    @test core.dmem[1] == TestPrograms.RETURN_ADDRESS

    # --- Move Accumulator --- #
    core =  test_core(TestPrograms.move_accumulator())
    runfor(core, 10)
    # Awkward conversion - take lower 16 bits of accumulator, convert to UInt16,
    # then reinterpret memory as Int16 to match the data types in DMEM.
    @test core.dmem[1] == reinterpret(Int16, UInt16(TestPrograms.ACCUMULATOR & 0xFFFF))

    # -- Bypass Registers -- #
    bypass_programs = (
        TestPrograms.bypass_1(),
        TestPrograms.bypass_2(),
        TestPrograms.bypass_3(),
    )
    for (index, program) in enumerate(bypass_programs)
        println("Testing Bypass Register: $index")
        core = test_core(program)
        runfor(core, 10)
        @test core.dmem[1] == 10
    end

    # Bypass 4
    println("Testing Bypass Register: 4")
    core = test_core(TestPrograms.bypass_4())
    runfor(core, 15)
    @test core.dmem[1] == 100
    @test core.dmem[2] == 100 * 10
    @test core.dmem[3] == 2 * 100 * 10

    # Bypass 5
    println("Testing Bypass Register: 5")
    core = test_core(TestPrograms.bypass_5())
    runfor(core, 15)
    @test core.dmem[1] == 100
    @test core.dmem[2] == 100 * 10
    @test core.dmem[3] == 2 * 100 * 10

    # -- Test pointer dereference -- #
    core = test_core(TestPrograms.pointers())
    runfor(core, 12)
    @test core.dmem[1] == TestPrograms.DMEM[128 + 1]
    @test core.dmem[2] == TestPrograms.DMEM[129 + 1]
    @test core.dmem[3] == TestPrograms.DMEM[130 + 1]
    @test core.dmem[4] == TestPrograms.DMEM[131 + 1]

    # -- Address Generator --
    core = test_core(TestPrograms.address_generators())
    runfor(core, 15)
    # Expect this to be the contents of DMEM at address 12 + 1.
    @test core.dmem[0 + 1] == TestPrograms.DMEM[12 + 1]
    # Since the address gneerator was not incremented, still expect it to have
    # dereferenced the same memory location.
    @test core.dmem[1 + 1] == TestPrograms.DMEM[12 + 1]
    # The previous instruction incremented the address generator. Since we used
    # an ADD instruction, we expect the result to be twice the value in memory
    # at address 14 + 1
    @test core.dmem[2 + 1] == 2 * TestPrograms.DMEM[14 + 1]
    # Since the ADD incremented the address generator, expect the result to be
    # the contents at DMEM 16 + 1
    @test core.dmem[3 + 1] == TestPrograms.DMEM[16 + 1]

    # -- Test Pointer Bypass -- #
    core = test_core(TestPrograms.pointer_bypass())
    runfor(core, 10)
    @test core.dmem[1] == TestPrograms.DMEM[12 + 1]

    # -- Run a fifo test -- #
    core = test_core(TestPrograms.fifo_read())
    runfor(core, 10)
    # Nothing should have been written back to memory, core should be stalled
    # on FIFO read.
    @test core.dmem[1] == 0

    fifo_value_1 = -1
    fifo_value_2 = -2
    # Write to the fifo.
    fifo = core.fifos[1]
    write(fifo, fifo_value_1)
    AsapSim.writeupdate!(fifo)
    # Run the core some more. The MOVI instruction should be stalled in stage 5
    # before wite back
    runfor(core, 4)
    @test core.dmem[1] == 0
    @test core.dmem[2] == 0
    # Write again to the fifo.
    write(fifo, fifo_value_2)
    AsapSim.writeupdate!(fifo)
    runfor(core, 10)

    @test core.dmem[1] == 10
    @test core.dmem[2] == fifo_value_1
    @test core.dmem[3] == fifo_value_2
    @test core.dmem[4] == fifo_value_2
end

################################################################################
# Branch Tests
################################################################################
@testset "Running Branch Tests" begin
    # Test branch predicted taken
    core = test_core(TestPrograms.branch_forward_taken_1())
    runfor(core, 20)
    # Make sure first MOVE worked
    @test core.dmem[1] == 1
    # These instructions should have been skipped.
    @test core.dmem[2] == 0
    @test core.dmem[3] == 0
    @test core.dmem[4] == 0
    @test core.dmem[5] == 0
    # This instruction should have been branced to and executed
    @test core.dmem[6] == 10

    # --
    core = test_core(TestPrograms.branch_forward_taken_2())
    runfor(core, 15)
    @test core.dmem[1] == 1
    @test core.dmem[2] == 1

    # -- Test branch predicted not-taken
    core = test_core(TestPrograms.branch_forward_untaken_1())
    runfor(core, 20)
    # Make sure first MOVE worked
    @test core.dmem[1] == 1
    # These instructions should have been skipped.
    @test core.dmem[2] == 0
    @test core.dmem[3] == 0
    @test core.dmem[4] == 0
    @test core.dmem[5] == 0
    # This instruction should have been branced to and executed
    @test core.dmem[6] == 10

    # --
    core = test_core(TestPrograms.branch_forward_untaken_2())
    runfor(core, 15)
    @test core.dmem[1] == 1
    @test core.dmem[2] == 1

    # ------------------ #
    # Branch Stress Test #
    # ------------------ #
    core = test_core(TestPrograms.branch_stress_1())
    runfor(core, 30)
    # Make sure the entries we want set are set
    for i in 1:5
        @test core.dmem[i] == 1
    end
    # Make sure nothing else was set.
    for i in 50:62
        @test core.dmem[i] == 0
    end

    # -------------------- #
    # Test Jump and Return #
    # -------------------- #
    core = test_core(TestPrograms.test_jump_1())
    runfor(core, 10)
    @test core.return_address == 3

    # Test that return address is handled correctly when called at the end
    # of a RPT block
    core = test_core(TestPrograms.test_jump_2())
    runfor(core, 20)

    # Make sure the return address was saved correctly
    @test core.return_address == 2
    @test core.dmem[1] == 1
    for i in 2:4
        @test core.dmem[i] == 0
    end

    # -- Test the retun instruction
    core = test_core(TestPrograms.test_return_1())
    runfor(core, 30)

    # Expect address 2 to be saved as the return address
    @test core.return_address == 2
    # Should have written 1 to addresses 0 - 2
    for i in 0:2
        @test core.dmem[i+1] == 1
    end
    # Should NOT have written to address 5 - 7
    for i in 5:7
        @test core.dmem[i+1] == 0
    end

end

@testset "Running Repeat Tests" begin
    core = test_core(TestPrograms.repeat_1())
    runfor(core, 20)
    @test core.dmem[1] == 10

    core = test_core(TestPrograms.repeat_2())
    runfor(core, 40)
    @test core.dmem[1] == 10
    @test core.dmem[2] == 10

    core = test_core(TestPrograms.repeat_3())
    @time runfor(core, 60)
    @test core.dmem[1] == 10
    @test core.dmem[2] == 10
    @test core.dmem[3] == 10

    # Test doing a function call.
    core = test_core(TestPrograms.repeat_4())
    @time runfor(core, 80)
    # Should have added 1 to dmem[0] 10 times. Since it should have started at
    # 1, the final value should be 11.
    @test core.dmem[1] == 11
    @test core.dmem[2] == 1

    # Test multiple repeat blocks.
    core = test_core(TestPrograms.repeat_5())
    @time runfor(core, 80)
    @test core.dmem[1] == 10
    @test core.dmem[2] == 10
end

################################################################################
# Destination Write Tests
################################################################################

@testset "Testing Destination Writes" begin
    # Test writing to hardware pointers.
    core = test_core(TestPrograms.pointer_write_1())
    runfor(core, 10)
    @test core.pointers[1] == 1
    @test core.pointers[2] == 2
    @test core.pointers[3] == 3
    @test core.pointers[4] == 4

    # Test writing to address generators
    core = test_core(TestPrograms.address_generator_write_1())
    runfor(core, 20)
    # AG 0
    @test core.address_generators[1].start == 10
    @test core.address_generators[1].stop == 20
    @test core.address_generators[1].stride == 3
    # AG 1
    @test core.address_generators[2].start == 0
    @test core.address_generators[2].stop == 255
    @test core.address_generators[2].stride == 255
    # AG 2
    @test core.address_generators[3].start == 100
    @test core.address_generators[3].stop == 0
    @test core.address_generators[3].stride == 1

    # Write to CX Mask
    core = test_core(TestPrograms.cxmask_write_1())
    runfor(core, 10)
    @test core.condexec[1].mask == signed(UInt16(0xFFFF))
    @test core.condexec[1].unary_op == :XOR
    @test core.condexec[1].early_kill == true

    # Test obuf mask write
    core = test_core(TestPrograms.obufmask_write_1())
    runfor(core, 10)
    for i in core.obuf_mask
        @test i == true
    end
end


text:
using MinFEM

mesh = unit_square(100)
L = asmLaplacian(mesh)
M = asmMassMatrix(mesh)

n=3
m=2
f(x) = ((n*pi)^2 + (m*pi)^2) *sin(n*x[1]*pi)*sin(m*x[2]*pi)
s = evaluateMeshFunction(mesh, f)

boundary = union(mesh.Boundaries[1001].Nodes,
                 mesh.Boundaries[1002].Nodes,
                 mesh.Boundaries[1003].Nodes,
                 mesh.Boundaries[1004].Nodes)

pde = PDESystem(A=L, b=M*s, bc=zeros(mesh.nnodes), DI=boundary)
solve(pde)

vtkfile = open_vtk_file(mesh, "poisson.vtu")
write_point_data(vtkfile, pde.state, "Y")
write_point_data(vtkfile, s, "S")
save_vtk_file(vtkfile)


text:
# PreprocessorOptions
function PrintStats(x::PreprocessorOptions)
    @check_ptrs x
    return clang_PreprocessorOptions_PrintStats(x)
end


text:
#! /usr/bin/julia

# Rosetta Code, Roman numerals/Encode

function romanencode(n::Integer)
    const DR = [["I", "X", "C", "M"] ["V", "L", "D", "MMM"]]
    rnum = ""
    if n > 4999 || n < 1
        throw(DomainError())
    end
    for (omag, d) in enumerate(digits(n))
        if d == 0
            omr = ""
        elseif d < 4
            omr = DR[omag, 1]^d
        elseif d == 4
            omr = DR[omag, 1]*DR[omag, 2]
        elseif d == 5
            omr = DR[omag, 2]
        elseif d < 9
            omr = DR[omag, 2]*(DR[omag, 1]^(d - 5))
        else
            omr = DR[omag, 1]*DR[(omag +1), 1]
        end
        rnum = omr*rnum
    end
    return rnum
end

testcases = Int64[1990, 2008, 1668]
for i in 1:12
    push!(testcases, rand(1:4999))
end
testcases = unique(testcases)

println("Test romanencode, arabic => roman:")
for i in testcases
    println(i, " => ", romanencode(i))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using ACE1pack, Test, LazyArtifacts

##

@testset "ACE1pack.jl" begin

    @testset "Read data" begin include("test_data.jl") end

    @testset "Basis" begin include("test_basis.jl") end

    @testset "Solver" begin include("test_solver.jl") end

    @testset "Fit ACE" begin include("test_fit.jl") end

    @testset "Read params" begin include("test_read_params.jl") end

end


text:
const neib_8 = @SMatrix [1.0 1 1; 1 0 1; 1 1 1]
const neib_8_inf = @SMatrix [1.0 1 1; 1 Inf 1; 1 1 1]

function buffer_array(A::AbstractMatrix{<:Real})
    oA = OffsetMatrix(fill(zero(eltype(A)), size(A) .+ 2), UnitRange.(0, size(A) .+ 1))
    # Update center
    oA[begin+1:end-1, begin+1:end-1] .= A
    # Set edges to mirror center
    oA[begin, begin+1:end-1] .= A[begin, :]
    oA[end, begin+1:end-1] .= A[end, :]
    oA[begin+1:end-1, begin] .= A[:, begin]
    oA[begin+1:end-1, end] .= A[:, end]
    # Set corners to mirror corners of center
    oA[begin, begin] = A[begin, begin]
    oA[begin, end] = A[begin, end]
    oA[end, begin] = A[end, begin]
    oA[end, end] = A[end, end]
    return oA
end

"""
roughness(dem::Matrix{<:Real})

Roughness is the largest inter-cell difference of a central pixel and its surrounding cell, as defined in Wilson et al (2007, Marine Geodesy 30:3-35).
"""
function roughness(dem::AbstractMatrix{<:Real})

    ex_dem = buffer_array(dem)
    roughness = similar(dem)
    x = @MMatrix zeros(3, 3)

    @inbounds for I ∈ CartesianIndices(size(roughness))
        patch = I-Δ:I+Δ
        rdata = view(ex_dem, patch)

        x .= rdata .- rdata[2, 2]
        roughness[I] = maximum(abs.(x))
    end
    roughness
end

"""
TPI(dem::Matrix{<:Real})

TPI stands for Topographic Position Index, which is defined as the difference between a central pixel and the mean of its surrounding cells (see Wilson et al 2007, Marine Geodesy 30:3-35).
"""
function TPI(dem::AbstractMatrix{<:Real})

    ex_dem = buffer_array(dem)
    tpi = similar(dem)
    x = @MMatrix zeros(3, 3)
    Δ = CartesianIndex(1, 1)

    @inbounds for I ∈ CartesianIndices(size(tpi))
        patch = I-Δ:I+Δ
        rdata = view(ex_dem, patch)

        x .= rdata .* neib_8
        tpi[I] = rdata[2, 2] - mean(x)
    end
    return tpi
end

"""
TRI(dem::Matrix{<:Real})

TRI stands for Terrain Ruggedness Index, which measures the difference between a central pixel and its surrounding cells.
This algorithm uses the square root of the sum of the square of the difference between a central pixel and its surrounding cells.
This is recommended for terrestrial use cases.
"""
function TRI(dem::AbstractMatrix{<:Real})

    ex_dem = buffer_array(dem)
    tri = similar(dem)
    x = @MMatrix zeros(3, 3)
    Δ = CartesianIndex(1, 1)

    @inbounds for I ∈ CartesianIndices(size(tri))
        patch = I-Δ:I+Δ
        rdata = view(ex_dem, patch)

        x .= (rdata .- rdata[2, 2]) .^ 2
        tri[I] = sqrt(sum(x))
    end
    return tri
end

"""
pitremoval(dem::Matrix{<:Real})

Remove pits from a DEM Array if the center cell of a 3x3 patch is `limit` lower or than the surrounding cells.
"""
function pitremoval(dem::AbstractMatrix{<:Real}, limit = 2.0)

    ex_dem = buffer_array(dem)
    tri = similar(dem)
    x = @MMatrix zeros(3, 3)
    Δ = CartesianIndex(1, 1)

    @inbounds for I ∈ CartesianIndices(size(tri))
        patch = I-Δ:I+Δ
        x .= view(ex_dem, patch)
        tri[I] = _pitremoval(x, limit)
    end
    return tri
end

@inline function _pitremoval(x, limit)
    A = vec(x)
    order = sortperm(A)
    ifelse(
        (order[1] == 5) &&
            ((A[order[2]] - A[order[1]]) >= limit),
        Inf,
        x[2, 2])
end


text:
using InteractionNetworkModels

d1 = LSP()

x = [1,1,2,1,2,1,2]
y = [1,1,2,1,3,4,5]

@time d1(x,y)

d2 = FastLSP(100)

@time d2(x,y)


text:
@testset " Maxwell 1D FEM solver " begin
    using GEMPIC
    using LinearAlgebra

    """
    Test 1D Maxwell spline finite element solver on a periodic grid

    `L_x` domain dimensions and M is an integer.

    ```math
    B_z(x,y,t) = \\cos(\\frac{2 M \\pi}{L_x} x)  \\cos(\\frac{2 M \\pi}{L_x} t)
    ```
    ```math
    E_y(x,y,t) = \\sin(\\frac{2 M \\pi}{L_x} x)  \\sin(\\frac{2 M \\pi}{L_x} t)
    ```

    """

    mode = 2
    eta1_min = 0.0
    eta1_max = 2π
    nc_eta1 = 256

    Lx = eta1_max - eta1_min

    delta_eta1 = Lx / nc_eta1

    mesh = OneDGrid(eta1_min, eta1_max, nc_eta1)

    deg = 3

    maxwell_1d = Maxwell1DFEM(mesh, deg)

    ex = zeros(Float64, nc_eta1)
    ey = zeros(Float64, nc_eta1)
    bz = zeros(Float64, nc_eta1)

    bz_exact = zeros(Float64, nc_eta1)
    ex_exact = zeros(Float64, nc_eta1)
    ey_exact = zeros(Float64, nc_eta1)
    rho = zeros(Float64, nc_eta1)
    sval = zeros(Float64, nc_eta1)

    cos_k(x) = cos(mode * 2 * pi * x / Lx)
    sin_k(x) = sin(mode * 2 * pi * x / Lx)

    # Test Poisson
    #-------------
    # Set exact solution
    for i in 1:nc_eta1
        xi = eta1_min + (i - 1) * delta_eta1
        ex_exact[i] = sin_k(xi) / (2.0 * mode * pi / Lx)
    end

    compute_rhs_from_function!(rho, maxwell_1d, cos_k, deg)

    compute_e_from_rho!(ex, maxwell_1d, rho)

    # Evaluate spline curve at grid points and compute error
    # Ex is a 1-form, i.e. one spline degree lower
    sval = eval_uniform_periodic_spline_curve(deg - 1, ex)
    err_ex = maximum(abs.(sval .- ex_exact))
    println(" error Poisson  $err_ex ")
    @test err_ex ≈ 0.0 atol = 1e-6

    # Test Ampere
    #-------------
    # Set time step
    dt = 0.5 * delta_eta1
    # Set exact solution
    for i in 1:nc_eta1
        xi = eta1_min + (i - 1) * delta_eta1
        ex_exact[i] = -cos_k(xi) * dt
    end

    compute_rhs_from_function!(rho, maxwell_1d, cos_k, deg - 1)
    fill!(ex, 0.0)
    compute_e_from_j!(ex, maxwell_1d, dt .* rho, 1)

    # Evaluate spline curve at grid points and compute error
    # Ex is a 1-form, i.e. one spline degree lower
    sval = eval_uniform_periodic_spline_curve(deg - 1, ex)
    err_ex2 = maximum(abs.(sval .- ex_exact))
    #println( " error Ampere  $err_ex2 ")
    @test err_ex2 ≈ 0.0 atol = 1e-6

    l2norm = l2norm_squared(maxwell_1d, ex, deg - 1)
    err_l2norm = l2norm - dt^2 * pi
    #println( " error l2 norm $err_l2norm ")

    # Test Maxwell on By and Ez
    #--------------------------
    # Set time stepping parameters

    time = 0.0
    dt = 0.5 * delta_eta1
    nstep = 10

    # Compute initial fields
    ex = 0.0 # 0-form -> splines of degree deg
    l2projection!(bz, maxwell_1d, cos_k, deg - 1) # 0-form -> splines of degree deg-1

    for istep in 1:nstep
        compute_b_from_e!(bz, maxwell_1d, 0.5 * dt, ey)
        compute_e_from_b!(ey, maxwell_1d, dt, bz)
        compute_b_from_e!(bz, maxwell_1d, 0.5 * dt, ey)

        time = time + dt

        for i in 1:nc_eta1
            xi = eta1_min + (i - 1) * delta_eta1
            ey_exact[i] = sin(mode * 2π * xi / Lx) * sin(mode * 2π * time / Lx)
            bz_exact[i] = cos(mode * 2π * xi / Lx) * cos(mode * 2π * time / Lx)
        end

        sval = eval_uniform_periodic_spline_curve(deg, ey)
        err_ey = norm(sval .- ey_exact)

        sval = eval_uniform_periodic_spline_curve(deg - 1, bz)
        err_bz = norm(sval .- bz_exact)

        @test err_ey ≈ 0.0 atol = 1e-2
        @test err_bz ≈ 0.0 atol = 1e-2
    end
end


text:
T = [.80 .05 .15 # sunny
     .20 .60 .20 # rainy
     .20 .30 .50] # foggy

E = [.10 .90 # sunny
     .80 .20 # rainy
     .30 .70] # foggy

x = [2, 1, 1] # 1: umbrella 2: no umbrella
states = Dict(1 => "sunny", 2 => "rainy", 3 => "foggy")

function viterbi(x, states, T::Array{Float64, 2}, E::Array{Float64, 2})
    n = length(x)
    p = length(states)
    prior = 1 / p # assume equal priors

    delta = zeros(n, p) # holds the most likely path
    gamma = convert(Array{Int64, 2}, zeros(n, p)) # holds the arguments for the most likely path

    delta[1, :] = prior * E[:, 2]

    for i in 2:n
        b = x[i] # current observation

        for j in 1:p
            a = T[:, j] # transition probability *to* hidden state j {sunny, rainy, foggy}
            prev = delta[i-1, :] # previous deltas

            cont = diag(a * prev)
            delta[i, j] = maximum(cont) * E[j, b]
            _, gamma[i, j] = findmax(cont)
        end
    end

    _, dlast = findmax(delta[n, :])
    sequence = [dlast]

    for i in n:-1:2
        _, prev = findmax(delta[i, :])
        unshift!(sequence, gamma[i, prev])
    end

    return delta, [states[i] for i in sequence]
end


text:
"""
    coefficienttype(optimizer::AbstractSDOptimizer)

Returns the `coefficienttype` that should be used for `optimizer`.
"""
function coefficienttype end

"""
    getnumberofconstraints(optimizer::AbstractSDOptimizer)

Returns the number of constraints of the model.
"""
function getnumberofconstraints end

"""
    getnumberofblocks(optimizer::AbstractSDOptimizer)

Returns the number of blocks of the block matrix.
"""
function getnumberofblocks end

"""
    getblockdimension(optimizer::AbstractSDOptimizer, blk)

Returns the dimension of the block `blk`.
"""
function getblockdimension end

"""
    init!(optimizer::AbstractSDOptimizer, blkdims::Vector{Int}, nconstrs::Integer)

Initialize the optimizer with nconstrs constraints and blkdims blocks.
"""
function init! end

"""
    getconstraintconstant(optimizer::AbstractSDOptimizer, constr::Integer)

Sets the entry `constr` of `b` to `val`.
"""
function getconstraintconstant end

"""
    setconstraintconstant!(optimizer::AbstractSDOptimizer, val, constr::Integer)

Get the entry `constr` of `b`.
"""
function setconstraintconstant! end

"""
    getobjectivecoefficients(optimizer::AbstractSDOptimizer)

Return the list of entries `blk`, `i`, `j` of the objective matrix.
"""
function getobjectivecoefficients end

"""
    setobjectivecoefficient!(optimizer::AbstractSDOptimizer, val, blk::Integer, i::Integer, j::Integer)

Set the entry `i`, `j` of the block `blk` of the objective matrix to `val`.
"""
function setobjectivecoefficient! end

"""
    getconstraintcoefficients(optimizer::AbstractSDOptimizer, constr::Integer)

Return the list of entries `blk`, `i`, `j` of the matrix of the constraint `constr`.
"""
function getconstraintcoefficients end

"""
    setconstraintcoefficient!(optimizer::AbstractSDOptimizer, val, constr::Integer, blk::Integer, i::Integer, j::Integer)

Set the entry `i`, `j` of the block `blk` of the matrix of the constraint `constr` to `val`.
"""
function setconstraintcoefficient! end

"""
    getX(optimizer::AbstractSDOptimizer)

Returns the solution X as a block matrix.
"""
function getX end

"""
    gety(optimizer::AbstractSDOptimizer)

Returns the solution y.
"""
function gety end

"""
    getZ(optimizer::AbstractSDOptimizer)

Returns the solution Z.
"""
function getZ end

"""
    getprimalobjectivevalue(optimizer::AbstractSDOptimizer)

Returns the primal objective value.
"""
function getprimalobjectivevalue end

"""
    getdualobjectivevalue(optimizer::AbstractSDOptimizer)

Returns the dual objective value.
"""
function getdualobjectivevalue end


text:
using DukascopyTicksReader
using DukascopyTicksReader: get, to_arrays, to_dataframe, to_timearray
using DukascopyTicksReader: _cache_dir, _destination_filename, _url

using Base.Test

@testset "DukascopyTicksReader" begin

    @testset "utils" begin
        source = DukascopyTicks()
        cache = CacheDirectory()
        ticker = "EURUSD"

        dt = DateTime(2016, 3, 28, 0, 40)

        @test _cache_dir(source, cache, ticker, dt) == joinpath(homedir(), "data", "dukascopy", "ticks",
                               "2016", "2016-03", "2016-03-28", "2016-03-28_000000")

        @test _url(source, ticker, dt) == "http://www.dukascopy.com/datafeed/EURUSD/2016/03/28/00h_ticks.bi5"

        @test _destination_filename(source, ticker, dt) == "EURUSD.bi5"
    end

    @testset "usage: get" begin
        source = DukascopyTicks()
        cache = CacheDirectory()
        ticker = "USDCHF"

        dt = DateTime(2016, 3, 28, 0, 40)
        dt2 = DateTime(2016, 4, 8, 0, 40)

        data = get(source, ticker, dt)
        #println(data)

        #a = to_arrays(data)

        @testset "DataFrame" begin
            df = to_dataframe(data)
            # println(df)
        end

        @testset "TimeArray" begin
            ta = to_timearray(data)
            # println(ta)
        end

        #data = get(source, ticker, dt, dt2)
        #println(data)

        #=
        tickers = ["EURUSD", "USDJPY"]
        data = get(source, tickers, dt, dt2)
        println(data)
        =#
    end

    @testset "usage: bulk download" begin
        skip_error = true
        tickers = ["AUDUSD", "USDCAD"]
        # tickers = ["AUDUSD", "USDCAD", "USDCHF", "EURUSD", "GBPUSD", "NZDUSD"]
        source = DukascopyTicks()
        download(source, tickers, DateTime(2016, 1, 5), DateTime(2016, 1, 8), skip_error=skip_error)
        # download(source, tickers, DateTime(2016, 1, 1), DateTime(Dates.today()), skip_error=skip_error)
        # download(source, tickers, DateTime(2014, 1, 1), DateTime(Dates.today()), skip_error=skip_error)
    end

end


text:
function serial_random_increments!(counters, indices)
    for i in indices
        @inbounds counters[i] += 1
    end
end

random_increments!(counters, indices) =
    random_increments!(counters, indices, Val{false}(), lastindex(indices), -1)

function random_increments!(
    counters,
    indices,
    tsx::ValBool,
    nbatches::Integer,
    ntries::Integer,
)
    k = firstindex(indices)
    batchbound = lastindex(indices) - nbatches
    # nabt = 0
    while k <= batchbound
        if tsx === Val(true)
            for _ in 1:ntries
                code = xbegin()
                if code == XBEGIN_STARTED
                    for _ in 1:nbatches
                        @inbounds counters[indices[k]] += 1
                        k += 1
                    end
                    xend()
                    @goto continue_outer
                end
            end
        end
        # nabt += nbatches

        for _ in 1:nbatches
            i = @inbounds indices[k]
            UnsafeAtomics.add!(pointer(counters, i), one(eltype(counters)), monotonic)
            k += 1
        end

        @label continue_outer
    end

    while k <= lastindex(indices)
        if tsx === Val(true)
            for _ in 1:ntries
                code = xbegin()
                if code == XBEGIN_STARTED
                    @inbounds counters[indices[k]] += 1
                    k += 1
                    xend()
                    @goto continue_outer2
                end
            end
        end
        # nabt += 1

        i = @inbounds indices[k]
        UnsafeAtomics.add!(pointer(counters, i), one(eltype(counters)), monotonic)
        k += 1

        @label continue_outer2
    end

    # @show nabt / length(indices)
end

function random_increments_adaptive!(counters, indices, ntries = 4)
    k = firstindex(indices)
    batchsize = 1
    while true
        k <= lastindex(indices) || return

        kbatch = min(lastindex(indices), k + batchsize - 1)
        for _ in 1:ntries
            code = xbegin()
            if code == XBEGIN_STARTED
                while k <= kbatch
                    @inbounds counters[indices[k]] += 1
                    k += 1
                end
                xend()
                batchsize += 1  # stretch batch on success
                batchsize = max(1, min(batchsize, lastindex(counters) - k + 1))
                @goto continue_outer
            end
        end

        batchsize = max(1, batchsize ÷ 2)  # shrink batch on failure

        i = @inbounds indices[k]
        UnsafeAtomics.add!(pointer(counters, i), one(eltype(counters)), monotonic)
        k += 1

        @label continue_outer
    end
end

function parallel_random_increments!(counters, indices, ntasks, args...)
    chuncked_indices = Iterators.partition(indices, cld(length(indices), ntasks))
    @assert length(chuncked_indices) == ntasks
    foreach_thread(chuncked_indices) do chunk
        random_increments!(counters, chunk, args...)
    end
end

function parallel_random_increments_adaptive!(counters, indices, ntasks, args...)
    chuncked_indices = Iterators.partition(indices, cld(length(indices), ntasks))
    @assert length(chuncked_indices) == ntasks
    foreach_thread(chuncked_indices) do chunk
        random_increments_adaptive!(counters, chunk, args...)
    end
end

function random_increments_problem_uniform(ncounters, nindices)
    counters = zeros(Int, ncounters)
    indices = rand(1:ncounters, nindices)
    return (; counters, indices)
end

function random_increments_problem_clusters(ncounters, nindices, clustersize)
    counters = zeros(Int, ncounters)
    indices = Vector{Int}(undef, cld(nindices, clustersize) * clustersize)
    @assert length(indices) >= nindices
    rand!(@view(indices[1:clustersize:end]), 1:ncounters)
    for i in 2:clustersize
        @view(indices[i:clustersize:end]) .=
            mod1.(@view(indices[i-1:clustersize:end]) .+ 1, ncounters)
    end
    resize!(indices, nindices)
    @assert all(in(eachindex(counters)), indices)
    return (; counters, indices)
end


text:
module Ising

using Erdos
using Random, Statistics
using UnPack
using Printf

export run_bp

getref(v::Vector, i::Integer) = pointer(v, i)

include("ising_bp.jl")
include("ising_tap.jl")

end #module


text:
push!(LOAD_PATH, "/Users/czhang/Desktop/Projects/dw_/julialib/")
import DimmWitted
DimmWitted.set_libpath("/Users/czhang/Desktop/Projects/dw_/libdw_julia")

######################################
# The following piece of code creates a
# synthetic data set:
#    - Data type is Cdouble
#    - Modle type is Array{Cdouble}
#
nexp = 100000
nfeat = 1024
examples = Array(Cdouble, nexp, nfeat+1)
for row = 1:nexp
	for col = 1:nfeat
		examples[row, col] = 1
	end
	if rand() > 0.8
		examples[row, nfeat+1] = 0
	else
		examples[row, nfeat+1] = 1
	end
end
model = Cdouble[0 for i = 1:nfeat]

sparse_example=sparse(examples)

######################################
# Define the loss function and gradient
# function for logistic regression
#
function loss(row, model::Array{Cdouble,1})
	const nfeat = length(model)
	const lastcol = nfeat + 1
	const nnz = length(row)
	label = 0.0
	if row[nnz].idx == lastcol
		label = row[nnz].data
	end
	d = 0.0
	for i = 1:nnz
		if row[i].idx != lastcol
			d = d + row[i].data*model[row[i].idx]
		end
	end
	return (-label * d + log(exp(d) + 1.0))
end

function grad(row, model::Array{Cdouble,1})
	const nfeat = length(model)
	const lastcol = nfeat + 1
	const nnz = length(row)
	label = 0.0
	if row[nnz].idx == lastcol
		label = row[nnz].data
	end
	d = 0.0
	for i = 1:nnz
		if row[i].idx != lastcol
			d = d + row[i].data*model[row[i].idx]
		end
	end
	d = exp(-d)
	Z = 0.00001 * (-label + 1.0/(1.0+d))
	for i = 1:nnz
		if row[i].idx != lastcol
		model[row[i].idx] = model[row[i].idx] - row[i].data * Z
		end
	end
	return 1.0
end

######################################
# Create a DimmWitted object using data
# and model. You do not need to specify
# the type, they are infer'ed by the
# open() function, which is parametric.
#
dw = DimmWitted.open(sparse_example, model,
                DimmWitted.MR_PERMACHINE,
                DimmWitted.DR_SHARDING,
                DimmWitted.AC_ROW)

######################################
# Register functions.
#
handle_loss = DimmWitted.register_row(dw, loss)
handle_grad = DimmWitted.register_row(dw, grad)
#DimmWitted.register_model_avg(dw, handle_loss, avg, true)
#DimmWitted.register_model_avg(dw, handle_grad, avg, true)

######################################
# Run 10 epoches.
#
for iepoch = 1:10
	rs = DimmWitted.exec(dw, handle_loss)
	println("LOSS: ", rs/nexp)
	rs = DimmWitted.exec(dw, handle_grad)
end

max_stars_repo_name: GloriaGo/dw

max_issues_repo_name: GloriaGo/dw

max_forks_repo_name: GloriaGo/dw


text:
#because both A and B can again be of type UnionTrunc, this thing can also represent union(alg1,alg2,alg3)
"
    Take the union of 2 truncation algorithms using &
"
struct UnionTrunc{A,B} <: Algorithm
    alg1 :: A
    alg2 :: B
end

function changebonds(state,alg::UnionTrunc)
    state = changebonds(state,alg.alg1);
    state = changebonds(state,alg.alg2);
    return state
end

function changebonds(state,H,alg::UnionTrunc,pars = params(state,H))
    (state,pars) = changebonds(state,H,alg.alg1,pars)
    (state,pars) = changebonds(state,H,alg.alg2,pars)
    return (state,pars)
end

Base.:&(alg1::Algorithm,alg2::Algorithm) = UnionTrunc(alg1,alg2)


text:
using JuliaFormatter

function main()
    perfect = true
    # note: keep in sync with `.github/workflows/format-check.yml`
    for d in ["src/", "test/", "format/"]
        @info "...linting $d ..."
        dir_perfect = format(d; style=YASStyle())
        perfect = perfect && dir_perfect
    end
    if perfect
        @info "Linting complete - no files altered"
    else
        @info "Linting complete - files altered"
        run(`git status`)
    end
    return nothing
end

main()


text:
export iauRv2m
"""
Form the r-matrix corresponding to a given r-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   w        double[3]      rotation vector (Note 1)

Returned:
   r        double[3][3]    rotation matrix

Notes:

   1. A rotation matrix describes a rotation through some angle about
   some arbitrary axis called the Euler axis.  The "rotation vector"
   supplied to This function has the same direction as the Euler
   axis, and its magnitude is the angle in radians.

   2. If w is null, the unit matrix is returned.

   3. The reference frame rotates clockwise as seen looking along the
   rotation vector from the origin.

This revision:  2015 January 30

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRv2m(double w[3], double r[3][3])
function iauRv2m(w::AbstractVector{<:Real})
   r = zeros(Float64, 3, 3)

   ccall((:iauRv2m, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
          w, r)

   return SMatrix{3,3}(r')
end


text:
"""
	mutable struct Multiblock <: AbstractSpecies

A species representing a linear polymer chain with variable chain length and number of blocks.
"""
mutable struct Multiblock <: AbstractSpecies
	monomers      :: Vector{Monomer}
	mids          :: Vector{Int}
	mid_block     :: Vector{Int}	       # Block idx  --> monomer id
	mid_map       :: Dict{Int,Vector{Int}} # Monomer id --> block idx, blocks which contain a given mid

	# Chain structure
	Ns            :: Int
    N             :: Int
    N_block       :: Vector{Int}
    Nref          :: Float64
    Nref_block    :: Vector{Float64}
	f_block       :: Vector{Float64}
	b_block       :: Vector{Float64}
	ds_block      :: Vector{Float64}
	begin_block   :: Vector{Int}

	# SCF fields
	Q             :: Float64
	q             :: Vector{FieldGrid{Float64}}
	qc            :: Vector{FieldGrid{Float64}}
	density       :: Dict{Int,FieldGrid{Float64}} # Stores density per monomer type
	density_block :: Dict{Int,FieldGrid{Float64}} # Stores density per block idx
	system        :: Option{FieldSystem}
end

function Multiblock(mon_block::AbstractVector{Monomer}, N_block::AbstractVector{<:Integer},
    b_block::AbstractVector{<:Real}, Ns::Integer = 100)
    @assert length(mon_block) == length(N_block) == length(b_block)

    # Setup monomer info and block-type mappings
    monomers = sort(unique(mon_block))
	mids = [mon.id for mon in monomers]

    mid_block = [mon.id for mon in mon_block]
    mid_map = Dict()
    for (iblk, mid) in enumerate(mid_block)
        if !haskey(mid_map, mid)
            mid_map[mid] = []
        end
        push!(mid_map[mid], iblk)
    end

    # Determine contour locations for each block
    N = sum(N_block)
    Nref_block = [mon_block[iblk].vol * N_block[iblk] for iblk = 1:length(mon_block)]
    Nref, begin_block, ds_block = chain_grid(Nref_block, Ns) # Use the Nref block lengths to discretize solver
    f_block = [Nr/Nref for Nr in Nref_block]

    # Check if values were adjusted
    if Ns != begin_block[end]-1
        @warn "Non-integer chain discretization. Adjusting `ds` and `Ns` appropriately."
        Ns = begin_block[end]-1
    end

    return Multiblock(monomers, mids, mid_block, mid_map,
        Ns, N, N_block, Nref, Nref_block, f_block, b_block, ds_block, begin_block,
        0.0, [], [], Dict(), Dict(), nothing
    )
end

function Homopolymer(mon::Monomer, N::Integer, b::Real, Ns::Integer = 100)
    return Multiblock([mon], [N], [b], Ns)
end

function Diblock(amon::Monomer, bmon::Monomer, N::Integer, f::Real, b1::Real, b2::Real, Ns::Integer = 100)
    N_block = trunc.(Int, [f*N, (1-f)*N])
    if sum(N_block) != N
        error("Non-integer block lengths for N = $N, f = $f.")
    end
    return Multiblock([amon, bmon], N_block, [b1, b2], Ns)
end

#==============================================================================#

function Base.show(io::IO, chain::Multiblock)
	if nblocks(chain) == 1
        @printf(io, "Homopolymer(mid = %d, N = %d)", chain.mids[1], chain.N)
    elseif nblocks(chain) == 2
        @printf(io, "Diblock(mids = %s, N = %d, f = [%.3f, %.3f])", chain.mids, chain.N, chain.f_block[1], chain.f_block[2])
    else
        @printf(io, "Multiblock(%d blocks, %d monomers, N = %d)", nmonomers(chain), nblocks(chain), chain.N)
    end
end

nblocks(chain::Multiblock) = length(chain.N_block)

function monomer_fraction(chain::Multiblock, mid::Integer)
	# Check if mid present
	if !hasmonomer(chain, mid); return 0.0; end
	# It has the mid, so calculate fraction
	f = 0.0
	for iblk in chain.mid_map[mid]
		f += chain.f_block[iblk]
	end
	return f
end

function setup!(chain::Multiblock, sys::FieldSystem)
    chain.system = sys

	# Setup propagators
	chain.Q = 0.0
	chain.q = [zeros(Float64, sys.dims) for n = 1:chain.Ns+1]
	chain.qc = [zeros(Float64, sys.dims) for n = 1:chain.Ns+1]

	# Setup density grids
	for mid in chain.mids
		chain.density[mid] = zeros(Float64, sys.dims)
	end
	for iblk = 1:nblocks(chain)
		chain.density_block[iblk] = zeros(Float64, sys.dims)
	end

	return nothing
end

#==============================================================================#

function density!(chain::Multiblock)
	@assert !isnothing(chain.system)
	sys = chain.system
    plan = sys.fftplan

	# Necessary class fields
	q, qc = chain.q, chain.qc
	solver = sys.solver
	ksq = sys.cell.ksq

	# Solve q by integrating forward along chain
	q[1] .= 1.0
	for iblk = 1:nblocks(chain)
		b = chain.b_block[iblk]
		ds = chain.ds_block[iblk]

		mid = chain.mid_block[iblk]
		mon = sys.monomers[mid]
		omega = sys.fields[mid]

		update!(solver, omega, ksq, mon.vol, b, ds)
		for s = chain.begin_block[iblk] : chain.begin_block[iblk+1]-1
			propagate!(solver, plan, q[s], q[s+1])
		end
	end

	# Solve qc by integrating reverse along the chain
    # If a homopolymer, short-cut the integration and copy the propagator in reverse
    if nblocks(chain) == 1
        for (s, qi) in enumerate(Iterators.reverse(q))
            @inbounds qc[s] .= qi
        end
    else
    	qc[end] .= 1.0
    	for iblk = nblocks(chain):-1:1
    		b = chain.b_block[iblk]
    		ds = chain.ds_block[iblk]

			mid = chain.mid_block[iblk]
			mon = sys.monomers[mid]
			omega = sys.fields[mid]

    		update!(solver, omega, ksq, mon.vol, b, ds)
    		for s = chain.begin_block[iblk+1] : -1 : chain.begin_block[iblk]+1
    			propagate!(solver, plan, qc[s], qc[s-1])
    		end
    	end
    end

	# Partition function
	chain.Q = sum(chain.q[end]) / ngrid(sys)

	# Integrate density w/ Simpson's rule
	# Reset density grids
	for rho in values(chain.density); rho .= 0.0; end
	for rho in values(chain.density_block); rho .= 0.0; end

	for iblk = 1:nblocks(chain)
		# Block start and end indices
		blk_bgn = chain.begin_block[iblk]
		blk_end = chain.begin_block[iblk+1]

		ds = chain.ds_block[iblk]
		rho = chain.density_block[iblk]

		# First and last
		@simd for i in eachindex(rho)
			@inbounds rho[i] += q[blk_bgn][i] * qc[blk_bgn][i]
		end
		@simd for i in eachindex(rho)
			@inbounds rho[i] += q[blk_end][i] * qc[blk_end][i]
		end

		# Odd indices
		for s = blk_bgn+1:2:blk_end-1
            qs, qcs = q[s], qc[s]
			@simd for i in eachindex(rho)
				@inbounds rho[i] += 4.0 * qs[i] * qcs[i]
			end
		end

		# Even indices
		for s = blk_bgn+2:2:blk_end-2
            qs, qcs = q[s], qc[s]
			@simd for i in eachindex(rho)
				@inbounds rho[i] += 2.0 * qs[i] * qcs[i]
			end
		end

		# Normalize by ds/3 for Simpson's rule
		rho .*= (ds / 3.0)
	end

    # Normalize by 1/(Q*N)
    for rho in values(chain.density_block)
        rho ./= chain.Q * chain.Nref
    end

	# Add each block density to appropriate monomer type
	for mid in chain.mids
		for iblk in chain.mid_map[mid]
			chain.density[mid] .+= chain.density_block[iblk]
		end
	end

	return nothing
end

function scfstress(chain::Multiblock)
	@assert !isnothing(chain.system)
	sys = chain.system
	cell = sys.cell
	plan = sys.fftplan
	q, qc = chain.q, chain.qc

    # Get temp storage grids from the FFT helper
    qk1, qk2, qtmp = kgrid(plan,1), kgrid(plan,2), kgrid(plan,3)

	# Calculate variation in partition function by looping over all chain contour segments
	dQ = zeros(nparams(cell))
	for iblk = 1:nblocks(chain)
		# Block start and end indices
		blk_bgn = chain.begin_block[iblk]
		blk_end = chain.begin_block[iblk+1]

		b = chain.b_block[iblk]
		ds0 = chain.ds_block[iblk]

		# Block start and end indices
		for s = blk_bgn:blk_end
			# Adjust prefactor ds for Simpson's rule quadrature
			ds = ds0
			if s != blk_bgn && s!= blk_end
				if s % 2 == 0
					ds *= 4.0
				else
					ds *= 2.0
				end
			end

            # FFT the propagators to k-space grids
            qs = q[s]; qcs = qc[s]
            mul!(qk1, plan.FT, qs)
            mul!(qk2, plan.FT, qcs)

            # Apply dksq for each parameter
            for k = 1:nparams(cell)
                dksq = cell.dksq[k]
                @. qtmp = qk1 * dksq

                # Perform summation over kgrid propagators
                # Since we use half-space real FFT, the complex conjugate components do not cancel
                # We obtain the full result by doubling all the real components
                #   of the k != 0 waves
                qint = real(qk2[1]) * real(qtmp[1])
                @simd for i = 2:length(qk2)
                    @inbounds qint += 2 * real(qk2[i]) * real(qtmp[i])
                end

                dQ[k] -= b^2 * ds * qint / 6.0 / ngrid(sys)^2
            end

		end
	end

	# Normalize stress after completion
	dQ ./= 3.0 # Simpson's rule factor
	dQ ./= chain.Q * chain.Nref # Stress formula normalization

    return dQ
end


text:
# Server may be sleepy and still waking up; wait until we can ask it for `/registries`
t_start = time()
@info("Waiting for PkgServer liveness")

function handle_http_error(e)
    # New HTTP error code in 0.9.13+
    if isa(e, HTTP.TimeoutRequest.ReadTimeoutError)
        return -Base.Libc.ECONNRESET
    end

    if !isa(e, HTTP.IOExtras.IOError)
        rethrow(e)
    end

    # Peel it to get the inner exception
    e = e.e

    # If it just can't connect, deal with it silently, otherwise rethrow
    if isa(e, Base.IOError)
        if e.code in (-Base.Libc.ECONNREFUSED, -Base.Libc.EPIPE, -Base.Libc.ECONNRESET)
            return e.code
        end
    end
    # I don't know how this occurs, but it does sometimes.  Probably an HTTP bug.
    if isa(e, Base.EOFError)
        return -Base.Libc.ECONNRESET
    end

    # If it's none of the "whitelisted" errors, rethrow it.
    rethrow(e)
end

while true
    # Try to get an HTTP 200 OK on /registries
    response_code = try
        HTTP.get("$(server_url)/registries"; retry = false, readtimeout=1).status
    catch e
        handle_http_error(e)
    end

    if response_code == 200
        break
    end

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

    # Sleep a bit between attempts
    sleep(0.2)
end

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

    # Test asking for that registry directly, unpacking it and verifying the treehash
    mktemp() do tarball_path, tarball_io
        response = HTTP.get("$(server_url)/registry/$(registry_uuid)/$(registry_treehash)"; response_stream=tarball_io)
        close(tarball_io)
        @test response.status == 200
        @test registry_treehash == Tar.tree_hash(open(pipeline(`cat $(tarball_path)`, `gzip -d`), read=true))
    end

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

    # Next, hit the `/meta` endpoint, ensure that the version it reports matches with what we expect:
    response = HTTP.get("$(server_url)/meta")
    @test response.status == 200
    meta = JSON3.read(String(response.body))
    @test haskey(meta, "pkgserver_version")
    @test meta["pkgserver_version"] == PkgServer.get_pkgserver_version()
    @test haskey(meta, "pkgserver_url")
    @test meta["pkgserver_url"] == "https://starfleet-central.pkg.julialang.org"
    @test meta["registry_update_task"] == "started"

    # Also hit the `/meta/siblings` endpoint
    response = HTTP.get("$(server_url)/meta/siblings")
    @test response.status == 200
    siblings = collect(JSON3.read(String(response.body)))
    @test "https://us-west.pkg.julialang.org" in siblings
    @test "https://in.pkg.julialang.org" in siblings
    @test "https://au.pkg.julialang.org" in siblings

    response = HTTP.get("$(server_url)/robots.txt")
    @test response.status == 200
    @test chomp(String(response.body)) == "User-agent: * Disallow: /"

    # Ensure that some random URL gets a 404
    @test_throws HTTP.ExceptionRequest.StatusError HTTP.get("$(server_url)/docs")

    # Test a dynamically-generated artifact TOML
    art_tree_hash = "4ed4e6caa3c4559f34d9eafd2f42e9863f83b573"
    art_name = "Xorg_xineramaproto"
    response = HTTP.get("$(server_url)/artifact/$(art_tree_hash)/$(art_name)")
    @test response.status == 200
    art_toml = TOML.parse(String(response.body))
    @test haskey(art_toml, art_name)
    @test haskey(art_toml[art_name], "git-tree-sha1")
    @test art_toml[art_name]["git-tree-sha1"] == art_tree_hash
    @test haskey(art_toml[art_name], "download")
end

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

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

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

    # Sleep until there's nothing left in the `temp` directory.
    t_start = time()
    function any_inprogress_files(temp_dir)
        for (root, dirs, files) in walkdir(temp_dir)
            for f in files
                if endswith(f, ".inprogress")
                    return true
                end
            end
        end
        return false
    end
    while any_inprogress_files(joinpath(cache_dir, "..", "temp"))
        sleep(0.1)
        if time() - t_start >= 10
            error("Timed out waiting for temp directory to finish")
        end
    end

    # Test that the resources we expect to exist are in fact cached on the server.
    # We just test that a few UUIDs and tree hashes exist.  Note we use the exact packages
    # listed above, otherwise we may not know what the treehash actually should be.
    pkg_uuid_treehashes = [
        # Gtk
        "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" => "ab2676a65345f8813be07675aa657464f98e6704",
        # DotNET
        "5091b313-875f-492b-8fe1-0f0d7d725ad8" => "58259eee27b838bfbf234cbfd043e3f090389f66",
        # FIGlet
        "3064a664-84fe-4d92-92c7-ed492f3d8fae" => "bfc6b52f75b4720581e3e49ae786da6764e65b6a"
    ]
    for (uuid, treehash) in pkg_uuid_treehashes
        @test isfile(joinpath(cache_dir, "package", uuid, treehash))
    end

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

    # Ensure that, when we hit `/meta` now, the server knows that it has a bunch of packages and artifacts:
    response = HTTP.get("$(server_url)/meta")
    @test response.status == 200
    meta = JSON3.read(String(response.body))
    @test haskey(meta, "julia_version")
    @test VersionNumber(meta["julia_version"]) >= v"1.3"
    @test haskey(meta, "pkgserver_version")
    @test meta["pkgserver_version"] == PkgServer.get_pkgserver_version()
end

@testset "Access Tracking" begin
    # Test that the `/meta/stats` endpoint works as expected
    response = HTTP.get("$(server_url)/meta/stats")
    @test response.status == 200
    stats = JSON3.read(String(response.body))
    @test haskey(stats, "packages_cached")
    @test stats["packages_cached"] >= 70
    @test haskey(stats, "artifacts_cached")
    @test stats["artifacts_cached"] >= 30

    # Find figlet fonts endpoint, hit it a couple of times, ensure that the LRU count goes up each time
    figlet_fonts_resource = "artifact/125ac0315d68bbb612f8c2189ea83401f73238f0"
    figlet_fonts_entry = stats["lru"][figlet_fonts_resource]
    @test figlet_fonts_entry["num_accessed"] > 0

    # Now, hit this a couple more times.  We'll hit it a few more times now, while it's still streaming,
    # to trigger the streaming codepaths
    for idx in 1:5
        HTTP.get("$(server_url)/$(figlet_fonts_resource)")
    end

    # Refresh the stats
    response = HTTP.get("$(server_url)/meta/stats")
    @test response.status == 200
    stats = JSON3.read(String(response.body))

    # Ensure that the counter went up, and that the timing has advanced: (these are no longer true with nginx in the way)
    #@test stats["lru"][figlet_fonts_resource]["num_accessed"] >= figlet_fonts_entry["num_accessed"] + 5
    #@test stats["lru"][figlet_fonts_resource]["last_accessed"] > figlet_fonts_entry["last_accessed"]

    # Test that payload bytes received and transmitted are both nonzero,
    # and that (since we've served the same resources multiple times) the
    # latter is greater than the former
    @test haskey(stats, "payload_bytes_received")
    @test stats["payload_bytes_received"] > 0
    @test haskey(stats, "payload_bytes_transmitted")
    @test stats["payload_bytes_transmitted"] > 0
    # This is not necessarily true with nginx in the way
    #@test stats["payload_bytes_transmitted"] > stats["payload_bytes_received"]
end

@testset "Skip-nonskip ambiguity testing" begin
    art_yskip_hash = "00499ee910a92cd27ecab7620029a802136c1048"
    art_nskip_hash = "2da0ddeae4275db146c85efe7310b1d3148938d1"

    # Hit a resource that contains an artifact that we know has empty directories within it.
    # Ensure that it exists both at its nskip and its yskip hashes, because it was originally
    # published under the yskip hash.
    mktemp() do tarball_path, tarball_io
        response = HTTP.get("$(server_url)/artifact/$(art_yskip_hash)"; response_stream=tarball_io)
        close(tarball_io)
        @test response.status == 200
        @test art_nskip_hash == Tar.tree_hash(open(pipeline(`cat $(tarball_path)`, `gzip -d`), read=true))

        # Also test that it's available at its nskip hash:
        @test HTTP.head("$(server_url)/artifact/$(art_yskip_hash)").status == 200
    end
end

@testset "Partial Content" begin
    # Example@0.5.3
    uuid = "7876af07-990d-54b4-ab0e-23690620f79a"
    treehash = "46e44e869b4d90b96bd8ed1fdcf32244fddfb6cc"
    content_url = "$(server_url)/package/$uuid/$treehash"
    cache_path = joinpath(cache_dir, "package", uuid, treehash)
    # Get full file
    full_resp = HTTP.get(content_url)
    @test full_resp.status == 200
    full_length = length(full_resp.body)
    @test full_length > 2048
    # Specifying both startbyte and stopbyte
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=0-1023"]);
    @test partial_resp.status == 206
    @test HTTP.header(partial_resp, "Content-Range") == "bytes 0-1023/$(full_length)"
    @test HTTP.header(partial_resp, "Content-Length") == "1024"
    @test partial_resp.body == full_resp.body[1:1024]
    # Specifying startbyte only
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=1024-"]);
    @test partial_resp.status == 206
    @test HTTP.header(partial_resp, "Content-Range") == "bytes 1024-$(full_length-1)/$(full_length)"
    @test HTTP.header(partial_resp, "Content-Length") == string(full_length - 1024)
    @test partial_resp.body == full_resp.body[1025:end]
    # Specifying suffix length only
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=-2047"]);
    @test partial_resp.status == 206
    @test HTTP.header(partial_resp, "Content-Range") == "bytes 638-2684/$(full_length)"
    @test HTTP.header(partial_resp, "Content-Length") == "2047"
    @test partial_resp.body == full_resp.body[end-2046:end]
    # Stopbyte larger than filesize
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=0-10000"]);
    @test partial_resp.status == 206
    @test HTTP.header(partial_resp, "Content-Range") == "bytes 0-$(full_length-1)/$(full_length)"
    @test HTTP.header(partial_resp, "Content-Length") == string(full_length)
    @test partial_resp.body == full_resp.body
    # Edgecase: startbyte larger than stopbyte
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=2-1"]; status_exception = false);
    @test partial_resp.status == 416

    # Specifying both startbyte and stopbyte but streaming rather than cached
    rm(cache_path; force=true)
    partial_resp = HTTP.get(content_url, ["Range"=>"bytes=0-1023"]);
    @test partial_resp.status == 206
    @test HTTP.header(partial_resp, "Content-Range") == "bytes 0-1023/$(full_length)"
    @test HTTP.header(partial_resp, "Content-Length") == "1024"
    @test partial_resp.body == full_resp.body[1:1024]
end


text:
# Mostly copied from https://github.com/IainNZ/Humanize.jl/blob/master/src/Humanize.jl#L27
function _datasize(bytes::Number)
    base = 1024.0
    nbytes = float(bytes)
    unit = base
    suffix = first(DATA_SUFFIX)

    for (i, s) in enumerate(DATA_SUFFIX)
        unit = base ^ i

        if nbytes < unit
            suffix = s
            break
        end
    end

    return @sprintf("%.1f%s", (base * nbytes / unit), suffix)
end