text:
# Compressed Matrix
# =================

"""
A data type to represent a count matrix with lossless compression.
"""
struct CMatrix{T<:Integer} <: AbstractMatrix{T}
    data::Vector{UInt8}
    size::Tuple{Int,Int}
end

CMatrix(counts::Matrix; kwargs...) = CMatrix{eltype(counts)}(blosc_compress(counts; kwargs...), size(counts))

Base.convert(::Type{CMatrix{T}}, M::AbstractMatrix) where {T} = CMatrix(Matrix{T}(M))

Base.size(M::CMatrix) = M.size
Base.IndexStyle(::Type{<:CMatrix}) = IndexLinear()
Base.getindex(M::CMatrix, i::Integer) = getitem(M.data, eltype(M), i)
Base.getindex(M::CMatrix, ::Colon, j::Integer) = getitems(M.data, eltype(M), (j - 1) * size(M, 1) + 1, size(M, 1))
function Base.getindex(M::CMatrix, ::Colon, cols::AbstractVector{<:Integer})
    n = length(cols)
    A = Matrix{eltype(M)}(undef, size(M, 1), n)
    for j in 1:n
        A[:,j] .= M[:,cols[j]]
    end
    return A
end

const BLOSC_MAX_OVERHEAD = 16  # BLOSC_MIN_HEADER_LENGTH

function blosc_compress(
        data::Matrix{T};
        level::Integer=9,
        shuffle::Bool=true,
        compressor::String="lz4hc",
        blocksize::Integer=size(data, 1) * sizeof(T),
        nthreads::Integer=1,
        mmap::Bool=false,
    ) where {T<:Integer}
    srcsize = sizeof(data)
    src = pointer(data)
    dstsize = srcsize + BLOSC_MAX_OVERHEAD
    dst = Vector{UInt8}(undef, dstsize)
    ret = ccall(
        (:blosc_compress_ctx, libblosc),
        Cint,
        (Cint, Cint, Csize_t,
         Csize_t, Ptr{Cvoid}, Ptr{Cvoid},
         Csize_t, Cstring,
         Csize_t, Cint),
        level, shuffle, sizeof(T),
        srcsize, src, dst,
        dstsize, compressor,
        blocksize, nthreads)
    @assert ret ≥ 0
    if mmap
        dst_mmap = Mmap.mmap(Vector{UInt8}, (ret,), shared=false)
        copyto!(dst_mmap, 1, dst, 1, ret)
        return dst_mmap
    else
        resize!(dst, ret)
        return dst
    end
end

function getitem(data::Vector{UInt8}, ::Type{T}, i::Integer) where T
    src = pointer(data)
    dst = Ref{T}()
    ret = ccall(
        (:blosc_getitem, libblosc),
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ref{T}),
        src, i - 1, 1, dst)
    @assert ret == sizeof(T)
    return dst[]
end

function getitems(data::Vector{UInt8}, ::Type{T}, i::Integer, n::Integer) where T
    src = pointer(data)
    dst = Vector{T}(undef, n)
    ret = ccall(
        (:blosc_getitem, libblosc),
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ref{T}),
        src, i - 1, n, dst)
    @assert ret == sizeof(T) * n
    return dst
end


text:
module Functions

using ..TopOpt: dim, whichdevice, CPU, GPU, TopOpt, PENALTY_BEFORE_INTERPOLATION
using ..TopOptProblems, ..FEA, ..CheqFilters
using ..Utilities, ForwardDiff, LinearAlgebra, Requires
using Parameters: @unpack
using TimerOutputs, Ferrite, StaticArrays
using StatsFuns, MappedArrays, LazyArrays
using ..TopOptProblems: getdh
using SparseArrays, Statistics, ChainRulesCore, Zygote
using Nonconvex: Nonconvex

export  Volume,
        Compliance,
        MeanCompliance,
        BlockCompliance,
        AbstractFunction,
        getfevals,
        getmaxfevals,
        maxedfevals,
        MicroVonMisesStress,
        MacroVonMisesStress,
        project,
        generate_scenarios,
        hutch_rand!,
        hadamard!

const to = TimerOutput()

abstract type AbstractFunction{T} <: Nonconvex.AbstractFunction end

include("compliance.jl")
include("volume.jl")
include("stress.jl")
include("trace.jl")
include("mean_compliance.jl")
include("block_compliance.jl")

end

ext: jl


text: abstract type Storage <: StaticInjection end


text:
using Pkg

dependencies = [
    "FFTW",
    "Parameters",
    "DifferentialEquations",
    "Plots",
    "Statistics",
    "PyCall"
]

Pkg.add(dependencies)


text:
using SparseArrays

function build_metric_induced_graph(metric, graph)
    n = nv(graph)

    metric_induced_graph = SimpleGraph(n)

    for u in 1:n
        for v in u:n
            if metric(u, v) <= 1 || has_edge(graph, u, v) || has_edge(graph, v, u)
                add_edge!(metric_induced_graph, u, v)
            end
        end
    end

    return metric_induced_graph
end

function build_simple_graph(graph)
    SimpleGraph(adjacency_matrix(graph) + transpose(adjacency_matrix(graph)))
end

function get_simple_digraph(graph)
    SimpleDiGraph(adjacency_matrix(graph))
end

function get_simple_weighted_graph(graph)
    SimpleWeightedGraph(convert(SparseMatrixCSC{Real, Int64}, max.(adjacency_matrix(graph), transpose(adjacency_matrix(graph)))))
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function molecule(::Molecule"ammonia")
  return """
N	0.0000000	0.0000000	0.1173470
H	0.0000000	0.9326490	-0.2738090
H	0.8076980	-0.4663250	-0.2738090
H	-0.8076980	-0.4663250	-0.2738090
  """
end


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

mutable struct IoK8sApiAdmissionregistrationV1beta1Webhook <: SwaggerModel
    admissionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: admissionReviewVersions
    clientConfig::Any # spec type: Union{ Nothing, IoK8sApiAdmissionregistrationV1beta1WebhookClientConfig } # spec name: clientConfig
    failurePolicy::Any # spec type: Union{ Nothing, String } # spec name: failurePolicy
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAdmissionregistrationV1beta1RuleWithOperations} } # spec name: rules
    sideEffects::Any # spec type: Union{ Nothing, String } # spec name: sideEffects
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiAdmissionregistrationV1beta1Webhook(;admissionReviewVersions=nothing, clientConfig=nothing, failurePolicy=nothing, name=nothing, namespaceSelector=nothing, rules=nothing, sideEffects=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("admissionReviewVersions"), admissionReviewVersions)
        setfield!(o, Symbol("admissionReviewVersions"), admissionReviewVersions)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("clientConfig"), clientConfig)
        setfield!(o, Symbol("clientConfig"), clientConfig)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("failurePolicy"), failurePolicy)
        setfield!(o, Symbol("failurePolicy"), failurePolicy)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("sideEffects"), sideEffects)
        setfield!(o, Symbol("sideEffects"), sideEffects)
        validate_property(IoK8sApiAdmissionregistrationV1beta1Webhook, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiAdmissionregistrationV1beta1Webhook

const _property_map_IoK8sApiAdmissionregistrationV1beta1Webhook = Dict{Symbol,Symbol}(Symbol("admissionReviewVersions")=>Symbol("admissionReviewVersions"), Symbol("clientConfig")=>Symbol("clientConfig"), Symbol("failurePolicy")=>Symbol("failurePolicy"), Symbol("name")=>Symbol("name"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("rules")=>Symbol("rules"), Symbol("sideEffects")=>Symbol("sideEffects"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiAdmissionregistrationV1beta1Webhook = Dict{Symbol,String}(Symbol("admissionReviewVersions")=>"Vector{String}", Symbol("clientConfig")=>"IoK8sApiAdmissionregistrationV1beta1WebhookClientConfig", Symbol("failurePolicy")=>"String", Symbol("name")=>"String", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("rules")=>"Vector{IoK8sApiAdmissionregistrationV1beta1RuleWithOperations}", Symbol("sideEffects")=>"String", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1beta1Webhook }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1beta1Webhook))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1beta1Webhook }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1beta1Webhook[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1beta1Webhook }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1beta1Webhook[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1beta1Webhook)
    (getproperty(o, Symbol("clientConfig")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

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


text:
using Gadfly

# wide-form plotting of matrices

set_default_plot_size(6inch, 4inch)

n = 20
m = 40
X = randn(m, n) * diagm(1:n)

plot(X, x=Row.index, y=Col.value, color=Col.index, Geom.line)


text:
using ReachabilityAnalysis, SparseArrays

using ReachabilityAnalysis: TaylorModelReachSet, AbstractLazyReachSet

const μ = 3.986e14 * 60^2
const r = 42164.0e3
const r² = r^2
const mc = 500.0
const n² = μ / r^3
const n = sqrt(n²)

const two_n = 2*n
const μ_r² = μ/r²

# columns correspond to x, y, vx, vy, t
const K₁ = [-28.8287 0.1005 -1449.9754 0.0046 0.0;
            -0.087 -33.2562 0.00462 -1451.5013 0.0]
const K₂ = [-288.0288 0.1312 -9614.9898 0.0 0.0;
            -0.1312 -288.0 0.0 -9614.9883   0.0]

const K₁mc = K₁/mc
const K₂mc = K₂/mc

function mymul!(v, A, x)
    @inbounds for ind = 1:length(v)
        v[ind] = zero(x[1])
        for jind = 1:length(x)
            v[ind] += A[ind, jind] * x[jind]
        end
    end
    return nothing
end

# dynamics in the 'approaching' mode
@taylorize function spacecraft_approaching!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₁mc, u)

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    # t' = 1
    du[5] = one(x)

    return du
end

# dynamics in the 'rendezvous attempt' mode
@taylorize function spacecraft_attempt!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    uxy = Vector{typeof(x)}(undef, 2)
    mymul!(uxy, K₂mc, u)

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    du[3] = (n²*x + two_n*vy) + ((μ_r² - μ_rc³*rx) + uxy[1])

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    du[4] = (n²*y - two_n*vx) - (μ_rc³*y - uxy[2])

    # t' = 1
    du[5] = one(x)

    return du
end

# dynamics in the 'aborting' mode
@taylorize function spacecraft_aborting!(du, u, p, t)
    x, y, vx, vy, t = u

    rx = r + x
    rx² = rx^2
    y² = y^2
    rc = sqrt(rx² + y²)
    rc³ = rc^3
    μ_rc³ = μ / rc³

    # x' = vx
    du[1] = vx

    # y' = vy
    du[2] = vy

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x)
    du[3] = (n²*x + two_n*vy) + (μ_r² - μ_rc³*rx)

    # vy' = n²y - 2n*vx - μ/(rc^3)y
    du[4] = (n²*y - two_n*vx) - μ_rc³*y

    # t' = 1
    du[5] = one(x)

    return du
end

function spacecraft(; X0 = Hyperrectangle([-900., -400., 2.5, 2.5, 0.],
                                          [25., 25., 2.5, 2.5, 0.]),
                         init=[(1, X0)],
                         abort_time=(120.0, 150.0))

    # variables
    x = 1   # x position
    y = 2   # y position
    vx = 3  # x velocity
    vy = 4  # y velocity
    t = 5   # time
    n = 4 + 1  # number of variables
    t_abort_lower, t_abort_upper = abort_time[1], abort_time[2]

    automaton = LightAutomaton(3)

    # mode 1 "approaching"
    invariant = HalfSpace(sparsevec([x], [1.], n), -100.) # x <= -100
    approaching = @system(x' = spacecraft_approaching!(x), dim:5, x ∈ invariant)

    # mode 2 ("rendezvous attempt")
    invariant = HalfSpace(sparsevec([x], [-1.], n), 100.) # x >= -100
    attempt = @system(x' = spacecraft_attempt!(x), dim:5, x ∈ invariant)

    # mode 3 "aborting"
    invariant = Universe(n)
    aborting = @system(x' = spacecraft_aborting!(x), dim: 5, x ∈ invariant)

    # transition "approach" -> "attempt"
    add_transition!(automaton, 1, 2, 1)
    guard = HalfSpace(sparsevec([x], [-1.], n), 100.) # x >= -100
    t1 = @map(x -> x, dim: n, x ∈ guard)

    # transition "approach" -> "abort"
    add_transition!(automaton, 1, 3, 2)
    guard_time = HPolyhedron([HalfSpace(sparsevec([t], [-1.], n), -t_abort_lower),  # t >= t_abort_lower
                              HalfSpace(sparsevec([t], [1.], n), t_abort_upper)])   # t <= t_abort_upper
    t2 = @map(x -> x, dim: n, x ∈ guard_time)

    # transition "attempt" -> "abort"
    add_transition!(automaton, 2, 3, 3)
    t3 = @map(x -> x, dim: n, x ∈ guard_time)

    H = HybridSystem(automaton=automaton,
                     modes=[approaching, attempt, aborting],
                     resetmaps=[t1, t2, t3])

    return InitialValueProblem(H, init)
end

# variables
const x, y, vx, vy, t = 1:5
const numvars = 5
const tan30 = tand(30)

LineOfSightCone = HPolyhedron([HalfSpace(sparsevec([x], [-1.], numvars), 100.),   # x >= -100
                    HalfSpace(sparsevec([x, y], [tan30, -1.], numvars), 0.),  # y >= x tan(30°)
                    HalfSpace(sparsevec([x, y], [tan30, 1.], numvars), 0.),   # -y >= x tan(30°)
                   ])

Target = Hyperrectangle(zeros(2), [0.2, 0.2]);

function line_of_sight(sol)
    all_idx = findall(x -> x == 2, location.(sol))  # attempt
    for idx in all_idx
        verif = all(set(R) ⊆ LineOfSightCone for R in sol[idx][2:end])
        !verif && return false
    end
    return true
end

function absolute_velocity(R::AbstractLazyReachSet)
    vx = 3
    vy = 4
    vx2 = set(overapproximate(project(R, vars=(vx,)), Interval)).dat
    vy2 = set(overapproximate(project(R, vars=(vy,)), Interval)).dat
    sqrt(vx2 + vy2)
end

function absolute_velocity(R::TaylorModelReachSet)
    Z = overapprximate(R, Zonotope)
    absolute_velocity(Z)
end

function velocity_constraint(sol)
    all_idx = findall(x -> x == 2, location.(sol)) # attempt
    for idx in all_idx
        # maximum velocity measured in m/min
        verif = all(absolute_velocity(R) < 0.055 * 60. for R in sol[idx])
        !verif && return false
    end
    return true
end

function target_avoidance(sol)
    all_idx = findall(x -> x == 3, location.(sol)) # aborting
    for idx in all_idx
        verif = all(is_intersection_empty(set(Projection(R, [x, y])), Target) for R in sol[idx])
        !verif && return false
    end
    return true
end


text:
using SeisRequests, Test
using Seis
using Dates: DateTime, Second, Millisecond
import HTTP
import StationXML

"""
Compare two Seis.Stations and return `true` if they are equal
apart from the `created` property of any StationXML in their `meta` fields.
"""
function stations_equal(s1, s2)
    # First test the StationXML apart from `created`
    if all(x->hasproperty(x.meta, :stationxml), (s1, s2))
        if !all(getproperty(s1.meta.stationxml, p) ==
                    getproperty(s2.meta.stationxml, p)
                for p in propertynames(s1.meta.stationxml) if p !== :created)
            return false
        end
    elseif any(x->hasproperty(x, :stationxml), (s1, s2))
        return false
    end
    # Remove stationxml field if any
    s1, s2 = deepcopy(s1), deepcopy(s2)
    s1.meta.stationxml = s2.meta.stationxml
    # Remove the request field if any
    # s1.meta.request = s2.meta.request = missing
    # Test remaining parts
    s1 == s2
end

@testset "get_stations" begin
    server = "http://service.example.com"
    SUCCESS = SeisRequests.CODE_SUCCESS

    @testset "parse_station_response" begin
        @testset "No response" begin
            request = FDSNStation(nodata=204)
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_station_response(Float64, request,
                response, server) == []

            @testset "Element type $T" for T in (Float32, Float64)
                @test SeisRequests.parse_station_response(T, request, response,
                    server) isa Vector{Seis.GeogStation{T}}
            end
            response.body = [1]
            @test (@test_logs (:warn,
                            "unexpected data in response reporting no data") SeisRequests.parse_station_response(
                            Float64, request, response, server) == [])
            request = FDSNStation(nodata=404)
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_station_response(Float64, request, response,
                server) == []
        end

        @testset "'Success' but no data" begin
            request = FDSNStation()
            response = HTTP.Messages.Response(SUCCESS)
            @test_throws ErrorException SeisRequests.parse_station_response(Float64,
                request, response, server)
        end

        @testset "Content-Type" begin
            @testset "Text" begin
                request = FDSNStation(network="AB", level="network", format="text")
                body = "AB|XXX|2000-01-01|3000-01-01|2\n"

                @testset "Wrong" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"text/plain"), body=body)
                    response_type_wrong = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"WRONG"), body=body)
                    @test (@test_logs (:warn, "content type of response is \"WRONG\", " *
                             "not \"text/plain\" as expected"
                        ) SeisRequests.parse_station_response(
                            Float64, request, response_type_wrong, server)) ==
                        SeisRequests.parse_station_response(Float64, request,
                            response_correct, server)
                end

                @testset "Empty" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"text/plain"), body=body)
                    response_type_empty = HTTP.Messages.Response(SUCCESS,
                        Dict(), body=body)
                    @test SeisRequests.parse_station_response(Float64,
                        request, response_type_empty, server) ==
                        SeisRequests.parse_station_response(Float64,
                            request, response_correct, server)
                end
            end
            @testset "XML" begin
                request = FDSNStation(network="AB", level="network")
                body = """<?xml version="1.0" encoding="UTF-8"?>
                    <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0">
                    <Source>Me</Source>
                    <Created>2000-01-01T00:00:00</Created>
                    <Network code="AN">
                    </Network>
                    </FDSNStationXML>"""
                @testset "Wrong" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"application/xml"), body=body)
                    response_type_wrong = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"WRONG"), body=body)
                    @test_logs (:warn, "content type of response is \"WRONG\", " *
                            "not \"application/xml\" as expected"
                        ) SeisRequests.parse_station_response(
                            Float64, request, response_type_wrong, server) ==
                        SeisRequests.parse_station_response(Float64, request,
                            response_correct, server)
                end
                @testset "Empty" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"application/xml"), body=body)
                    response_type_empty = HTTP.Messages.Response(SUCCESS, body)
                    @test SeisRequests.parse_station_response(
                            Float64, request, response_type_empty, server) ==
                        SeisRequests.parse_station_response(Float64, request,
                            response_correct, server)
                end
            end
        end

        @testset "Text format" begin
            @testset "Network level" begin
                request = FDSNStation(network="?N", level="network",
                    format="text")
                body = """
                    # Network|Description|StartTime|EndTime|TotalStations
                    # Arbitrary comment
                    AN|Network 1|2000-01-01T00:00:00.123|2002-02-03T04:05:06.789|99
                    # Another arbitrary comment
                    BN|Network 2|3000-01-02T03:04:05.678||101
                    """
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"text/plain"), body=body)
                    out = SeisRequests.parse_station_response(T, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{T}}
                    @test length(out) == 2
                    @test out[1] == Seis.GeogStation{T}(net="AN",
                        meta=Dict{Symbol,Any}(:network_description=>"Network 1",
                            :startdate=>DateTime(2000, 1, 1, 0, 0, 0, 123),
                            :enddate=>DateTime(2002, 2, 3, 4, 5, 6, 789),
                            :network_number_stations=>99,
                            :server=>server))
                    @test out[2] == Seis.GeogStation{T}(net="BN",
                        meta=Dict{Symbol,Any}(:network_description=>"Network 2",
                            :startdate=>DateTime(3000, 1, 2, 3, 4, 5, 678),
                            :network_number_stations=>101,
                            :server=>server))
                end
            end
            @testset "Station level" begin
                request = FDSNStation(network="IU", station="ANMO", level="station",
                    format="text")
                body = """
                    IU|ANMO|34.9459|-106.4572|1850.0|Albuquerque, New Mexico, USA|1989-08-29T00:00:00|1995-07-14T00:00:00
                    IU|ANMO|34.9459|-106.4572|1850.0|Albuquerque, New Mexico, USA|1995-07-14T00:00:00|
                    """
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Messages.Response(SUCCESS, body)
                    out = SeisRequests.parse_station_response(T, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{T}}
                    @test length(out) == 2
                    @test out[1] == Seis.GeogStation{T}(net="IU", sta="ANMO",
                        lat=34.9459, lon=-106.4572, elev=1850, meta=Dict(
                            :sitename=>"Albuquerque, New Mexico, USA",
                            :startdate=>DateTime(1989, 8, 29),
                            :enddate=>DateTime(1995, 7, 14), :server=>server))
                    @test out[2] == Seis.GeogStation{T}(net="IU", sta="ANMO",
                        lat=34.9459, lon=-106.4572, elev=1850, meta=Dict(
                            :sitename=>"Albuquerque, New Mexico, USA",
                            :startdate=>DateTime(1995, 7, 14), :server=>server))
                end
            end
            @testset "Channel level" begin
                request = FDSNStation(network="?N", station="FA*", channel="*",
                    level="channel", format="text")
                body = """
                    IU|COLA|20|HNE|64.873599|-147.8616|200.0|0.0|90.0|0.0|Kinemetrics FBA-23 Low-GainSensor|53687.1|1.0|M/S**2|80.0|2005-09-28T22:00:00|2009-07-08T22:00:00
                    IU|COLA|20|HNN|64.873599|-147.8616|200.0|0.0|0.0|0.0|Kinemetrics FBA-23 Low-GainSensor|53687.1|1.0|M/S**2|80.0|2005-09-28T22:00:00|
                    """
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Messages.Response(SUCCESS, body)
                    out = SeisRequests.parse_station_response(Float64, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{Float64}}
                    @test length(out) == 2
                    out_true = [Seis.GeogStation{Float64}(net="IU", sta="COLA",
                            loc="20", cha="HNE", lat=64.873599, lon=-147.8616, elev=200,
                            dep=0, azi=90, inc=90, meta=Dict(
                                :sensor_description=>"Kinemetrics FBA-23 Low-GainSensor",
                                :scale=>53687.1, :scale_frequency=>1.0,
                                :scale_units=>"M/S**2", :sample_rate=>80.0,
                                :startdate=>DateTime(2005, 9, 28, 22),
                                :enddate=>DateTime(2009, 7, 8, 22),
                                :server=>server)),
                        Seis.GeogStation{Float64}(net="IU", sta="COLA",
                            loc="20", cha="HNN", lat=64.873599, lon=-147.8616, elev=200,
                            dep=0, azi=0, inc=90, meta=Dict(
                                :sensor_description=>"Kinemetrics FBA-23 Low-GainSensor",
                                :scale=>53687.1, :scale_frequency=>1.0,
                                :scale_units=>"M/S**2", :sample_rate=>80.0,
                                :startdate=>DateTime(2005, 9, 28, 22),
                                :server=>server))]
                    for (o, ot) in zip(out, out_true)
                        # `collect(propertynames(o))` for compatibility with Julia v1.2
                        @testset "Field: $f" for f in filter(x->x∉(:pos, :meta), collect(propertynames(o)))
                            @test getproperty(o, f) === getproperty(ot, f)
                        end
                        @testset "Key: $k" for (k, v) in o.meta
                            @test o.meta[k] === ot.meta[k]
                        end
                    end
                end
            end
        end

        @testset "StationXML format" begin
            sxml = StationXML.FDSNStationXML(source="Me", created=DateTime(1000),
                schema_version="1.1")
            net = StationXML.Network(code="AN", start_date=DateTime(2000), end_date=DateTime(3000))
            sta = StationXML.Station(code="FAKE",
                longitude=1, latitude=2, elevation=3, site=StationXML.Site(name="Fake site"),
                start_date=DateTime(2000), end_date=DateTime(2001))
            cha1 = StationXML.Channel(code="UHT", longitude=10, latitude=20, elevation=30,
                depth=40, location_code="", start_date=DateTime(2002), azimuth=100,
                dip=10)
            cha2 = StationXML.Channel(code="UHU", longitude=10, latitude=20, elevation=30,
                depth=0, location_code="00", start_date=DateTime(2003))

            @testset "Network level" begin
                request = FDSNStation(network="?N", level="network")
                sxml′ = deepcopy(sxml)
                push!(sxml′.network, net)
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Messages.Response(SUCCESS, Dict(),
                        body=string(StationXML.xmldoc(sxml′)))
                    out = SeisRequests.parse_station_response(T, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{T}}
                    @test length(out) == 1
                    @test out[1] == Seis.GeogStation{T}(net="AN", meta=Dict(
                            :startdate=>DateTime(2000),
                            :enddate=>DateTime(3000),
                            :stationxml=>sxml′,
                            :server=>server,
                            :request=>request))
                end
            end

            @testset "Station level" begin
                request = FDSNStation(network="?N", station="FAKE", level="station")
                sxml′ = deepcopy(sxml)
                push!(sxml′.network, deepcopy(net))
                push!(sxml′.network[1].station, deepcopy(sta))
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Response(SUCCESS, Dict(),
                        body=string(StationXML.xmldoc(sxml′)))
                    out = SeisRequests.parse_station_response(T, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{T}}
                    @test length(out) == 1
                    @test out[1] == Seis.GeogStation{T}(net="AN", sta="FAKE",
                        lon=1, lat=2, elev=3, meta=Dict(
                            :startdate=>DateTime(2000), :enddate=>DateTime(2001),
                            :stationxml=>sxml′, :server=>server, :request=>request))
                end
            end

            @testset "Channel level" begin
                request = FDSNStation(network="?N", station="FAKE", channel="*",
                    level="channel")
                sxml′ = deepcopy(sxml)
                push!(sxml′.network, net)
                push!(sxml′.network[1].station, deepcopy(sta))
                append!(sxml′.network[1].station[1].channel, deepcopy([cha1, cha2]))
                @testset "Eltype $T" for T in (Float32, Float64)
                    response = HTTP.Messages.Response(SUCCESS, Dict(),
                        body=string(StationXML.xmldoc(sxml′)))
                    out = SeisRequests.parse_station_response(T, request, response,
                        server)
                    @test out isa Vector{Seis.GeogStation{T}}
                    @test length(out) == 2
                    # Remove the other channel from the StationXML, as done by
                    # SeisRequests.filter_stationxml
                    sxml′_filtered = deepcopy(sxml′)
                    pop!(sxml′_filtered.network[1].station[1].channel)
                    @test out[1] == Seis.GeogStation{T}(net="AN", sta="FAKE",
                        loc="", cha="UHT", lon=10, lat=20, elev=30, dep=0.04,
                        azi=100, inc=100,
                        meta=Dict(:startdate=>DateTime(2002), :stationxml=>sxml′_filtered,
                            :server=>server, :request=>request))
                end
            end
        end
    end

    @testset "Real requests" begin
        @testset "Verbosity" begin
            @test_logs (:info, "Request status: Successful request, results follow"
                ) get_stations(network="GB", station="J*", channel="?HZ",
                    level="station", verbose=true)
            @test_logs get_stations(network="GB", station="JSA", level="station",
                verbose=false)
        end

        @testset "Eltype" begin
            @testset "$T" for T in (Float32, Float64)
                @test get_stations(network="GB", station="J*", channel="?HZ",
                    level="channel", format="text", T=T,
                    verbose=false) isa Vector{Seis.GeogStation{T}}
            end
            @testset "Default eltype" begin
                @test get_stations(network="GB", level="network",
                    verbose=false) isa Vector{Seis.GeogStation{Float64}}
            end
        end

        # Request with only keyword arguments forwarded to FDSNStation
        @testset "Simple" begin
            stas = get_stations(network="GB", station="J*", channel="?HZ",
                    level="channel", format="text", verbose=false)
            @test length(stas) == 2
            @test stas[1] == Seis.GeogStation{Float64}(net="GB", sta="JSA", loc="",
                cha="BHZ", lon=-2.171698, lat=49.187801, elev=39.0,
                azi=0, inc=0, dep=0, meta=Dict(
                    :startdate=>DateTime(2007, 9, 6),
                    :sensor_description=>"TR-240", :sample_rate=>50.0,
                    :scale_units=>"M/S", :scale_frequency=>1,
                    :scale=>4.79174e8, :server=>"IRIS"))
        end

        @testset "Event" begin
            evt = Seis.Event(lon=-90, lat=45, time=DateTime(2000))
            @testset "No period" begin
                stas = get_stations(evt, network="US", maxradius=10, verbose=false)
                stas′ = get_stations(network="US", longitude=evt.lon, latitude=evt.lat,
                    maxradius=10, starttime=evt.time, endtime=evt.time, verbose=false)
                @test length(stas) == length(stas′)
                # Use custom comparison function because StationXML 'created' date
                # will be different regardless
                for (s, s′) in zip(stas, stas′)
                    @test stations_equal(s, s′)
                end
                # Stations definitely were active when the event occurred
                @test all(x -> x.meta.startdate <= evt.time <=
                    coalesce(x.meta.enddate, typemax(DateTime)), stas)
                # Stations are within 10 degrees of event
                @test all(x -> distance_deg(evt, x) <= 10, stas)
            end

            @testset "Date/time range" begin
                stas = get_stations(evt, -0.1, 60, network="US", maxradius=10,
                    verbose=false)
                stas′ = get_stations(evt, -Millisecond(100), Second(60), network="US",
                    maxradius=10, verbose=false)
                @test all(stations_equal.(stas, stas′))
            end
        end
    end
end


text:
Base.@kwdef mutable struct ToyEnv <: AbstractEnvironment
    s::Int = 0
    a::Int = 0
    t::Int = 0
    max_length::Int = typemax(Int)
    step_hook = identity
    reward_scale::Float64 = 1
end

LyceumBase.statespace(::ToyEnv) = VectorShape(Int, 1)
LyceumBase.getstate!(s, e::ToyEnv) = s .= e.s
LyceumBase.setstate!(e::ToyEnv, s) = (e.s = s[]; e)

LyceumBase.observationspace(::ToyEnv) = VectorShape(Int, 1)
LyceumBase.getobservation!(o, e::ToyEnv) = o .= e.t

LyceumBase.actionspace(::ToyEnv) = VectorShape(Int, 1)
LyceumBase.getaction!(a, e::ToyEnv) = a .= e.a
LyceumBase.setaction!(e::ToyEnv, a) = (e.a = a[]; e)

LyceumBase.getreward(s, a, o, e::ToyEnv) = s[] * e.reward_scale

function LyceumBase.reset!(e::ToyEnv)
    e.s = e.a = e.t = 0
    return e
end

function LyceumBase.randreset!(rng::AbstractRNG, e::ToyEnv)
    e.s = rand(rng, -5:5)
    e.t = 0
    return e
end

function LyceumBase.step!(e::ToyEnv)
    e.s += e.a
    e.t += 1
    e.step_hook(e)
    return e
end

LyceumBase.isdone(s, o, e::ToyEnv) = e.t >= e.max_length

Base.time(e::ToyEnv) = e.t
LyceumBase.timestep(e::ToyEnv) = e.dt

function busyloop(dt::Real)
    t0 = time()
    while time() - t0 < dt
        ccall(:jl_gc_safepoint, Cvoid, ())
    end
end


text:
module EquivalenceScales
    #
    # This module implements a set of standard-ish Equivalence Scales, mainly taken from
    #  https://www.gov.scot/binaries/content/documents/govscot/publications/statistics/2020/01/equivalence-scales/documents/equivalence-scales-rationales-uses-and-assumptions/equivalence-scales-rationales-uses-and-assumptions/govscot%3Adocument/paper%2Bdiscussing%2Bequivalence%2Bscales%2Band%2Bunderlying%2Bassumptions.pdf
    #
    # As part of this, it describes a simple interface that can be used by any data structure than
    # contains records of people with ages and relationships.
    #
    export
        EQ_P_Type,
        EQ_Person,
        EQScales,
        eq_dependent_child,
        eq_head,
        eq_other_adult,
        eq_spouse_of_head,
        get_equivalence_scales,
        get_age,
        eq_rel_to_hoh

    @enum EQ_P_Type eq_head eq_spouse_of_head eq_other_adult eq_dependent_child
    @enum Scales oecd oxford mcclements square_root per_capita
    struct EQ_Person
        age :: Int
        eqtype :: EQ_P_Type
    end

    struct EQScales{T<:Real}
        oxford :: T
        oecd_bhc :: T
        oecd_ahc :: T
        mcclements_bhc :: T
        mcclements_ahc :: T
        square_root :: T
        per_capita :: T
    end

    function eq_rel_to_hoh( p )::EQ_P_Type
        p.eqtype
    end

    function get_age( p )::Int
        p.age
    end

    """
    More Julian version, perhaps?
    """
    function onescale( T::Type, scale :: Scales, perss , before_hc :: Bool ) :: T
        s = zero(T)
        n = size(perss)[1]
        add = zero(T)
        eq = zero(T)
        if scale == per_capita
            return n
        elseif scale == square_root
            return sqrt(n)
        elseif scale in [oxford,oecd]
            for p in perss
                rel = eq_rel_to_hoh(p)
                if rel == eq_head
                    eq += 1
                elseif n == 1
                    eq += 1
                    println( "only 1 person; non head rel=$rel")
                else
                    if get_age(p) <= 14
                        add = if scale == oxford
                            0.5
                        elseif scale == oecd
                            before_hc ? 0.3 : 0.34
                        end
                    else
                        add = if scale == oxford
                            0.7
                        elseif scale == oecd
                            before_hc ? 0.5 : 0.72
                        end
                    end
                    eq += add
                end
            end # pers loop
            @assert eq >= 1 "eq is $eq num people $n scale $scale"
        elseif scale == mcclements
            num_extra_adults = 0
            for p in perss
                age = get_age(p)
                rel = eq_rel_to_hoh( p )
                if rel == eq_head
                    eq += 1
                elseif rel == eq_spouse_of_head
                    eq += before_hc ? 0.64 : 0.82
                elseif rel == eq_other_adult
                    num_extra_adults += 1
                    if num_extra_adults == 1
                        eq += before_hc ? 0.75 : 0.82
                    elseif num_extra_adults == 2
                        eq += before_hc ? 0.69 : 0.82
                    else
                        eq += before_hc ? 0.59 : 0.73
                    end
                elseif rel == eq_dependent_child
                    if age in 0:1
                        eq += before_hc ? 0.148 : 0.13
                    elseif age in 2:4
                        eq += before_hc ?  0.295 : 0.33
                    elseif age in 5:7
                        eq += before_hc ? 0.344 : 0.38
                    elseif age in 8:10
                        eq += before_hc ? 0.377 : 0.42
                    elseif age in 11:12
                        eq += before_hc ? 0.41 : 0.47
                    elseif age in 13:15
                        eq += before_hc ? 0.443 : 0.51
                    elseif age in 16:21
                        eq += before_hc ? 0.59 : 0.69
                    end
                end # dependent child
            end # mcclements
        end
        return eq
    end

    #
    # more julian? anything with get_age, get_relation_to_hoh
    #

    """
    Make a struct full of Equivalence Scales; see:
    https://www.gov.scot/binaries/content/documents/govscot/publications/statistics/2020/01/equivalence-scales/documents/equivalence-scales-rationales-uses-and-assumptions/equivalence-scales-rationales-uses-and-assumptions/govscot%3Adocument/paper%2Bdiscussing%2Bequivalence%2Bscales%2Band%2Bunderlying%2Bassumptions.pdf
    perss is some collection of things each of which can respond to
    `get_age` and `eq_get_rel`
    """
    function get_equivalence_scales( T :: Type, perss ) :: EQScales{T}
        v_oecd_bhc = onescale( T, oecd, perss, true )
        v_oecd_ahc = onescale( T, oecd, perss, false )
        v_oxford = onescale( T, oxford, perss, false ) # bhc irrelevant
        v_mcclements_bhc = onescale( T, mcclements, perss, true )
        v_mcclements_ahc = onescale( T, mcclements, perss, false )
        v_square_root = onescale( T, square_root, perss, false )
        v_per_capita = onescale( T, per_capita, perss, false )
        return EQScales{T}(
            v_oxford,
            v_oecd_bhc,
            v_oecd_ahc,
            v_mcclements_bhc,
            v_mcclements_ahc,
            v_square_root,
            v_per_capita )
    end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module TimeSeriesIO

# package code goes here
export TimeArray, DataFrame

using DataFrames
using TimeSeries
using DataStructures: CircularBuffer

import TimeSeries.TimeArray

include("timearray.jl")
include("dataframe.jl")

end # module


text:
"""
Run policy on a history loaded in HistoryReplayer and return a new history
"""
function replay(hrpl::HistoryReplayer, policy::Policy, rng::AbstractRNG)::SimHistory

    hrec = HistoryRecorder(max_steps = selling_period_end(hrpl), rng = rng)
    h = simulate(hrec, hrpl, policy)

    return h
end

"""
Get NamedTuple of metrics:
    :r revenue,
    :u utilization,
    :nₛ number of sold products,
    :nᵣ number of non-empty requests
 for given SimHistory
"""
function get_metrics(m::PMDP, h::AbstractSimHistory)::NamedTuple
    revenue = sum(h[:r])
    sold_product_sizes::Array{Int64,1} = [sum(product(m, e.s)) for e in h if e.s.c!=e.sp.c]
    n_sales = length(sold_product_sizes)
    n_requests = length([e for e in h if e.s.iₚ!=PMDPs.empty_product_id(m)])
    utilization = sum(sold_product_sizes)
    @assert utilization == sum(h[1].s.c - h[end].sp.c) string("utilization ", utilization, " is not ", h[1].s.c, " - ", h[end].sp.c, h)
    (r = revenue, u = utilization, nₛ = n_sales, nᵣ = n_requests)
end

"""
Return DataFrame of evaluation metrics from running given policies on a given history of requests.

In addition to the metrics, save the name of the policy, hash of the request sequence and seed used in the replay.

requests::AbstractSimHistory or similar
"""
function eval_policy(mdp::PMDP, requests::Union{AbstractSimHistory, AbstractArray{<:NamedTuple}}, policies::NamedTuple,
              rng::AbstractRNG)::DataFrame
    hrpl = HistoryReplayer(mdp, requests)

    metrics = DataFrame()

    for (name, policy) in pairs(policies)
        m = NamedTuple()
        error = nothing
        try
            h, stats... = @timed replay(hrpl, policy, rng)
            m = (;get_metrics(mdp, h)..., stats...)

        catch err
           error = err
        end
        push!(metrics, (m..., name=name, sequence=hash(requests), error=error, replay_rng_seed=rng.seed))
    end
    return metrics
end

"""
Return DataFrame of evaluation metrics of given tuple of policies on a sequence of request sequences.

requests_sequences::AbstractArray{<:AbstractSimHistory} or similar shape
"""
function eval_policy(mdp::PMDP, request_sequences::AbstractArray{<:AbstractSimHistory},
              policies::NamedTuple, rng::AbstractRNG)::DataFrame
    metrics = DataFrame()
    @showprogress 1 "Evaluating traces..." for sequence in request_sequences
        mₛ = eval_policy(mdp, sequence, policies, rng)
        metrics = vcat(metrics, mₛ, cols=:union)
    end
    return metrics
end


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function prepare_beta_vec(epmat::EPMat{T}, beta_vec::AbstractVector{T}) where T
    N = length(epmat.v)
    isempty(beta_vec) && return spzeros(T, N)
    beta_len = length(beta_vec)
    beta_len != N && error("Dim missmatch, beta_vec lenght != N, ($(beta_len)) != ($N)")
    return sparse(beta_vec)
end

function prepare_beta_vec(epmat::EPMatT0{T}, beta_vec::AbstractVector{T}) where T
    M = length(epmat.vd) # Number of dependent variables
    N = M + length(epmat.vi)
    isempty(beta_vec) && return spzeros(T, N)
    beta_len = length(beta_vec)
    beta_len != N && error("Dim missmatch, beta_vec lenght != N, ($(beta_len)) != ($N)")
    return sparse(beta_vec[epmat.usperm]) # see echelonize
end


text:
module GJRand

export GJRandRNG

using Random: Random, AbstractRNG, Sampler, SamplerType, RandomDevice

import Random: rand, seed!

mutable struct GJRandRNG <: AbstractRNG
    a::UInt64
    b::UInt64
    c::UInt64
    d::UInt64
    GJRandRNG(seed::UInt64) = seed!(new(), seed)
    GJRandRNG() = GJRandRNG(rand(RandomDevice(),UInt64))
end

function seed!(gjrand::GJRandRNG,seed::UInt64)
    gjrand.a=seed
    gjrand.b=UInt64(0)
    gjrand.c=UInt64(2000001)
    gjrand.d=UInt64(0)
    for j=1:14
        advance!(gjrand)
    end
    gjrand
end
function rotate(x::UInt64,k::UInt64)
        return ((x << k) | (x >> (64 - k)))
end

function advance!(gjrand::GJRandRNG)
        gjrand.b += gjrand.c;
        gjrand.a =  rotate(gjrand.a, UInt64(32))
        gjrand.c ⊻= gjrand.b;
        gjrand.d += UInt64(0x55aa96a5);
        gjrand.a += gjrand.b;
        gjrand.c =  rotate(gjrand.c, UInt64(23))
        gjrand.b ⊻= gjrand.a;
        gjrand.a += gjrand.c;
        gjrand.b =  rotate(gjrand.b, UInt64(19))
        gjrand.c += gjrand.a;
        gjrand.b += gjrand.d;
    end

function rand(gjrand::GJRandRNG, ::SamplerType{UInt64})
    advance!(gjrand);
    return UInt64(gjrand.a);
end
for T = [Bool, Base.BitInteger64_types...]
    T === UInt64 && continue
    @eval rand(rng::GJRandRNG, ::SamplerType{$T}) = rand(rng, UInt64) % $T
end
Random.rng_native_52(::GJRandRNG) = UInt64

end


text:
"""
	struct IndexSubGrid{G,I,T,N} <: AbstractSubGrid{T,N}

An IndexSubGrid is a subgrid corresponding to a certain range of indices of the
underlying grid.
It is assumed to be an 1D grid.
"""
struct IndexSubGrid{G,I,T,N} <: AbstractSubGrid{T,N}
	supergrid  :: G
	subindices :: I
	domain 	   :: Domain

	function IndexSubGrid{G,I,T,N}(supergrid::AbstractGrid{T,N}, subindices, domain=Interval(first(supergrid), last(supergrid))) where {G,I,T,N}
		@assert length(subindices) <= length(supergrid)
		new(supergrid, subindices, domain)
	end
end

IndexSubGrid(grid::AbstractGrid{T,N}, i, domain=Interval(first(grid), last(grid))) where {T,N} =
    IndexSubGrid{typeof(grid),typeof(i),T,N}(grid, i, domain)

name(g::IndexSubGrid) = "Index-based subgrid"

subindices(g::IndexSubGrid) = g.subindices

similar_subgrid(g::IndexSubGrid, g2::AbstractGrid) = IndexSubGrid(g2, subindices(g))

length(g::IndexSubGrid) = length(subindices(g))

size(g::IndexSubGrid) = (length(g),)

eachindex(g::IndexSubGrid) = eachindex(subindices(g))

# The speed of this routine is the main reason why supergrid and subindices
# are typed fields, leading to extra type parameters.
unsafe_getindex(g::IndexSubGrid, idx) = unsafe_getindex(g.supergrid, g.subindices[idx])

function mask(g::IndexSubGrid)
    mask = falses(size(supergrid(g)))
	for i in g.subindices
		mask[i] = true
	end
    mask
end

support(g::IndexSubGrid) = g.domain

# Check whether element grid[i] (of the underlying grid) is in the indexed subgrid.
issubindex(i, g::IndexSubGrid) = in(i, subindices(g))

# getindex(grid::AbstractGrid, i::Range) = IndexSubGrid(grid, i)

getindex(grid::AbstractGrid, i::AbstractArray{Int}) = IndexSubGrid(grid, i)


text:
# ---
# title: 1345. Jump Game IV
# id: problem1345
# author: zhwang
# date: 2022-01-22
# difficulty: Hard
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/jump-game-iv/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`, you are initially positioned at the first
# index of the array.
#
# In one step you can jump from index `i` to index:
#
#   * `i + 1` where: `i + 1 < arr.length`.
#   * `i - 1` where: `i - 1 >= 0`.
#   * `j` where: `arr[i] == arr[j]` and `i != j`.
#
# Return _the minimum number of steps_ to reach the **last index** of the array.
#
# Notice that you can not jump outside of the array at any time.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [100,-23,-23,404,100,23,23,23,3,404]
#     Output: 3
#     Explanation: You need three jumps from index 0 --> 4 --> 3 --> 9. Note that index 9 is the last index of the array.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [7]
#     Output: 0
#     Explanation: Start index is the last index. You don't need to jump.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [7,6,9,6,9,6,9,7]
#     Output: 1
#     Explanation: You can jump directly from index 0 to index 7 which is last index of the array.
#
#
# **Example 4:**
#
#
#
#     Input: arr = [6,1,9]
#     Output: 2
#
#
# **Example 5:**
#
#
#
#     Input: arr = [11,22,7,7,7,7,7,7,7,22,13]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 5 * 10^4`
#   * `-10^8 <= arr[i] <= 10^8`
#
#
## @lc code=start
using LeetCode

function min_jumps_1345(arr::Vector{Int})::Int
    ## data initialize
    n, step = length(arr), 0
    nums = Dict{Int,Set}(x => Set{Int}() for x in Set{Int}(arr))
    for (i, val) in enumerate(arr)
        push!(nums[val], i)
    end
    ## BFS Search
    up_nums, up_inds, inds = Set{Int}([arr[1]]), Set{Int}([1]), Set{Int}([0, 1, n + 1])
    while !(n in up_inds)
        down_inds = Set{Int}()
        ## indexs by numbers
        for num in up_nums
            union!(down_inds, nums[num])
        end
        ## nearby indexs
        for i in up_inds
            push!(down_inds, i - 1)
            push!(down_inds, i + 1)
        end
        ## remove extra indexs
        for i in down_inds
            i in inds && delete!(down_inds, i)
        end
        up_inds, up_nums = down_inds, Set{Int}(arr[i] for i in down_inds)
        union!(inds, up_inds)
        step += 1
    end
    return step
end

## @lc code=end


text:
# a simple LP example
#
#   maximize x
#
#   s.t. x >= 2
#        x <= 1
#        x >= 0
#
#   conflict: the two constraints

@testset "IIS" begin
    env = CPLEX.Env()

    model = CPLEX.Model(env, "iis")
    CPLEX.set_sense!(model, :Max)

    # add variables
    CPLEX.add_var!(model, 1.0, 0., Inf)  # x

    # add constraints
    CPLEX.add_constr!(model, [1.], '>', 2.)
    CPLEX.add_constr!(model, [1.], '<', 1.)

    # compute conflict
    c = CPLEX.c_api_getconflict(model)

    @test c !== nothing
    @test c.stat == CPLEX.CPX_STAT_CONFLICT_MINIMAL
    @test c.nrows == 2
    @test c.rowind[1] == 0
    @test c.rowind[2] == 1
    @test c.ncols == 0

    # gc()  # test finalizers
end


text:
const Field = Composable
const ComposableLike = Union{Composable, Tuple{Vararg{Pair}}}

"""
    UsageError(binding, value)
    UsageError(message::String, value)

Exception type for incorrect usage of some type or function `binding`.
Defaults to throwing the docstring as a message, but you can supply a
custom `message`. Any arguments can be supplied as `value`.
"""
immutable UsageError <: Exception
    message::String
    value
    function UsageError(message::String, value)
        new(message, value)
    end
end
function UsageError(binding, value)
    message = stringmime("text/plain", Docs.doc(binding))
    UsageError(message, value)
end

function Base.showerror(io::IO, e::UsageError)
    println(io, "Usage Error:")
    for line in split(e.message, r"\r?\n")
        println(io, "   ", line)
    end
    print(io, "Supplied value: ", e.value)
end

"""
    Fields(typ)

Gets a tuple of all field types of a composable
"""
Fields(x::Type{<:Composable}) = ()
Fields(x::Composable) = Fields(typeof(x))

function field2string(field)
    tname = typename(field)
    mod = tname.module
    fname = tname.name
    mname = module_name(mod)
    # this is kinda terrible.. But we need nice names, and uniqueness of keyes.
    mname_str = mod == Main || module_parent(mod) == Main ? "" : string(mname, "_")
    string(mname_str, fname)
end
function field2symbol(field)
    name_str = field2string(field)
    Symbol(lowercase(name_str))
end

"""
    @needs composable: Field1, Field2, ...
    @needs composable: Field1, name2 = Field2, ...

Asserts that a `composable` has the listed fields. Optionally assign a
variable binding to chosen fields; this variable has scope within the
containing block (e.g., function definition).
"""
macro needs(expr::Expr)
    comp_fields = @match expr begin
        (composable_:f1_)  => (composable, Any[f1])
        (composable_:f1_, ftail__)  => (composable, Any[f1, ftail...,])
        (composable_:var1_ = field1_)  => (composable, Any[:($var1 = $field1),])
        (composable_:var1_ = field1_, ftail__)  => (composable, Any[:($var1 = $field1), ftail...,])
    end
    if comp_fields == nothing
        throw(UsageError("Usage: `@needs composable: Field1, name = Field2, ...`", expr))
    end
    composable, fields = comp_fields
    checks = Expr(:block)
    checks.args = map(fields) do field
        field_expr = @match field begin
            (var_ = field_) => begin
                quote
                    haskey($composable, $field) || error("$($composable) doesn't contain field $($field)")
                    $var = $composable[$field]
                end
            end
            field_ => begin
                :(haskey($((composable)), $((field))) || error("$($composable) doesn't contain field $($field)"))
            end
        end
        if field_expr == nothing
            throw(UsageError("Usage: `@needs composable: Field1, name = Field2, ...`", expr))
        end
        field_expr
    end
    esc(checks)
end

"""
An partially initialized composed type for type Parent
"""
immutable Partial{Parent, T}
    val::T
end
Partial{P, T}(::Type{P}, x::T) = Partial{P, T}(x)
Partial{P}(::Type{P}) = Partial(P, ComposedDict())

"""
    default(FieldType)
    default(parent, FieldType)

Construct a default for `FieldType`.  Can be specialized for different
`parent` objects.  If you do not specialize it for type `FieldType`, this
calls `FieldType()`. Values optionally holds a partially initialised Partial object, e.g.
when called like Partial(A => x, B => c), A and B will be supplied via values.
This way we can construct better defaults!
"""
function default{P <: ComposableLike, F <: Field}(
        ::Type{F}, ::Type{P},
    )
    default(F, Partial(P))
end

default{F <: Field}(::Type{F}, ::Partial) = default(F)
default{F <: Field}(::Type{F}) = F()

"""
    convert(ComposableType, FieldType, x)

Converts field `FieldType` in `x` to the type of field `FieldType` in
`ComposableType`. `x` may be supplied as another composable (thus
extracting `x[FieldType]`), or as a tuple of `FieldType => value` pairs.
"""
function convertfor{F <: Field, P}(::Type{F}, x::Partial{P}, value)
    convert(cfieldtype(P, F), value)
end

"""
    fieldindex(ComposableType, FieldType)

Calculates the index into a struct type from a Composable type. The
return is a tuple of `Val{i}` instances, where `i` is the index of
`FieldType` in `ComposableType`. If `FieldType` is not a field of
`ComposableType`, then `i` is 0.

# Example

```jldocstring
@composed type Composed
    Scale
    Position
end
fieldindex(Composed, Position) == (Val{2}(),)
```
"""
fieldindex{T <: Composable, F <: Field}(::Type{T}, ::Type{F}) = Val{0}()
fieldindex{T <: Composable, F <: Field}(::T, ::Type{F}) = fieldindex(T, F)

Base.@inline haskey{T <: Composable, N}(::Type{T}, field::Val{N}) = N > 0
Base.@pure haskey{T <: Composable, F <: Field}(c::Type{T}, field::Type{F}) = haskey(T, fieldindex(T, F))
Base.@pure haskey{T <: Composable, F <: Field}(c::T, field::Type{F}) = haskey(T, fieldindex(T, F))

# Would like to extend Base.fieldtype, but that is an Builtin function which can't
# be extended.
"""
    cfieldtype(ComposableType, FieldType)

Return the type of `FieldType` in a composed type `ComposedType`. The default
is `Any`. This is automatically specialized by `@composed`.
"""
cfieldtype{T <: Composable, F <: Field}(ct::T, ::Type{F}) = cfieldtype(T, fieldindex(T, F), F)
cfieldtype{T <: Composable, F <: Field}(ct::Type{T}, ::Type{F}) = cfieldtype(T, fieldindex(T, F), F)
function cfieldtype{T <: Composable, N}(::Type{T}, ::Val{N}, field = N)
    if N === 0
        error("Can't get $field from type $T")
    else
        fieldtype(T, N)
    end
end
"""
    fieldconvert(ComposableType, FieldType, value)

Converts `value` to the type needed for field `FieldType` in type `ComposableType`.
"""
function fieldconvert{T <: Composable}(ct::T, field, value)
    convert(cfieldtype(T, field), value)
end

@propagate_inbounds function _setindex!{N}(ct::Composable, val, field::Val{N})
    setfield!(ct, N, fieldconvert(ct, field, val))
end

@propagate_inbounds function _setindex!{F <: Field}(ct::Composable, val, ::Type{F})
    idx = fieldindex(ct, F)
    if !isa(idx, Val{0})
        _setindex!(ct, val, idx)
    else
        throw(BoundsError(ct, F))
    end
end

@propagate_inbounds function setindex!{F <: Field}(ct::Composable, value, field::Type{F})
    _setindex!(ct, value, field)
end
@propagate_inbounds function getindex{N}(ct::Composable, field::Val{N})
    getfield(ct, N)
end
@propagate_inbounds function getindex{F <: Field}(ct::Composable, ::Type{F})
    idx = fieldindex(ct, F)
    if !isa(idx, Val{0})
        return getindex(ct, idx)
    else
        throw(BoundsError(ct, F))
    end
end

Base.length(x::Composable) = length(Fields(x))
Base.start(x::Composable) = (Fields(x), 1)
function Base.next(x::Composable, state)
    fields, idx = state
    field = fields[idx]
    (field => x[field], (fields, idx + 1))
end
Base.done(x::Composable, state) = state[2] > length(state[1])

"""
Default empty constructor.
Will initialize fields with `default(T, field)`
"""
function (t::Type{T}){T <: Composable}()
    fields = Fields(T)
    if isempty(fields) # we're at a leaf field without an empty constructor defined
        # TODO think of good error handling, that correctly advises the user
        error("No default for $t")
    end
    p = Partial(T)
    T(map(field-> default(field, p), fields)...)
end

# Constructor from partial composed type (tuple of pairs)
"""
Partially initialization constructor.
You can supply values for fields as pairs, and the rest will be filled in by
`default(T, field, pairs)`.
"""
function (::Type{T}){T <: Composable}(c::ComposableLike)
    partial1 = Partial(T, c)
    conv_fields = ComposedDict(Dict{DataType, Any}(map(c) do fieldval
        Field, val = fieldval
        if haskey(T, Field)
            Field => convertfor(Field, partial1, val) # converts might be expensive or a noop, so lets do it first
        else
            fieldval
        end
    end))
    partialconv = Partial(T, conv_fields)
    fields = map(Fields(T)) do Field
        x = if haskey(conv_fields, Field)
            get(conv_fields, Field)
        else
            default(Field, partialconv)
        end
        if !haskey(partialconv, Field)
            partialconv[Field] = x
        end
        x
    end
    T(fields...)
end

"""
    @field FieldType

Creates a type `FieldType` as a component of a Composable type.
"""
macro field(expr)
    usage = """
        "Must be @field Sym [<: OptionalSuperType] = default_body.
        Found: $expr
    """
    name_supertyp_default = @match expr begin
        (name_ <: supertyp_ = default_) => name, supertyp, default
        (name_ <: supertyp_) => name, supertyp, nothing
        (name_ = default_) => name, FieldTraits.Field, default
        (name_) => name, FieldTraits.Field, nothing
    end
    if name_supertyp_default == nothing
        error(usage)
    end
    name, supertyp, default = map(esc, name_supertyp_default)
    default_expr = if default != nothing
        :(FieldTraits.default{T <: $name}(::Type{T}) = $default)
    else
        :()
    end
    quote
        Base.@__doc__ immutable $name <: $supertyp end
        $default_expr
    end
end

function add_fields!(composedfields, result, fields, T_args, ftype_funcs, idx = 1)
    for field in composedfields
        sym, Field, T = field, field, field
        if is_linenumber(field)
            #push!(result, field)
            continue
        # untyped field which will get a parameter
        elseif isa(field, Symbol) || isa(field, GlobalRef) ||
                (isa(field, Expr) && field.head == :(.))
            Field = eval(current_module(), field)
            sym = Symbol(field2string(Field))
            T = sym
            push!(T_args, Field => T)
        elseif isa(field, Expr) # typed field with no parameter
            if field.head == :(::)
                sym, T = field.args
                Field = eval(current_module(), sym)
                push!(ftype_funcs, Field => T)
            elseif field.head == :(.) || isa(field, GlobalRef) # untyped field which will get a parameter
                push!(T_args, field => field)
                Field = eval(current_module(), field)
            else
                match = @match field begin
                    (<: Composable_) => begin
                        T = eval(current_module(), Composable)
                        idx = add_fields!(Fields(T), result, fields, T_args, idx)
                        continue
                    end
                end
                match == nothing && error("Must be `field`, `field::T`, or `<: T`. Found: $field")
            end
        elseif isa(field, DataType) || (isdefined(Base, :UnionAll) && isa(field, UnionAll))
            sym = Symbol(field2string(field))
            T = sym; Field = field
            push!(T_args, field => T)
        else
            error("Unsupported expr: $field $(typeof(field))")
        end
        fname = field2symbol(Field)
        push!(fields, Field)
        push!(result, :($fname::$T))
        idx += 1
    end
    idx
end

"""
Usage:
```Julia
@composed type Name [<: OptionalSuperType]
    <: InheritFieldsFromThisType
    FieldTrait1
    FieldTrait2::OptionalStrictType
end
```
"""
macro composed(expr::Expr)
    composed_type(expr)
end

function composed_type(expr::Expr, additionalfields = [], supertyp = Composable)
    name_supertyp_fields, mutable = @match expr begin
        type name_ <: supertyp_
            fields__
        end => (name, supertyp, fields), true
        type name_
            fields__
        end => (name, supertyp, fields), true

        immutable name_ <: supertyp_
            fields__
        end => (name, supertyp, fields), false
        immutable name_
            fields__
        end => (name, supertyp, fields), false
    end
    if name_supertyp_fields == nothing
        throw(UsageError(composed, expr))
    end
    name, supertyp, fields = name_supertyp_fields
    composedfields = [additionalfields; fields...]
    typedfields = []; fields = []; FT_args = []; ftype_funcs = []
    add_fields!(composedfields, typedfields, fields, FT_args, ftype_funcs)
    idxfuncs = Expr(:block)
    for (i, Field) in enumerate(fields)
        push!(idxfuncs.args,
            :(FieldTraits.fieldindex(::Type{<: $name}, ::Type{$Field}) = Val{$i}())
        )
    end
    T_args = last.(FT_args)
    fieldfuncs = quote
        FieldTraits.Fields(::Type{$name}) = ($(fields...),)
    end
    typ = Expr(
        :type, mutable,
        :($(name){$(T_args...)} <: $supertyp),
        Expr(:block, typedfields...)
    )
    expr = quote
        Base.@__doc__ $typ
        $(fieldfuncs)
        $(idxfuncs)
    end
    esc(expr)
end

# A dictionary wrapper supporting the ComposedApi
immutable ComposedDict{V} <: Composable
    data::Dict{DataType, V}
end
ComposedDict() = ComposedDict(Dict{DataType, Any}())

Base.start(x::ComposedDict) = start(x.data)
Base.next(x::ComposedDict, state) = next(x.data, state)
Base.done(x::ComposedDict, state) = done(x.data, state)
Base.length(x::ComposedDict) = length(x.data)

function getindex{F <: Field}(cd::ComposedDict, ::Type{F})
    # TODO search!
    getindex(cd.data, F)
end

function setindex!{F <: Field}(cd::ComposedDict, val, ::Type{F})
    # TODO search!
    setindex!(cd.data,val, F)
end
function setindex!{F <: Field}(cd::Partial, val, ::Type{F})
    # TODO search!
    setindex!(cd.val, val, F)
end
function get!{F <: Field}(cd::ComposedDict, ::Type{F}, default)
    get!(cd.data, F, default)
end

Fields(::Type{ComposedDict}) = ()
Fields(::ComposedDict) = ()

# helper to use tuples as partial Composed Traits
function Base.getindex{F <: Field}(x::Tuple{}, ::Type{F})
    error("Tuple doesn't have field $F")
end
# @pure be legal and helpful here?
@propagate_inbounds function Base.getindex{N, F <: Field}(x::Tuple{Vararg{Pair, N}}, ::Type{F})
    a, b = first(x)
    a <: F && return b
    getindex(Base.tail(x), F)
end

haskey{F <: Field}(cd::ComposedDict, ::Type{F}) = haskey(cd.data, F)

function haskey{F <: Field}(x::Partial, ::Type{F})
    haskey(x.val, F)
end
haskey{F <: Field}(x::Tuple{}, ::Type{F}) = false

function haskey{N, F <: Field}(x::Tuple{Vararg{Pair, N}}, ::Type{F})
    a, b = first(x)
    a <: F && return true
    haskey(Base.tail(x), F)
end

##################################################################
# get
#

# implements fallback to default
# implements fallback to default
function get{F <: Field}(p::Partial, ::Type{F})
    get(()-> default(F, p), p.val, F)
end
function get!{F <: Field}(cd::Partial, ::Type{F}, default...)
    x = get(cd, F, default...)
    cd.val[F] = x
    x
end

function get{T <: ComposableLike, F <: Field}(x::T, ::Type{F})
    get(()-> default(F, Partial(T, x)), x, F)
end
# implements to return a default value
function get{F <: Field}(x::ComposableLike, ::Type{F}, default)
    get(()-> default, x, F)
end

get{F <: Field}(Func, x::Tuple{}, ::Type{F}) = Func()
function get{N, F <: Field}(Func, x::Tuple{Vararg{Pair, N}}, ::Type{F})
    a, b = first(x)
    a <: F && return b
    get(Func, Base.tail(x), F)
end

function get{F <: Field}(Func, x::Composable, ::Type{F})
    haskey(x, F) && return x[F]
    Func()
end

function (==)(x::Composable, y::Composable)
    f1, f2 = Fields(x), Fields(y)
    f1 != f2 && return false
    for f in f1
        x[f] == y[f] || return false
    end
    true
end


text:
#
# GENERALIZED ARGMAX FUNCTION
# Implements the EVT-I business...

errmsg = "Please specify a log sum behavior. Usage is
`generalized_argmax!(y,x,sig, SumNormalize)` or
`generalized_argmax!(y,x,sig,SumRaw)`"

generalized_argmax!(y,x,sig)=error(errmsg)

function generalized_argmax!(y, x, sig, logsumspec)
    # `logsumspec` controls whether we subtract by greatest element of x or not
    # when computing stuff

    expsum = vexpsum(x, sig, logsumspec)  # Note: might be a number or a tuple
    # (depending on `logsumspec`)

    @inbounds for i in eachindex(x)
        y[i] = donorm(x[i], expsum, sig, logsumspec)
    end

    return genlog(expsum, sig, logsumspec)
end

struct SumNormalize end
struct SumRaw end

genlog(ss,sig,::Type{SumNormalize})=ss[2]+sig*log(ss[1])
genlog(ss,sig,::Type{SumRaw})=sig*log(ss)

vexpsum(x, sig, ::Type{SumNormalize})=vexpsum_norm(x,sig)
vexpsum(x, sig, ::Type{SumRaw})=vexpsum_raw(x,sig)

function vexpsum_norm(x::AbstractArray, sig)
    s = zero(eltype(x))
    maxel = maximum(x)
    @inbounds for i in eachindex(x)
        xnorm = x[i]-maxel
        s += exp(xnorm/sig)
    end
    return s, maxel
end

function vexpsum_raw(x::AbstractArray, sig)
    s = zero(eltype(x))
    @inbounds for i in eachindex(x)
        s += exp(x[i]/sig)
    end
    return s
end

donorm(xx, ss, sig, ::Type{SumNormalize})=exp((xx - ss[2])/sig)/ss[1]
donorm(xx, ss, sig, ::Type{SumRaw})=exp(xx/sig)/ss


text:
# The following reproduces the example introduced in Section 3.2.3 of Lütkepohl (2006)

using VectorAutoregressions
using CSV, DataFrames
using Dates

exdir = pkgdir(VectorAutoregressions) * "/examples";
datadf = DataFrame(CSV.File(exdir * "/lutkepohl_data.csv"));
first(datadf, 6)

# Take the subset from 1960q1 to 1978q4, and use log-difference for each series
datadf = datadf[datadf.date.<Date(1979, 3, 1), :];
Y = DataFrame(
    investment = diff(log.(datadf.invest)),
    income = diff(log.(datadf.income)),
    consumption = diff(log.(datadf.cons)),
);
first(Y, 6)

# ## Calculate lag selection criteria
ls = LagSelection(Y, 8)
# ## Estimate the VAR
v = VAR(Y, ls.selection["AIC"])
# ## Check stability of the VAR
sc = StabilityCheck(v)
# ## Check for autocorrelation in the residuals
pt = PortmanteauTest(v, 12)
bg1 = LMCorrTest(v, 4)
bg2 = LMCorrTest(v, 4, smallsample = true)


text:
# This file was generated, do not modify it. # hide
ŷ = predict(clf, rows=test)

auc(ŷ, y[test])


text:
"""
Example usage of Agtor which distributes runs across available processors.
"""

using Dates, DataFrames, CSV
using Agtor
using Distributed, FileIO, Glob
using JLD2

using Profile, BenchmarkTools, OwnTime, Logging

addprocs(2, exeflags="--project=.")

@everywhere begin

    using Dates, DataFrames, CSV

    using Agtor

    # Start month/day of growing season
    # Represents when water allocations are announced.
    const gs_start = (5, 1)

    """
    Allocation callback for example.

    Any interactions with other models could be defined in these callbacks
    """
    function allocation_callback!(zone, dt_i)
        # Resetting allocations for example run
        if monthday(dt_i) == gs_start
            # Example annual water allocations
            allocs = (surface_water=150.0, groundwater=40.0)

            update_available_water!(zone, allocs)
        end
    end

end

"""
Model setup for example run.
"""
function setup_zone(data_dir::String="test/data/")
    zone_spec_dir::String = "$(data_dir)zones"
    zone_params::Dict{Symbol, Dict} = load_spec(zone_spec_dir)

    collated_specs::Array = []
    agparams = collect_agparams!(zone_params[:Zone_1], collated_specs; ignore=[:crop_spec])

    climate_data::String = "$(data_dir)climate/farm_climate_data.csv"
    climate::Climate = load_climate(climate_data)

    return create(zone_params[:Zone_1], climate), agparams
end

"""
Run example scenarios (for a single farm area) by distributing across available cores.
Results are saved to a file on completion, based on scenario id.
"""
function test_scenario_run(data_dir::String="test/data/", result_dir::String="")::Nothing
    z1, agparams = setup_zone(data_dir)
    z1.manager = BaseManager("test")

    scen_data = joinpath(data_dir, "scenarios", "sampled_params.csv")
    samples = DataFrame!(CSV.File(scen_data))

    # Each row represents a scenario to run.
    # We distribute these across available cores.
    @sync @distributed (hcat) for row_id in 1:nrow(samples)
        tmp_z = deepcopy(z1)
        update_model!(tmp_z, samples[row_id, :])
        res = run_model(tmp_z, run_timestep!, allocation_callback!)

        pth = joinpath(result_dir, "sampled_params_batch_run_distributed_$(row_id).jld2")
        save_results!(pth, string(row_id), res)
    end

    # Collate all results into a single JLD2 file
    fn_pattern = joinpath(result_dir, "sampled_params_batch_run_distributed_*.jld2")
    collate_results!(fn_pattern, "dist_collated.jld2")

    return
end

@btime test_scenario_run()

# Loading and displaying results
saved_results = load("dist_collated.jld2")
@info saved_results

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
function rowselect(df::DataFrame, s::Expr)
  opera = String(s.args[1])
  if !occursin(".", opera)
    opera = "."*opera
  end
  opera = Symbol(opera)
  obs = s.args[2]
  val = s.args[3]
  df[eval(Expr(:call, opera, df[!, obs], val)), :]
end


text:
@testset "Cubes" begin
  cube_json_string = JSON.parse("""
                                {
                                  "@odata.context": "\$metadata#Cubes(Dimensions)/\$entity",
                                  "Name": "Test",
                                  "Rules": "SKIPCHECK;\\n\\n#['test_dim':'test1'] = N: 1;",
                                  "LastSchemaUpdate": "2020-10-09T06:15:58.169Z",
                                  "LastDataUpdate": "2020-10-09T06:15:58.169Z",
                                  "Attributes": {
                                    "Caption": "Test"
                                  },
                                  "Dimensions": [
                                    {
                                      "Name": "Dim 1",
                                      "UniqueName": "[Dim 1]",
                                      "Attributes": {
                                        "Caption": "Dim 1"
                                      }
                                    },
                                    {
                                      "Name": "Dim 2",
                                      "UniqueName": "[Dim 2]",
                                      "Attributes": {
                                        "Caption": "Dim 2"
                                      }
                                    }
                                  ]
                                }
                                """)

  dim_1 = JSON.parse("""
                    {
                      "Name": "Dim 1",
                      "UniqueName": "[Dim 1]",
                      "Attributes": {
                        "Caption": "Dim 1"
                      }
                    }
  """)

  dim_2 = JSON.parse("""
                    {
                      "Name": "Dim 2",
                      "UniqueName": "[Dim 2]",
                      "Attributes": {
                        "Caption": "Dim 2"
                      }
                    }
  """)

  attributes = Dict("Caption" => "Test")

  cube_result = Cube(
    "Test",
    "SKIPCHECK;\n\n#['test_dim':'test1'] = N: 1;",
    Dates.DateTime("2020-10-09T06:15:58.169"),
    Dates.DateTime("2020-10-09T06:15:58.169"),
    attributes,
    [Dimension(dim_1), Dimension(dim_2)],
    nothing,
    nothing,
    nothing,
    nothing,
    nothing,
  )

  @test cube_result.Name == "Test"
  @test cube_result.Rules == "SKIPCHECK;\n\n#['test_dim':'test1'] = N: 1;"
  @test cube_result.Dimensions isa Vector
  @test name(cube_result) == "Test"

  cube_single_arg_constructor = Cube("Another Test")
  @test name(cube_single_arg_constructor) == "Another Test"

  cube_json = Cube(cube_json_string)

  @test name(cube_json) == name(cube_result)
  @test cube_json == cube_result

  cube_as_dict = Cube(Dict(
    "LastDataUpdate" => "2020-10-09T17:30:34.063Z",
    "Rules" => nothing,
    "Attributes" => Dict{String,Any}("Caption" => "}Capabilities"),
    "LastSchemaUpdate" => "2020-10-09T17:30:34.063Z",
    "DrillthroughRules" => nothing,
    "Name" => "}Capabilities",
    "@odata.etag" => "W/\"98079d38abe094b7ecaaf37b3519525df0abb891\"",
  ))

  @test name(cube_as_dict) == "}Capabilities"

end


text:
module LinAlgUtils

using SparseArrays
using LinearAlgebra
using StaticArrays
using Rotations

export normalizecolumns
export invertSym3x3
export average_rotmat_l2
export evec_accuracy_report, check_linops, iseigenpairaccurate, eigensystemresiduals
export grassmanian_metric_oper, grassmanian_metric_angle, grassmanian_metric_chordal

# divide all columns of a sparse matrix by their l2 norm.
# Return the new matrix and the vector of column norms
function normalizecolumns(A)
    # get the column sums of A
    S = sqrt.(sum(abs2.(A), dims=1))
    S = reshape(S, (:))

    # get the nonzero entries in A.
    # ei is row index, ej is col index, ev is the value in A
    ei, ej, ev = findnz(A)

    # get the number or rows and columns in A
    m,n = size(A)

    # create a new normalized matrix. For each nonzero index (ei,ej), its new value will be
    # the old value divided by the sum of that column, which can be obtained by S[ej]
    A = sparse(ei,ej,ev./S[ej],m,n)
    return A, S
end

"""
Compute and return the inverse of a Symmetric 3x3 matrix
(this has been carefully profiled / tweaked for performance)
"""
function invertSym3x3(A::Array{Float64,2}) :: SMatrix{3,3, Float64}

    A = SMatrix{3,3}(A) # 6x speedup from this
    detinv = 1 / det(A)

    a = (A[2,2]*A[3,3] - A[2,3]*A[3,2]) * detinv
    b = (A[2,3]*A[3,1] - A[2,1]*A[3,3]) * detinv
    c = (A[2,1]*A[3,2] - A[2,2]*A[3,1]) * detinv

    #d = (A[1,3]*A[3,2] - A[1,2]*A[3,3]) * detinv   # it's symmetric
    e = (A[1,1]*A[3,3] - A[1,3]*A[3,1]) * detinv
    f = (A[1,2]*A[3,1] - A[1,1]*A[3,2]) * detinv

    #g = (A[1,2]*A[2,3] - A[1,3]*A[2,2]) * detinv
    #h = (A[1,3]*A[2,1] - A[1,1]*A[2,3]) * detinv
    i = (A[1,1]*A[2,2] - A[1,2]*A[2,1]) * detinv

    # building array at end: 10% speedup with StaticArrays, slowdown with
    # dynamic memory Arrays.
    return @SMatrix [a b c;
                     b e f;
                     c f i]
end

"""
Compute the average rotation matrix under the l2 geodesic distance. This
algorithm only is guaranteed to work in the rotations lie in a π/2 ball. (This
is generally the case because a π/2 ball covers all but a measure 0 subset of
SO(3).)
"""
function average_rotmat_l2(rots::Vector{RotMatrix{3, Float64}}) :: RotMatrix{3, Float64}
    # Add up all of the rotation matrices using an R^3x3 sum.
    Rsum = SMatrix{3,3}(rots[1])
    for k = 2:length(rots)
        Rsum = Rsum + rots[k]
    end

    # Do an SVD factorization
    fact = svd(Rsum, full=true)

    # form the average matrix
    Ravg = fact.U * fact.Vt

    # correct for cheirality
    if det(Ravg) < 0
        Ravg = fact.U * Diagonal([1.0, 1.0, -1.0]) * fact.Vt
    end

    return Ravg
end

function evec_accuracy_report(operator, evecs, evals; constraints=nothing)

    println("Eigen accuracy report.")
    println("    evecs: ", typeof(evecs), size(evecs))
    N = size(evecs, 2)
    ortho = norm(evecs'*evecs - Matrix{Float64}(I, N, N))
    println("    evecs ortho: |V'V-I|: ", ortho)

    if constraints !== nothing
        Nc = size(constraints, 2)
        println("    constraints: ", typeof(constraints), size(constraints))
        ortho = norm(constraints'*constraints - Matrix{Float64}(I, Nc, Nc))
        println("    constraints ortho: |C'C-I|: ", ortho)
    end

    for (k, (v, λ)) in enumerate(zip(eachcol(evecs), evals))
        v ./= norm(v)
        println("$k. --->")
        # If v is a large eval of covariance
        # it's a small eval of J'J.
        product = operator * v
        λ_actual = norm(product)

        product ./= λ_actual

        println("    stat1: |AX - λX| : ", norm(λ*v - λ_actual * product))
        println("    stat2: |AX| / λ  : ", λ_actual / λ )

        θ = rad2deg(acos(clamp(dot(v, product), -1, 1)))
        println("    stat3: θ deg     : ", θ)
        println("    stat4: dot       : ", dot(v, product))

        println("    supplied λ:    ", λ)
        println("    actual λ=|Av|: ", λ_actual)
        println("    actual λ=v'Av: ", v'*operator*v)

        if constraints !== nothing
            println("    nullspace mass |v'C|: $(norm(v'*constraints))")
        end
    end
end

"""
Perform just a small part of the evec_accuracy_report computation: check the
accuracy of a single eigenvector under the |Av-λv| score. Return True if
the accuracy is satisfactory.
"""
function iseigenpairaccurate(operator, v, λ; ϵ=1e-6) :: Bool
    v ./= norm(v)
    return norm(λ*v - operator*v) < ϵ
end

"""
Return the residuals of each given eigenpair under the |Av-λv| score.
Normalize by a reference scale (Fro norm / num elements) to get a relative measure.
"""
function eigensystemresiduals(operator, evecs, evals) :: Vector{Float64}
    residuals = []
    for (k, (v, λ)) in enumerate(zip(eachcol(evecs), evals))
        v ./= norm(v)
        push!(residuals,  norm(λ*v - operator*v))
    end
    reference = norm(operator, 2) / prod(size(operator))
    return residuals / reference
end

"""
Check if two linear operators are the same.

Allow for them to be in different forms (i.e., dense matrix, SparseMatrixCSC,
    LinearMap, ...). Just assume that they support size() and mult!()

Rather than going element by element (which would take forever, given that we
    don't even have element-wise access) do a randomized search. Generate
    random vectors, apply the operator, and take the norm in the range-space.

Returns the maximum norm-difference across all trials.

Throws if operator dimensions are not identical.
"""
function check_linops(operatorA, operatorB; numtrials=1000)
    # I guess size() doesn't exist for some LinearMaps?
    N = size(operatorA, 1)
    M = size(operatorA, 2)
    # N == size(operatorB, 1) || throw("The operator dimensions are not the same.")
    # M == size(operatorB, 2) || throw("The operator dimensions are not the same.")

    diffs = SharedArray{Float64}(numtrials)
    Threads.@threads for i = 1:numtrials
        x = rand(M)
        diffs[i] = norm(operatorA*x - operatorB*x)
    end
    return maximum(diffs)
end

"""
Measure the distance between two linear subspaces. V and W are both n-by-k
matrices representing elements of the Grassmanian Gr(k, mathbb{R}^n).

Orthogonality is *not* assumed: V represents the subspaces spanned by the
columns of V, and likewise for W.

Return the metric d(V, W) = || P_V - P_W ||
    where P_V is the canonical orthogonal projection onto the column space of V
    where P_W is the canonical orthogonal projection onto the column space of W
    where || ⋅ || is the operator norm
Because the orthogonal projectors are real symmetric, the operator norm is the
largest eigenvalue of the operator.
"""
function grassmanian_metric_oper(V, W)
    size(V) == size(W) || throw("Dimensions don't match!")

    # orthnormalize V and W using a skinny QR
    function orthonormalize(A)
        fact = qr(A)
        return Matrix(fact.Q)
    end
    QV = LinearMap(orthonormalize(V))
    QW = LinearMap(orthonormalize(W))

    # create a matrix-free linear operator for the projection operation
    # PV - PW = (QV*QV') - (QW*QW')
    D = QV * transpose(QV) - QW * transpose(QW)

    # run a matrix-free eigensolver to find the largest eigenvalue
    λ, ϕ, nconv, niter, nmult, resid = eigs(D, nev=1, which=:LM, maxiter=100)
    println("λ: ", λ)
    println("nconv: ", nconv)
    println("niter: ", niter)
    println("nmult: ", nmult)

    return abs(λ[1])
end

"""
Measure the distance between two linear subspaces. V and W are both n-by-k
matrices representing elements of the Grassmanian Gr(k, mathbb{R}^n).

Orthogonality is *not* assumed: V represents the subspace spanned by the
columns of V, and likewise for W.

Return the cosine of the largest principle angle:
        cos θ = max(v in V, w in W, |v|=1, |w|=1)  v'w

Procedure: orthonormalize V and W, then find largest singular value of V'W

Note: this version is a lot faster than the operator norm. I'm seeing about .4ms
    for a 2000x5 problem.
"""
function grassmanian_metric_angle(V, W)
    size(V) == size(W) || throw("Dimensions don't match!")

    # orthnormalize V and W using a skinny QR
    function orthonormalize(A)
        fact = qr(A)
        return Matrix(fact.Q)
    end
    QV = orthonormalize(V)
    QW = orthonormalize(W)

    # assume dims are small: find all singular values with a dense solver.
    # (if #columns is not small, it would be better to use an iterative solver.)
    σ = svdvals(QV' * QW) # returned in descending order
    return σ[1]
end

"""
Measure the distance between two linear subspaces. V and W are both n-by-k
matrices representing elements of the Grassmanian Gr(k, mathbb{R}^n).

Orthogonality is *not* assumed: V represents the subspace spanned by the
columns of V, and likewise for W.

Return the sum of the sins of the principal angles.

Procedure: orthonormalize V and W, then find the singular values of V'W

Citation: Huang, Qui, Calderbank: The Role of Principal Angles in
Subspace Classification
"""
function grassmanian_metric_chordal(V, W)
    size(V) == size(W) || throw("Dimensions don't match!")

    # orthnormalize V and W using a skinny QR
    function orthonormalize(A)
        fact = qr(A)
        return Matrix(fact.Q)
    end
    QV = orthonormalize(V)
    QW = orthonormalize(W)

    # assume dims are small: find all singular values with a dense solver.
    # (if #columns is not small, it would be better to use an iterative solver.)
    σ = svdvals(QV' * QW) # returned in descending order
    # σ's are cos(θ). We want to return ∑sin^2(θ)
    return sum(1 .- σ.^2)
end

end # module

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
using RCall
using Distributions
using StatsBase
using Iterators
using Gurobi
using JLD
using Test
using EmpiricalBayes

using KernelDensity
nreps =20
k_max = 12 #in total 25*16 replications

true_dists = [MixtureModel([ Normal(-0.3,.2), Normal(0,0.9)],[0.8, 0.2]), EmpiricalBayes.ash_flattop ]
xs = collect(-3:0.5:3)
minimax_σs = [0.2, 0.5]
ms = [20_000, 100_000]

combs = []

for (true_dist, x, σ, m) in product(true_dists, xs, minimax_σs, ms)
    push!(combs, Dict(:dist=>true_dist, :x=>x, :σ=>σ, :m=>m) )
end

#i = parse(Int64, ARGS[1])
i=1
comb = combs[i]

m = comb[:m]
σ = comb[:σ]
true_dist = comb[:dist]
x = comb[:x]

marginal_grid = collect(range(-7.0, stop=7.0, length=1001));
prior_grid = collect(range(-3, stop=3, length=121));
marginal_h = marginal_grid[2]-marginal_grid[1];

d_true = NormalConvolutionProblem(true_dist, marginal_grid)
ds = MixingNormalConvolutionProblem(Normal, σ, prior_grid, marginal_grid);

target = PosteriorTarget(LFSRNumerator(x))

Xs_train = rand(d_true,Int(m/2))

Xs_test = rand(d_true, Int(m/2))

dry_run = CEB_ci(Xs_train, Xs_test, ds, target)

m_train = length(Xs_train)
m_test = length(Xs_test)

marginal_grid = ds.marginal_grid
marginal_h = ds.marginal_h

M_bd = marginal_h/sqrt(2*pi)

f_const = BinnedMarginalDensity(M_bd, marginal_grid, marginal_h)
# Train
# Estimate the numerator
M_max_num = MinimaxCalibrator(ds, f_const, m_train, target.num;
             C=Inf);

# Estimate the denominator
M_max_denom = MinimaxCalibrator(ds, f_const, m_train, target.denom;
            C=Inf);

f_nb = fit(BinnedMarginalDensityNeighborhood, Xs_train, ds)

C_bias = f_nb.C_bias

f_nb.f_kde

CEB_ci(Xs_train, Xs_test, ds, target,
            M_max_num,
            M_max_denom,
            C_bias;
            C=:auto, conf=0.1)

tmp1 = sinc_kde(Xs_train, marginal_grid)

marginal_grid2 = collect(range(-6.0, stop=6.0, length=1001));
tmp2= sinc_kde(Xs_train, marginal_grid2)

marginal_h
tmp3 = kde(Xs_train, -7:(marginal_h-1e-10):7)

marginal_h2 = marginal_grid2[2]-marginal_grid2[1]
tmp4 = kde(Xs_train, -6:marginal_h2:6)


text:
immutable FloatRect
    left::Cfloat
    top::Cfloat
    width::Cfloat
    height::Cfloat
end

immutable IntRect
    left::Cint
    top::Cint
    width::Cint
    height::Cint
end

function contains(rect::FloatRect, x::Real, y::Real)
    rectPtr = Ref{FloatRect}(rect)
    return ccall((:sfFloatRect_contains, libcsfml_graphics), Int32, (Ref{FloatRect}, Cfloat, Cfloat,), (rectPtr), x, y) == 1
end

function contains(rect::IntRect, x::Int, y::Int)
    rectPtr = Ref{FloatRect}(rect)
    return ccall((:sfIntRect_contains, libcsfml_graphics), Int32, (Ref{IntRect}, Cint, Cint,), (rectPtr), x, y) == 1
end

function intersects(rect1::FloatRect, rect2::FloatRect)
    rect1Ptr = Ref{FloatRect}(rect1)
    rect2Ptr = Ref{FloatRect}(rect2)
    return Bool(ccall((:sfFloatRect_intersects, libcsfml_graphics), Int32, (Ref{FloatRect}, Ref{FloatRect}, Ptr{FloatRect},), (rect1Ptr), (rect2Ptr), C_NULL))
end

function intersects(rect1::IntRect, rect2::IntRect)
    rect1Ptr = Ref{FloatRect}(rect1)
    rect2Ptr = Ref{FloatRect}(rect2)
    return ccall((:sfIntRect_intersects, libcsfml_graphics), Int32, (Ref{IntRect}, Ref{IntRect}, Ptr{IntRect},), (rect1Ptr), (rect2Ptr), C_NULL) == 1
end


text:
module FractionalDelayFilters

using LinearAlgebra, FFTW

include("../src/Thirian/fractional_shift.jl")

export FiniteSignalType, resetsignal!, applyshift!, evalsignal

end # module

max_stars_repo_name: RoyCCWang/FractionalDelayFilters

max_issues_repo_name: RoyCCWang/FractionalDelayFilters

max_forks_repo_name: RoyCCWang/FractionalDelayFilters


text:
using KitBase

u = -5:0.2:5 |> collect
prim = [1.0, 0.0, 1.0]

KitBase.∂maxwellian(u[1], prim)
KitBase.∂maxwellian(u, prim)


text:
highlighter_exists(config::CthulhuConfig) =
    Sys.which(config.highlighter.exec[1]) !== nothing

__init__() = CONFIG.enable_highlighter = highlighter_exists(CONFIG)

function highlight(io, x, lexer, config::CthulhuConfig)
    _print = endswith(x, '\n') ? print : println
    config.enable_highlighter || return _print(io, x)
    if !highlighter_exists(config)
        @warn "Highlighter command $(config.highlighter.exec[1]) does not exist."
        return _print(io, x)
    end
    cmd = `$(config.highlighter) $lexer`
    open(pipeline(cmd; stdout=io, stderr=stderr), "w") do io
        _print(io, x)
    end
end

function cthulhu_llvm(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig,
                      dump_module = false)
    @static if VERSION >= v"1.5.0-DEV.393"
        dump = InteractiveUtils._dump_function_linfo_llvm(
            mi, params.world,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            dump_module,
            optimize, debuginfo ? :source : :none, Base.CodegenParams())
    else
        dump = InteractiveUtils._dump_function_linfo(
            mi, params.world, #=native=# false,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            dump_module, #=syntax=# config.asm_syntax,
            optimize, debuginfo ? :source : :none)
    end
    highlight(io, dump, "llvm", config)
end

function cthulhu_native(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    @static if VERSION >= v"1.5.0-DEV.393"
        dump = InteractiveUtils._dump_function_linfo_native(
            mi, params.world,
            #=wrapper=# false, #=syntax=# config.asm_syntax,
            debuginfo ? :source : :none)
    else
        dump = InteractiveUtils._dump_function_linfo(
            mi, params.world, #=native=# true,
            #=wrapper=# false, #=strip_ir_metadata=# true,
            #=dump_module=# false, #=syntax=# config.asm_syntax,
            optimize, debuginfo ? :source : :none)
    end
    highlight(io, dump, "asm", config)
end

function cthulhu_ast(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    ast = definition(Expr, meth)
    if ast!==nothing
        if !config.pretty_ast
            dump(io, ast; maxdepth=typemax(Int))
        else
            show(io, ast)
            # Meta.show_sexpr(io, ast)
            # Could even highlight the above as some kind-of LISP
        end
    else
        @info "Could not retrieve AST. AST display requires Revise.jl to be loaded." meth
    end
end

function cthulhu_source(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    src, line = definition(String, meth)
    highlight(io, src, "julia", config)
end

cthulhu_warntype(args...) = cthulhu_warntype(stdout, args...)
function cthulhu_warntype(io::IO, src, rettype, debuginfo)
    if VERSION < v"1.1.0-DEV.762"
    elseif VERSION < v"1.2.0-DEV.229"
        lineprinter = Base.IRShow.debuginfo[debuginfo]
    else
        debuginfo = Base.IRShow.debuginfo(debuginfo)
        lineprinter = Base.IRShow.__debuginfo[debuginfo]
    end

    lambda_io::IOContext = io
    if src.slotnames !== nothing
        slotnames = Base.sourceinfo_slotnames(src)
        lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
        VERSION >= v"1.2" && show_variables(io, src, slotnames)
    end
    print(io, "Body")
    InteractiveUtils.warntype_type_printer(io, rettype, true)
    println(io)
    if VERSION < v"1.1.0-DEV.762"
        Base.IRShow.show_ir(lambda_io, src, InteractiveUtils.warntype_type_printer)
    else
        Base.IRShow.show_ir(lambda_io, src, lineprinter(src), InteractiveUtils.warntype_type_printer)
    end
    return nothing
end

function cthulu_typed(io::IO, debuginfo_key, CI, rettype, mi, iswarn)
    println(io)
    println(io, "│ ─ $(string(Callsite(-1, MICallInfo(mi, rettype))))")

    if iswarn
        cthulhu_warntype(io, CI, rettype, debuginfo_key)
    elseif VERSION >= v"1.1.0-DEV.762"
        show(io, CI, debuginfo = debuginfo_key)
    else
        show(io, MIME"text/plain"(), CI=>rettype)
    end
    println(io)
end

function show_variables(io, src, slotnames)
    println(io, "Variables")
    slottypes = src.slottypes
    for i = 1:length(slotnames)
        print(io, "  ", slotnames[i])
        if isa(slottypes, Vector{Any})
            InteractiveUtils.warntype_type_printer(io, slottypes[i], true)
        end
        println(io)
    end
    println(io)
end

# These are standard code views that don't need any special handling,
# This namedtuple maps toggle::Symbol to function
const codeviews = (;
    llvm=cthulhu_llvm,
    native=cthulhu_native,
    ast=cthulhu_ast,
    source=cthulhu_source,
)

"""
    Cthulhu.Bookmark

A `Cthulhu.Bookmark` remembers a method marked by `b` key during a descent.
It can be used with the following functions:

* `descend(::Bookmark)`, `descend_code_typed(::Bookmark)`,
  `descend_code_warntype(::Bookmark)`: continue the descent.
* `code_typed(::Bookmark)`, `code_warntype([::IO,] ::Bookmark)`: show typed IR
* `code_llvm([::IO,] ::Bookmark)`: pretty-print LLVM IR
* `code_native([::IO,] ::Bookmark)`: pretty-print native code
"""
struct Bookmark
    mi::MethodInstance
    params::CompilerParams
end

"""
    Cthulhu.BOOKMARKS :: Vector{Bookmark}

During a descent, methods can be "bookmarked" by pressing `b` key.  It
pushes a [`Cthulhu.Bookmark`](@ref) into `Cthulhu.BOOKMARKS`.  This can be
used to, e.g., continue descending by `descend(Cthulhu.BOOKMARKS[end])`.
See [`Cthulhu.Bookmark`](@ref) for other usages.
"""
const BOOKMARKS = Bookmark[]

# Default `show` is broken for `Core.Compiler.Params`. Trying not invoke it.
Base.show(io::IO, b::Bookmark) =
    print(io, "Cthulhu.Bookmark(", b.mi, ", ::", CompilerParams, ")")

# Turn off `optimize` and `debuginfo` for default `show` so that the
# output is smaller.
function Base.show(io::IO, ::MIME"text/plain", b::Bookmark;
                   optimize = false, debuginfo = :none, iswarn=false)
    CI, rt = InteractiveUtils.code_typed(b, optimize = optimize)
    if get(io, :typeinfo, Any) === Bookmark  # a hack to check if in Vector etc.
        print(io, Callsite(-1, MICallInfo(b.mi, rt)))
        print(io, " (world: ", b.params.world, ")")
        return
    end
    println(io, "Cthulhu.Bookmark (world: ", b.params.world, ")")
    cthulu_typed(io, debuginfo, CI, rt, b.mi, iswarn)
end

function InteractiveUtils.code_typed(b::Bookmark; optimize = true)
    (CI, rt, slottypes) = do_typeinf_slottypes(b.mi, optimize, b.params)
    preprocess_ci!(CI, b.mi, optimize, CONFIG)
    return CI => rt
end

InteractiveUtils.code_warntype(b::Bookmark; kw...) =
    InteractiveUtils.code_warntype(stdout, b; kw...)
function InteractiveUtils.code_warntype(io::IO, b::Bookmark; debuginfo = :source, kw...)
    CI, rt = InteractiveUtils.code_typed(b; kw...)
    cthulhu_warntype(io, CI, rt, debuginfo)
end

InteractiveUtils.code_llvm(b::Bookmark) = InteractiveUtils.code_llvm(stdout, b)
InteractiveUtils.code_llvm(io::IO, b::Bookmark; optimize = true, debuginfo = :source,
                           dump_module = false, config = CONFIG) =
    cthulhu_llvm(io, b.mi, optimize, debuginfo == :source, b.params, config, dump_module)

InteractiveUtils.code_native(b::Bookmark; kw...) =
    InteractiveUtils.code_native(stdout, b; kw...)
InteractiveUtils.code_native(io::IO, b::Bookmark; optimize = true, debuginfo = :source,
                             config = CONFIG) =
    cthulhu_native(io, b.mi, optimize, debuginfo == :source, b.params, config)


text:
module TestTransformVariables

include("preamble.jl")
using TransformVariables

@testset "setting(asℝ₊)" begin
    l = (@lens _.y[2]) ∘ setting(asℝ₊)
    obj = (x=0, y=(0, 1, 2))
    @test get(obj, l) == 0
    @test set(obj, l, -1).y[2] == exp(-1)
end

@testset "getting(as𝕀)" begin
    l = (@lens _.y[2]) ∘ getting(as𝕀)
    obj = (x=0, y=(0, -Inf, 2))
    @test get(obj, l) == 0
    @test set(obj, l, 0.5).y[2] == 0.0
end

end  # module


text:
# following the tutorial on https://julialang.org/blog/2021/01/precompile_tutorial/
# we implement some precompile solutions

precompile(DiscretizedOperator, (Ulam, PwMap))
precompile(DiscretizedOperator, (Hat, PwMap))
precompile(powernormbounds, (Ulam, PwMap, Int64, Int64, IntegralPreservingDiscretizedOperator))
precompile(powernormbounds, (Hat, PwMap, Int64, Int64, NonIntegralPreservingDiscretizedOperator))
precompile(invariant_vector, (Ulam, IntegralPreservingDiscretizedOperator))
precompile(invariant_vector, (Hat, NonIntegralPreservingDiscretizedOperator))


text:
struct UniformContext <: ContextDist
    lims::Matrix{Float64}
    dist::Vector{Uniform}
end

function UniformContext(lims::Matrix{Float64})
    dist = [Uniform(lims[d, 1], lims[d, 2]) for d in 1:size(lims, 1)]
    return UniformContext(lims, dist)
end

function sample(T::Int64, cdist::UniformContext)
    X = Matrix{Float64}(undef, T, size(cdist.lims, 1))
    for d in 1:size(cdist.lims, 1)
        X[:, d] = rand(cdist.dist[d], T)
    end
    return X
end