text:
"""
Implementation of the leader side of the Paxos algorithm.
"""
module Leaders

export lead

using Sockets

using ..Ballots
using ..Configurations
using ..Ledgers
using ..Nodes
using ..Protocols
using ..Transports.Common
using ..Utils

mutable struct Leader
  id::NodeID
  address::String
  ledger::Ledger
  requests::Channel
  listener::Union{Listener,Nothing}
end

Leader(id::NodeID, address::String, ledger::Ledger) =
  Leader(id, address, ledger, Channel(), nothing)

function serve(leader::Leader, transport::Transport) end

function lead(leaderID::NodeID, ledger::Ledger, cfg::Configuration, transport::Transport)
  @sync begin
    cluster = Cluster(50, cfg, Connections(transport))
    leader = Leader(leaderID, memberAddress(cfg, leaderID), ledger)
    server = @async serve(leader, transport)
    try
      while isopen(leader.requests)
        clientRequests = readAvailable(leader.requests)
        clients, ballots = recordRequests(ledger, leader.id, clientRequests)
        try
          choices = prepareBallots(cluster, ledger, ballots)
          promises = proposeBallots(cluster, ledger, choices)
          acceptances = acceptBallots(cluster, ledger, promises)
          reportResults(clients, ballots, choices, acceptances)
        catch ex
          result =
            isa(ex, NoQuorumAvailableException) ? Result.noQuorum : Result.requestFailed
          reportFailures(clients, ballots, result)
        end
      end
    finally
      close(server)
    end
  end
end

function recordRequests(
  leader::Leader,
  leaderID::NodeID,
  clientRequests::Vector{Tuple{Messenger,Request}},
)::Tuple{Dict{InstanceID,Messenger},Dict{InstanceID,Ballot}}
  ledger = leader.ledger
  clients = Dict{InstanceID,Messenger}()
  ballots = Dict{InstanceID,Ballot}()
  for clientRequest in clientRequests
    client, request = clientRequest
    ballot = addEntry(ledger, leaderID, request)
    instanceID = ballot.number.instanceID
    clients[instanceID] = client
    ballots[instanceID] = ballot
  end
  clients, ballots
end

"""
Execute the prepare phase, and collect results. For each instance, find the request with the highest
sequence number, and use that request for the instance--or use the one originally assigned by the leader
if no other such request found.
"""
function prepareBallots(
  cluster::Cluster,
  leader::Leader,
  ballots::Dict{InstanceID,Ballot},
)::Dict{InstanceID,Ballot}
  ledger = leader.ledger
  ballotNumbers = map(values(ballots)) do ballot
    ballot.number
  end
  votes::Vector{Ballot} = prepare(cluster, ballotNumbers)
  # concatenate the returned ballots with the requested,
  # and choose the ballot with highest sequence number in each instance
  votesByInstance = groupBy(votes) do ballot
    ballot.number.instanceID
  end
  chosenBallots = Dict{InstanceID,Ballot}()
  for instanceID in keys(ballots)
    ballot = ballots[instanceID]
    votes = votesByInstance[instanceID]
    maxVote = maxBy(votes) do vote
      vote.number.sequenceNumber
    end
    if maxVote === nothing
      # great, we get to use our request!
      chosenBallots[instanceID] = ballot
    else
      # TODO if we see any vote with an equal or higher number, then odds
      # are good that the leader has fallen behind. Ultimately, ballots shouldn't
      # succeed, but may need a "catchup" procedure if falls too far behind
      # and can't fulfill enough client requests
      if maxVote.number.sequenceNumber < ballot.number.sequenceNumber
        # we have to use the decree from the earlier, highest vote
        chosenBallots[instanceID] = Ballot(ballot.number, maxVote.request)
      end
    end
  end
  prepare!(ledger, chosenBallots)
end

"""
Execute the propose phase, and collect results. For each instance, find the ballot number
chosen in the prepare phase, and then make sure that there was a quorum that voted for that
ballot number. Once found, build a map to associate instance IDs to the approved ballot numbers.
"""
function proposeBallots(
  cluster::Cluster,
  leader::Leader,
  choices::Dict{InstanceID,Ballot},
)::Dict{InstanceID,BallotNumber}
  ledger = leader.ledger
  promises::Vector{Promise} = propose(cluster, collect(values(choices)))
  votes = Dict{InstanceID,Vector{NodeID}}()
  for promise in promises
    memberID = promise.memberID
    for ballotNumber in promise.ballotNumbers
      instanceID = ballotNumber.instanceID
      choice = choices[instanceID].number
      # only include the vote if it matches the leaeder's choice;
      # should simplify quorum verification
      if ballotNumber == choice
        members = get!(votes, instanceID, Vector())
        push!(members, memberID )
      end
    end
  end
  approved = Dict{InstanceID,BallotNumber}()
  for instanceID in keys(votes)
    voters = map(first, votes[instanceID])
    if isquorum(cluster.configuration, Set(voters))
      approved[instanceID] = choices[instanceID].number
    end
  end
  promise!(ledger, approved)
end

function acceptBallots(cluster::Cluster, leader::Leader, promises::Vector{BallotNumber})
  ledger = leader.ledger
  acceptances = accept(cluster, promises)
  accept!(ledger, acceptances)
end

function reportResults(
  clients::Dict{InstanceID,Messenger},
  ballots::Dict{InstanceID,Ballot},
  choices::Dict{InstanceID,Ballot},
  acceptances::Vector{BallotNumber},
)
  for acceptance in acceptances
    instanceID = acceptance.instanceID
    acceptedRequestID = choices[instanceID].request.id
    requestID = ballots[instanceID].request.id
    if acceptedRequestID == requestID
      respond(clients[instanceID], requestID, instanceID, Result.requestAccepted)
    else
      respond(clients[instanceID], requestID, instanceID, Result.requestObsolete)
    end
    # let's clear out the client--we've responded
    delete!(clients, instanceID)
  end
  # any clients left...tell them they are obsolete
  for instanceID in keys(clients)
    client = clients[instanceID]
    requestID = ballots[instanceID].request.id
    respond(client, requestID, instanceID, Result.requestObsolete)
  end

end

function reportFailures(
  clients::Dict{InstanceID,Messenger},
  ballots::Dict{InstanceID,Ballot},
  result::Result,
)
  for instanceID in keys(clients)
    client = clients[instanceID]
    requestID = ballots[instanceID].request.id
    respond(client, requestID, instanceID, result)
  end
end

end


text:
# Test helper functions for infinite variable macro
@testset "Infinite Helpers" begin
    # _check_rhs
    @testset "_check_rhs" begin
        # test with reversed sides
        arg1 = :(data[i]); arg2 = :(x[i=1:2])
        @test InfiniteOpt._check_rhs(arg1, arg2) == (arg2, arg1)
        # test with normal case that shouldn't be swapped
        arg1 = :((x[i=1:2])(t)); arg2 = :(data[i])
        @test InfiniteOpt._check_rhs(arg1, arg2) == (arg1, arg2)
        # test reversed case that stays reversed because cannot be distinguished
        arg1 = :(data(i)); arg2 = :((x[i=1:2])(t))
        @test InfiniteOpt._check_rhs(arg1, arg2) == (arg1, arg2)
    end
    # _less_than_parse
    @testset "_less_than_parse" begin
        # test with reversed sides
        arg1 = :(data[i]); arg2 = :(x[i=1:2])
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:($(arg2) >= $(arg1)),
                                                           nothing)
        # test normal reference with parameter tuple
        arg1 = :(x[i=1:2](t)); arg2 = :(data[i])
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x[i=1:2] <= $(arg2)),
                                                           :((t,)))
        # test normal with parameter tuple
        arg1 = :(x(t)); arg2 = :(data[i])
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x <= $(arg2)),
                                                           :((t,)))
        # test normal without tuple
        arg1 = :(x); arg2 = :(data[i])
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x <= $(arg2)),
                                                           nothing)
        # test normal without tuple
        arg1 = :(x); arg2 = 1
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x <= $(arg2)),
                                                           nothing)
    end
    # _greater_than_parse
    @testset "_greater_than_parse" begin
        # test with reversed sides
        arg1 = :(data[i]); arg2 = :(x[i=1:2])
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:($(arg2) <= $(arg1)),
                                                              nothing)
        # test normal reference with parameter tuple
        arg1 = :(x[i=1:2](t)); arg2 = :(data[i])
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x[i=1:2] >= $(arg2)),
                                                              :((t,)))
        # test normal with parameter tuple
        arg1 = :(x(t)); arg2 = :(data[i])
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x >= $(arg2)),
                                                              :((t,)))
        # test normal without tuple
        arg1 = :(x); arg2 = :(data[i])
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x >= $(arg2)),
                                                              nothing)
        # test normal without tuple
        arg1 = :(x); arg2 = 1
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x >= $(arg2)),
                                                              nothing)
    end
    # _less_than_parse (number on lhs)
    @testset "_less_than_parse (reversed)" begin
        # test with reference
        arg1 = 1; arg2 = :(x[i=1:2])
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:($(arg2) >= $(arg1)),
                                                           nothing)
        # test with reference and parameter tuple
        arg1 = 1; arg2 = :(x[i=1:2](t))
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x[i=1:2] >= $(arg1)),
                                                           :((t,)))
        # test normal without tuple
        arg1 = 1; arg2 = :(x)
        @test InfiniteOpt._less_than_parse(arg1, arg2) == (:(x >= $(arg1)),
                                                           nothing)
    end
    # _greater_than_parse (number on lhs)
    @testset "_greater_than_parse (reversed)" begin
        # test with reference
        arg1 = 1; arg2 = :(x[i=1:2])
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:($(arg2) <= $(arg1)),
                                                              nothing)
        # test with reference and parameter tuple
        arg1 = 1; arg2 = :(x[i=1:2](t))
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x[i=1:2] <= $(arg1)),
                                                              :((t,)))
        # test normal without tuple
        arg1 = 1; arg2 = :(x)
        @test InfiniteOpt._greater_than_parse(arg1, arg2) == (:(x <= $(arg1)),
                                                              nothing)
    end
    # _equal_to_parse
    @testset "_equal_to_parse" begin
        # test with reversed sides
        arg1 = :(data[i]); arg2 = :(x[i=1:2])
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:($(arg2) == $(arg1)),
                                                          nothing)
        # test normal reference with parameter tuple
        arg1 = :(x[i=1:2](t)); arg2 = :(data[i])
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x[i=1:2] == $(arg2)),
                                                          :((t,)))
        # test normal with parameter tuple
        arg1 = :(x(t)); arg2 = :(data[i])
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x == $(arg2)),
                                                          :((t,)))
        # test normal without tuple
        arg1 = :(x); arg2 = :(data[i])
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x == $(arg2)),
                                                          nothing)
        # test normal without tuple
        arg1 = :(x); arg2 = 1
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x == $(arg2)),
                                                          nothing)
    end
    # _equal_to_parse (number on lhs)
    @testset "_equal_to_parse (reversed)" begin
        # test with reference
        arg1 = 1; arg2 = :(x[i=1:2])
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:($(arg2) == $(arg1)),
                                                          nothing)
        # test with reference and parameter tuple
        arg1 = 1; arg2 = :(x[i=1:2](t))
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x[i=1:2] == $(arg1)),
                                                          :((t,)))
        # test normal without tuple
        arg1 = 1; arg2 = :(x)
        @test InfiniteOpt._equal_to_parse(arg1, arg2) == (:(x == $(arg1)),
                                                          nothing)
    end
    # _parse_parameters (call)
    @testset "_parse_parameters (call)" begin
        # test less than parse
        expr = :(x(t, x) <= 1)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                            expr.args) == (:(x <= 1), :((t, x)))
        # test greater than parse
        expr = :(x[1:2](t) >= 1)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                         expr.args) == (:(x[1:2] >= 1), :((t,)))
        # test equal to parse
        expr = :(x(t) == d)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                            expr.args) == (:(x == d), :((t,)))
        # test only variable parse
        expr = :(x(t))
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                            expr.args) == (:(x), :((t,)))
        # test invalid use of operator
        expr = :(x(t) in 1)
        @test_throws ErrorException InfiniteOpt._parse_parameters(error,
                                                                  Val(expr.head),
                                                                  expr.args)
    end
    # _parse_parameters (comparison)
    @testset "_parse_parameters (compare)" begin
        # test with parameters
        expr = :(0 <= x(t, x) <= 1)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                       expr.args) == (:(0 <= x <= 1), :((t, x)))
        # test with parameters and references
        expr = :(0 <= x[1:2](t, x) <= 1)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                  expr.args) == (:(0 <= x[1:2] <= 1), :((t, x)))
        # test without parameters
        expr = :(0 <= x <= 1)
        @test InfiniteOpt._parse_parameters(error, Val(expr.head),
                                         expr.args) == (:(0 <= x <= 1), nothing)
    end
end

# Test the infinite variable macro
@testset "Infinite" begin
    # initialize model and infinite parameters
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= t <= 1)
    @infinite_parameter(m, -1 <= x[1:2] <= 1)
    # test single variable definition
    @testset "Single" begin
        # test basic defaults
        vref = InfiniteVariableRef(m, 1)
        @test @infinite_variable(m, parameter_refs = t) == vref
        @test name(vref) == "noname(t)"
        # test more tuple input and variable details
        vref = InfiniteVariableRef(m, 2)
        @test @infinite_variable(m, parameter_refs = (t, x), base_name = "test",
                                 binary = true) == vref
        @test name(vref) == "test(t, x)"
        @test is_binary(vref)
        # test nonanonymous with simple single arg
        vref = InfiniteVariableRef(m, 3)
        @test @infinite_variable(m, a(x)) == vref
        @test name(vref) == "a(x)"
        # test nonanonymous with complex single arg
        vref = InfiniteVariableRef(m, 4)
        @test @infinite_variable(m, 0 <= b(x) <= 1) == vref
        @test name(vref) == "b(x)"
        @test lower_bound(vref) == 0
        @test upper_bound(vref) == 1
        # test nonanonymous with reversed single arg
        vref = InfiniteVariableRef(m, 5)
        @test @infinite_variable(m, 0 <= c(t)) == vref
        @test name(vref) == "c(t)"
        @test lower_bound(vref) == 0
        # test multi-argument expr 1
        vref = InfiniteVariableRef(m, 6)
        @test @infinite_variable(m, d(t) == 0, Int, base_name = "test") == vref
        @test name(vref) == "test(t)"
        @test fix_value(vref) == 0
    end
    # test array variable definition
    @testset "Array" begin
        # test anonymous array
        vrefs = [InfiniteVariableRef(m, 7), InfiniteVariableRef(m, 8)]
        @test @infinite_variable(m, [1:2], parameter_refs = t) == vrefs
        @test name(vrefs[1]) == "noname(t)"
        # test basic param expression
        vrefs = [InfiniteVariableRef(m, 9), InfiniteVariableRef(m, 10)]
        @test @infinite_variable(m, e[1:2], parameter_refs = (t, x)) == vrefs
        @test name(vrefs[2]) == "e[2](t, x)"
        # test comparison without params
        vrefs = [InfiniteVariableRef(m, 11), InfiniteVariableRef(m, 12)]
        @test @infinite_variable(m, 0 <= f[1:2] <= 1,
                                 parameter_refs = (t, x)) == vrefs
        @test name(vrefs[2]) == "f[2](t, x)"
        @test lower_bound(vrefs[1]) == 0
        @test upper_bound(vrefs[2]) == 1
        # test comparison with call
        vrefs = [InfiniteVariableRef(m, 13), InfiniteVariableRef(m, 14)]
        @test @infinite_variable(m, 0 <= g[1:2](t) <= 1) == vrefs
        @test name(vrefs[1]) == "g[1](t)"
        @test lower_bound(vrefs[1]) == 0
        @test upper_bound(vrefs[2]) == 1
        # test fixed
        vrefs = [InfiniteVariableRef(m, 15), InfiniteVariableRef(m, 16)]
        @test @infinite_variable(m, h[i = 1:2](t) == ones(2)[i]) == vrefs
        @test name(vrefs[1]) == "h[1](t)"
        @test fix_value(vrefs[1]) == 1
        # test containers
        vrefs = [InfiniteVariableRef(m, 17), InfiniteVariableRef(m, 18)]
        vrefs = convert(JuMP.Containers.SparseAxisArray, vrefs)
        @test @infinite_variable(m, [1:2](t),
                                 container = SparseAxisArray) == vrefs
        @test name(vrefs[1]) == "noname(t)"
    end
    # test errors
    @testset "Errors" begin
        # test model assertion errors
        m2 = Model()
        @test_throws AssertionError @infinite_variable(m2, parameter_refs = t)
        @test_throws AssertionError @infinite_variable(m2, i(t))
        @test_throws AssertionError @infinite_variable(m2, i, Int)
        @test_throws AssertionError @infinite_variable(m2, i(t), Bin)
        # test double specification
        @test_macro_throws ErrorException @infinite_variable(m, i(t),
                                                             parameter_refs = x)
        # test undefined parameter error
        @test_macro_throws ErrorException @infinite_variable(m, i, Int)
        # test invalid keyword arguments
        @test_macro_throws ErrorException @infinite_variable(m, i(t),
                                                           parameter_values = 1)
        @test_macro_throws ErrorException @infinite_variable(m, i(t),
                              infinite_variable_ref = InfiniteVariableRef(m, 1))
        @test_macro_throws ErrorException @infinite_variable(m, i(t), bad = 42)
        # test name duplication
        @test_macro_throws ErrorException @infinite_variable(m, a(t), Int)
    end
end

# Test the point variable macro
@testset "Point" begin
    # initialize model, parameters, and infinite variables
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= t <= 1)
    @infinite_parameter(m, -1 <= x[1:2] <= 1)
    @infinite_variable(m, 0 <= z(t, x) <= 1, Int)
    @infinite_variable(m, z2[1:2](t) == 3)
    # test single variable definition
    @testset "Single" begin
        # test simple anon case
        vref = PointVariableRef(m, 4)
        @test @point_variable(m, infinite_variable_ref = z,
                              parameter_values = (0, [0, 0])) == vref
        @test infinite_variable_ref(vref) == z
        pt = convert(JuMP.Containers.SparseAxisArray, zeros(2))
        @test parameter_values(vref) == (0, pt)
        @test is_integer(vref)
        @test lower_bound(vref) == 0
        # test anon with changes to fixed
        vref = PointVariableRef(m, 5)
        @test @point_variable(m, infinite_variable_ref = z, lower_bound = -5,
                          parameter_values = (0, [0, 0]), binary = true) == vref
        @test infinite_variable_ref(vref) == z
        @test parameter_values(vref) == (0, pt)
        @test !is_integer(vref)
        @test is_binary(vref)
        @test lower_bound(vref) == -5
        # test regular with alias
        vref = PointVariableRef(m, 6)
        @test @point_variable(m, z(0, [0, 0]), z0, Bin) == vref
        @test infinite_variable_ref(vref) == z
        @test parameter_values(vref) == (0, pt)
        @test is_binary(vref)
        @test lower_bound(vref) == 0
        @test name(vref) == "z0"
        # test regular with semi anon
        vref = PointVariableRef(m, 7)
        @test @point_variable(m, z(0, [0, 0]), base_name = "z0",
                              binary = true) == vref
        @test infinite_variable_ref(vref) == z
        @test parameter_values(vref) == (0, pt)
        @test is_binary(vref)
        @test lower_bound(vref) == 0
        @test name(vref) == "z0"
    end
    # test array variable definition
    @testset "Array" begin
        # test anon array with one infvar
        vrefs = [PointVariableRef(m, 8), PointVariableRef(m, 9)]
        @test @point_variable(m, [1:2], infinite_variable_ref = z,
                              parameter_values = (0, [0, 0])) == vrefs
        @test infinite_variable_ref(vrefs[1]) == z
        pt = convert(JuMP.Containers.SparseAxisArray, zeros(2))
        @test parameter_values(vrefs[2]) == (0, pt)
        @test is_integer(vrefs[1])
        @test lower_bound(vrefs[2]) == 0
        # test anon array with different inf vars
        vrefs = [PointVariableRef(m, 10), PointVariableRef(m, 11)]
        @test @point_variable(m, [i = 1:2], infinite_variable_ref = z2[i],
                              parameter_values = 0) == vrefs
        @test infinite_variable_ref(vrefs[1]) == z2[1]
        @test infinite_variable_ref(vrefs[2]) == z2[2]
        @test parameter_values(vrefs[2]) == (0,)
        @test fix_value(vrefs[2]) == 3
        @test name(vrefs[1]) == "z2[1](0)"
        # test array with same infvar
        vrefs = [PointVariableRef(m, 12), PointVariableRef(m, 13)]
        @test @point_variable(m, z(0, [0, 0]), a[1:2], Bin) == vrefs
        @test infinite_variable_ref(vrefs[1]) == z
        @test parameter_values(vrefs[2]) == (0, pt)
        @test is_binary(vrefs[1])
        @test lower_bound(vrefs[2]) == 0
        @test name(vrefs[1]) == "a[1]"
        # test test array with differnt infvars
        vrefs = [PointVariableRef(m, 14), PointVariableRef(m, 15)]
        @test @point_variable(m, z2[i](0), b[i = 1:2] >= -5) == vrefs
        @test infinite_variable_ref(vrefs[1]) == z2[1]
        @test infinite_variable_ref(vrefs[2]) == z2[2]
        @test parameter_values(vrefs[2]) == (0,)
        @test lower_bound(vrefs[2]) == -5
        @test name(vrefs[1]) == "b[1]"
        # test semi anon array
        vrefs = [PointVariableRef(m, 16), PointVariableRef(m, 17)]
        @test @point_variable(m, z2[i](0), [i = 1:2], lower_bound = -5) == vrefs
        @test infinite_variable_ref(vrefs[1]) == z2[1]
        @test infinite_variable_ref(vrefs[2]) == z2[2]
        @test lower_bound(vrefs[2]) == -5
        @test name(vrefs[1]) == "z2[1](0)"
    end
    # test errors
    @testset "Errors" begin
        # test model assertion errors
        m2 = Model()
        @test_throws AssertionError @point_variable(m2, infinite_variable_ref = z,
                                                 parameter_values = (0, [0, 0]))
        @test_throws AssertionError @point_variable(m2, z(0, [0, 0]), bob, Bin)
        @test_throws AssertionError @point_variable(m2, [1:2],
                      infinite_variable_ref = z, parameter_values = (0, [0, 0]))
        # test double specification
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]), bob,
                                                      infinite_variable_ref = z)
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]), bob,
                                                 parameter_values = (0, [0, 0]))
        # test the adding expressions to infvar
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]) >= 0,
                                                          bob, Bin)
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]) == 0,
                                                          bob, Bin)
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]) in 0,
                                                          bob, Bin)
        # test other syntaxes
        @test_macro_throws ErrorException @point_variable(m,
                                               0 <= z(0, [0, 0]) <= 1, bob, Bin)
        @test_macro_throws ErrorException @point_variable(m, [1:2], Bin,
                      infinite_variable_ref = z, parameter_values = (0, [0, 0]))
        @test_macro_throws ErrorException @point_variable(m, bob,
                     infinite_variable_ref = z, parameter_values = (0, [0, 0]))
        # test redefinition catch
        @test_macro_throws ErrorException @point_variable(m, z(0, [0, 0]), z0,
                                                          Bin)
    end
end

# Test the hold variable macro if no bounds are provided
# The bounded case is tested after @BDconstraint in constraints.jl
@testset "Hold" begin
    # initialize model
    m = InfiniteModel()
    @testset "@hold_variable (Not Bounded)" begin
        # test regular
        vref = HoldVariableRef(m, 1)
        @test @hold_variable(m, x >= 1, Bin) == vref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test anan
        vref = HoldVariableRef(m, 2)
        @test @hold_variable(m, binary = true, lower_bound = 1,
                               base_name = "x") == vref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test array
        vrefs = [HoldVariableRef(m, 3), HoldVariableRef(m, 4)]
        @test @hold_variable(m, y[1:2] == 2, Int) == vrefs
        @test name(vrefs[1]) == "y[1]"
        @test fix_value(vrefs[2]) == 2
        @test is_integer(vrefs[1])
        # test errors
        @test_throws AssertionError @hold_variable(Model(), z >= 1, Bin)
        @test_macro_throws ErrorException @hold_variable(m, x >= 1, Bin)
    end
    # test _make_interval_set
    @testset "_make_interval_set" begin
        expected = :(IntervalSet(0, 1))
        @test InfiniteOpt._make_interval_set(error, :([0, 1])) == expected
        @test_throws ErrorException InfiniteOpt._make_interval_set(error, :([0]))
    end
    # test _process_parameter
    @testset "_process_parameter" begin
        @test InfiniteOpt._process_parameter(:x) == :(InfiniteOpt._make_vector(x))
    end
    # test _make_bound_pair (in)
    @testset "_make_bound_pair (in)" begin
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:in)) == :($(tp) => $(set))
    end
    # test _make_bound_pair (==)
    @testset "_make_bound_pair (==)" begin
        expr = :(t == 0)
        set = :(IntervalSet(0, 0))
        tp = :(InfiniteOpt._make_vector(t))
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:(==))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison with <=)
    @testset "_make_bound_pair (<=)" begin
        expr = :(0 <= t <= 1)
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(<=)),
                                           Val(:(<=))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison with >=)
    @testset "_make_bound_pair (>=)" begin
        expr = :(1 >= t >= 0)
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(>=)),
                                           Val(:(>=))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison fallback)
    @testset "_make_bound_pair (Fallback)" begin
        expr = :(t <= 0)
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr,
                                                                 Val(:(<=)))
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr,
                                                         Val(:(<=)), Val(:(>=)))
    end
    # test _make_bound_pair (wrapper)
    @testset "_make_bound_pair (Wrapper)" begin
        # test :call expression
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:in)) == :($(tp) => $(set))
        # test :compariosn expression
        expr = :(0 <= t <= 1)
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(<=)),
                                           Val(:(<=))) == :($(tp) => $(set))
        # test error
        expr = :([2, 3])
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr)
    end
    # test _parse_parameter_bounds (Vector)
    @testset "_parse_parameter_bounds (Vector)" begin
        args = [:(t in [0, 1]), :(0 <= t <= 1)]
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        dict = :(Dict($(dict_arg), $(dict_arg)))
        bounds = :(ParameterBounds($dict))
        @test InfiniteOpt._parse_parameter_bounds(error, args) == bounds
    end
    # test _parse_parameter_bounds (Expr)
    @testset "_parse_parameter_bounds (Expr)" begin
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        dict = :(Dict($(dict_arg)))
        bounds = :(ParameterBounds($dict))
        @test InfiniteOpt._parse_parameter_bounds(error, expr) == bounds
    end
    # test _extract_bounds (:call)
    @testset "_extract_bounds (:call)" begin
        # test single anonymous bound
        args = [:in, :t, :([0, 1])]
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(Dict($(dict_arg))))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:call)) == (nothing, bounds)
    end
    # test _extract_bounds (:tuple)
    @testset "_extract_bounds (:tuple)" begin
        args = [:(t in [0, 1]), :(0 <= t <= 1)]
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(Dict($(dict_arg), $(dict_arg))))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:tuple)) == (nothing, bounds)
    end
    # test _extract_bounds (:comparison)
    @testset "_extract_bounds (:comparison)" begin
        args = [0, :(<=), :t, :(<=), 1]
        set = :(IntervalSet(0, 1))
        tp = :(InfiniteOpt._make_vector(t))
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(Dict($(dict_arg))))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:comparison)) == (nothing, bounds)
    end
    # test @hold_variable with bounds
    @testset "@hold_variable (Bounded)" begin
        # prepare parameters
        @infinite_parameter(m, 0 <= par <= 10)
        @infinite_parameter(m, 0 <= pars[1:2] <= 10)
        # test regular call
        vref = HoldVariableRef(m, 5)
        @test @hold_variable(m, xb >= 1, Bin,
                             parameter_bounds = (par == 0)) == vref
        @test name(vref) == "xb"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        @test m.has_hold_bounds
        @test m.vars[5].parameter_bounds == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test supports(par) == [0]
        # test regular tuple
        vref = HoldVariableRef(m, 6)
        @test @hold_variable(m, yb, parameter_bounds = (pars == 0,
                                                       par in [0, 1])) == vref
        @test name(vref) == "yb"
        dict = Dict(pars[1] => IntervalSet(0, 0), pars[2] => IntervalSet(0, 0),
                    par => IntervalSet(0, 1))
        @test m.vars[6].parameter_bounds == ParameterBounds(dict)
        @test supports(pars[1]) == [0]
        # test comparison
        vref = HoldVariableRef(m, 7)
        @test @hold_variable(m, zb, parameter_bounds = (0 <= par <= 1),
                             base_name = "bob") == vref
        @test name(vref) == "bob"
        @test m.vars[7].parameter_bounds == ParameterBounds(Dict(par => IntervalSet(0, 1)))
        # test unrecognized format
        @test_macro_throws ErrorException @hold_variable(m, parameter_bounds = par)
        # test container specification
        vrefs = [HoldVariableRef(m, 8), HoldVariableRef(m, 9)]
        vrefs = convert(JuMPC.SparseAxisArray, vrefs)
        @test @hold_variable(m, [i = 1:2], parameter_bounds = (pars[i] == 0),
                             container = SparseAxisArray) == vrefs
        @test name(vrefs[1]) == ""
        @test m.vars[9].parameter_bounds == ParameterBounds(Dict(pars[2] => IntervalSet(0, 0)))
        # test wrong model type
        @test_macro_throws ErrorException @hold_variable(Model(),
                                                  parameter_bounds = (par == 0))
    end
end

ext: jl


text:
import OrderedCollections: OrderedDict
import Highlights: Lexers

using Highlights: Highlights
using Highlights.Format

# ------------------------------- highlighting ------------------------------- #
highlight_regexes = OrderedDict(
    :number => (r"(?<group>(?<![a-zA-Z0-9_#])\d+(\.\d*)?+([eE][+-]?\d*)?)",),
    :operator =>
        (r"(?<group>(?<!\{)\/)", r"(?<group>(?![\:\<])[\+\-\*\%\^\&\|\!\=\>\<\~])"),
    :string => (r"(?<group>[\"\']{3}|[\'\"]{1}(\n|.)*?[\"\']{3}|[\'\"]{1})",),
    :code => (r"(?<group>([\`]{3}|[\`]{1})(\n|.)*?([\`]{3}|[\`]{1}))",),
    :expression => (r"(?<group>\:\(+.+[\)])",),
    :symbol => (r"(?<group>(?<!\:)(?<!\:)\:\w+)",),
    :emphasis_light => (r"(?<group>[\[\]\(\)])", r"(?<group>@\w+)"),
    :type => (r"(?<group>\:\:[\w\.]*)", r"(?<group>\<\:\w+)"),
)

"""
    highlight(text::AbstractString, theme::Theme)

Highlighs a text introducing markup to style semantically
relevant segments, colors specified by a theme object
"""
function highlight(text::AbstractString; theme::Theme = term_theme[])
    has_ansi(text) && return text

    # highlight with regexes
    for (symb, rxs) in pairs(highlight_regexes)
        markup = getfield(theme, symb)
        open, close = "{$markup}", "{/$markup}"
        for rx in rxs
            text = replace(text, rx => SubstitutionString(open * s"\g<0>" * close))
        end
    end

    return remove_markup(apply_style(text))
end

"""
    highlight(text::AbstractString, theme::Theme, like::Symbol)

Hilights an entire text as if it was a type of semantically
relevant text of type :like.
"""
function highlight(text::AbstractString, like::Symbol; theme::Theme = term_theme[])
    markup = getfield(theme, like)
    return apply_style(
        do_by_line((x) -> "{" * markup * "}" * x * "{/" * markup * "}", chomp(text)),
    )
end

# shorthand to highlight objects based on type
highlight(x::Union{UnionAll,DataType}; theme::Theme = term_theme[]) =
    highlight(string(x), :type; theme = theme)

highlight(x::Number; theme::Theme = term_theme[]) =
    highlight(string(x), :number; theme = theme)

highlight(x::Function; theme::Theme = term_theme[]) =
    highlight(string(x), :func; theme = theme)

highlight(x::Symbol; theme::Theme = term_theme[]) =
    highlight(string(x), :symbol; theme = theme)

highlight(x::Expr; theme::Theme = term_theme[]) =
    highlight(string(x), :expression; theme = theme)

highlight(x::AbstractVector; theme::Theme = term_theme[]) =
    highlight(string(x), :number; theme = theme)

highlight(x; theme = term_theme[]) = string(x)  # capture all other cases

# ------------------------------ Highlighters.jl ----------------------------- #

"""
    Format.render(io::IO, ::MIME"text/ansi", tokens::Format.TokenIterator)

custom ANSI lexer for Highlighters.jl
"""
function Format.render(io::IO, ::MIME"text/ansi", tokens::Format.TokenIterator)
    for (str, id, style) in tokens
        fg = style.fg.active ? map(Int, (style.fg.r, style.fg.g, style.fg.b)) : ""
        bg = style.bg.active ? map(Int, (style.bg.r, style.bg.g, style.bg.b)) : nothing

        bold = style.bold ? "bold" : ""
        italic = style.italic ? "italic" : ""
        underline = style.underline ? "underline" : ""
        bg = isnothing(bg) ? "" : "on_$(bg)"
        markup = "$fg $(bg) $(bold) $(italic) $(underline)"
        if length(strip(markup)) > 0
            print(io, "{$markup}$str{/$markup}")
        else
            print(io, str)
        end
    end
end

"""
    highlight_syntax(code::AbstractString; style::Bool=true)

Highlight Julia code syntax in a string.
"""
function highlight_syntax(code::AbstractString; style::Bool = true)
    txt = sprint(
        Highlights.highlight,
        MIME("text/ansi"),
        escape_brackets(code),
        Lexers.JuliaLexer,
        CodeTheme;
        context = stdout,
    )
    txt = unescape_brackets(txt)
    style && (txt = apply_style(txt))

    return do_by_line(rstrip, remove_markup(txt))
end

"""
    load_code_and_highlight(path::AbstractString, lineno::Int; δ::Int=3, width::INt=120)

Load a file, get the code and format it. Return styled text
"""
function load_code_and_highlight(path::AbstractString, lineno::Int; δ::Int = 3)
    η = countlines(path)
    @assert lineno > 0 "lineno must be ≥1"
    @assert lineno ≤ η "lineno $lineno too high for file with $(η) lines"

    lines = read_file_lines(path, lineno - 9, lineno + 10)

    linenos = first.(lines)
    lines = [ln[2] for ln in lines]
    code = split(highlight_syntax(join(lines); style = false), "\n")

    # clean
    clean(line) = replace(line, "    {/    }" => "")
    codelines = clean.(code)  # [10-δ:10+δ]
    linenos = linenos  # [10-δ:10+δ]

    # make n lines match
    if lineno ≤ δ
        codelines = clean.(code)[1:(lineno + δ)]
        linenos = linenos[1:(lineno + δ)]
    elseif η - lineno ≤ δ
        codelines = clean.(code)[(end - δ):end]
        linenos = linenos[(end - δ):end]
    else
        codelines = clean.(code)[(10 - δ):(10 + δ)]
        linenos = linenos[(10 - δ):(10 + δ)]
    end

    # format
    _len = textlen ∘ lstrip
    dedent = 100
    for ln in codelines
        if _len(ln) > 1
            dedent = min(dedent, textlen(ln) - _len(ln))
        end
    end
    dedent = dedent < 1 ? 1 : dedent

    cleaned_lines = []
    for (n, line) in zip(linenos, codelines)
        # style
        symb, color = if n == lineno
            "{red bold}❯{/red bold}", "white"
        else
            " ", "grey39"
        end

        # end
        line = textlen(line) > 1 ? lpad(line[dedent:end], 8) : line
        push!(cleaned_lines, symb * " {$color}$n{/$color} " * line)
    end

    return join(cleaned_lines, "\n")
end


text:
"""
Makes a coordinate given latitude and longitude
#Example:

```julia
Coordinate(37.5, 67.91666)
```
"""
struct Coordinate
    latitude
    longitude
end
"""
The mapping between the coordinates of earth and indices of the image is needed to convert from one to another. In this package such mapping is referred to as a coordinate system. To define a coordinate system, you need to provide coordinates of the top-left and bottom-right pixels, and the height and the width of the image.

Suppose the image we plan to load had 8000 rows and 7100 columns. The top-left coordinate (37.5, 67.91666) is and bottom-right (4.166, 97.5). We can use these to create a the mapping.

```julia
top_left     = Coordinate(37.5, 67.91666)
bottom_right = Coordinate(4.166, 97.5)
height = 8000
width  = 7100
my_coordinate_system = CoordinateSystem(top_left, bottom_right, height, width)
```

Now functions can use this mapping to go from one to another.
"""
struct CoordinateSystem
    top_left::Coordinate
    bottom_right::Coordinate
    height::Int64
    width::Int64
end
"""
To convert from latitude to row number of an image, use the ```lat_to_row``` function.
# Example:
```julia
lat_to_row(my_coordinate_system, 19.6)
```
"""
function lat_to_row(geometry::CoordinateSystem, x)
    top_left        = geometry.top_left.latitude
    bottom_right    = geometry.bottom_right.latitude
    height          = geometry.height
    return Int(round(((x - top_left) * height / (bottom_right - top_left))))
end
"""
To convert from longitude to column number of an image use the ```long_to_column``` function.
# Example:
```julia
long_to_column(my_coordinate_system, 73.33)
```
"""
function long_to_column(geometry::CoordinateSystem, x)
    top_left        = geometry.top_left.longitude
    bottom_right    = geometry.bottom_right.longitude
    width           = geometry.width
    return Int(round(((x - top_left) * width/(bottom_right - top_left))))
end
"""
To convert from row number of an image to latitude use the ```row_to_lat``` function.
# Example:
```julia
row_to_lat(my_coordinate_system, 100)
```
"""
function row_to_lat(geometry::CoordinateSystem, x)
    top_left        = geometry.top_left.latitude
    bottom_right    = geometry.bottom_right.latitude
    height          = geometry.height
    return  ((bottom_right - top_left)/height * x + top_left)
end
"""
To convert from column number of an image to longitude use the ```column_to_long``` function.
# Example:
```julia
column_to_long(my_coordinate_system, 100)
```
"""
function column_to_long(geometry::CoordinateSystem, x)
    top_left        = geometry.top_left.longitude
    bottom_right    = geometry.bottom_right.longitude
    width           = geometry.width
    return((bottom_right - top_left)/width * x + top_left)
end

"""
To convert a coordinate from earth's coordinate system to image's row number and column number, use the ```coordinate_to_image``` function.
# Example:
```julia
coordinate_to_image(my_coordinate_system, Coordinate(19.49907, 72.721252))
```
"""
function coordinate_to_image(geometry::CoordinateSystem, x::Coordinate)
    return [lat_to_row(geometry, x.latitude), long_to_column(geometry, x.longitude)]
end
"""
To convert a coordinate from image's row number and column number to earth's coordinate system use the ```image_to_coordinate``` function.
# Example:
```julia
image_to_coordinate(my_coordinate_system, [4320, 1153])
```
"""
function image_to_coordinate(geometry::CoordinateSystem, x)
    return Coordinate(row_to_lat(geometry, x[1]), column_to_long(geometry, x[2]))
end

"""
Translates a coordinate system given the new top-left and bottom-right coordinates
# Example:
```julia
MUMBAI = translate_geometry(INDIA_COORDINATE_SYSTEM, Coordinate(19.49907, 72.721252), Coordinate(18.849475, 73.074187))
```
"""
function translate_geometry(geometry::CoordinateSystem, top_left::Coordinate, bottom_right::Coordinate)
    top_left_row = lat_to_row(geometry, top_left.latitude)
    bottom_right_row = lat_to_row(geometry, bottom_right.latitude)
    bottom_right_col = long_to_column(geometry, bottom_right.longitude)
    top_left_col = long_to_column(geometry, top_left.longitude)
    height = round(bottom_right_row - top_left_row)
    width = round(bottom_right_col - top_left_col)
    return CoordinateSystem(top_left, bottom_right, height, width)
end

global TILE1_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(75, -180), Coordinate(0, -60), 18000, 28800)
global TILE2_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(75, -60), Coordinate(0, 60), 18000, 28800)
global TILE3_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(75, 60), Coordinate(0, 180), 18000, 28800)
global TILE4_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(0, -180), Coordinate(-65, -60), 15600, 28800)
global TILE5_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(0, -60), Coordinate(-65, 60), 15600, 28800)
global TILE6_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(0, 60), Coordinate(-65, 180), 15600, 28800)
global FULL_COORDINATE_SYSTEM  = CoordinateSystem(Coordinate(75, -180), Coordinate(-65, 180), 33600, 86400)

global INDIA_COORDINATE_SYSTEM = CoordinateSystem(Coordinate(37.5, 67.91666), Coordinate(4.166, 97.5), 8000, 7100)
global MUMBAI_COORDINATE_SYSTEM = translate_geometry(INDIA_COORDINATE_SYSTEM, Coordinate(19.49907, 72.721252), Coordinate(18.849475, 73.074187))


text:
@testset "Errors" begin
    @testset begin
        config = DockerConfig(; base_image = "foo", platform = :foo)
        mktempdir() do build_directory
            @test_throws ArgumentError DockerSandbox._generate_dockerfile(config)
        end
    end

    @testset begin
        configs = [
            DockerConfig(;
                base_image = "julia:latest",
                read_only_maps = Dict("/" => "/foo"),
            ),
            DockerConfig(;
                base_image = "julia:latest",
                read_write_maps = Dict("/" => "/foo"),
            ),
        ]
        for config in configs
            with_container() do container
                code = """
                println("123")
                """
                @test_throws ArgumentError run(container, config, `julia -e $(code)`)
            end
        end
    end
end


text:
# Adaptation of BayesLogit C code by Nicholas Polson, James Scott, Jesse Windle, 2012-2019

const FOURPISQ = 4 * __PI^2
const HALFPISQ = 0.5 * __PI * __PI
const __TRUNC = 0.64
const __TRUNC_RECIP = 1 / __TRUNC

struct PolyaGammaDevRoye1Sampler <:
       Distributions.Sampleable{Distributions.Univariate,Distributions.Continuous}
    h::Float64
    z::Float64
end

struct PolyaGammaDevRoyeSampler <:
       Distributions.Sampleable{Distributions.Univariate,Distributions.Continuous}
    h::Float64
    z::Float64
end

## Calculate coefficient n in density of PG(1.0, 0.0), i.e. J* from Devroye.

function _a(n::Int64, x::Float64)
    K = (n + 0.5) * __PI
    y = 0.0
    if (x > __TRUNC)
        y = K * _exp_c(-K * K * x / 2)
    elseif (x > 0)
        expnt = -1.5 * (_log_c(0.5 * __PI) + _log_c(x)) + _log_c(K) - 2.0 * (n + 0.5)^2 / x
        y = _exp_c(expnt)
        #y = 1 / (__PI * x / 2)^(1.5) * K * 1 /(_exp_c((n + 0.5)^2 / x)^2)
    end
    return y
end

function _mass_texpon(z::Float64)
    fz = __PI^2 / 8 + z^2 / 2
    b = sqrt(1.0 / __TRUNC) * (__TRUNC * z - 1)
    a = -sqrt(1.0 / __TRUNC) * (__TRUNC * z + 1)

    x0 = _log_c(fz) + fz * __TRUNC
    xb = x0 - z + _log_c(Distributions.cdf(Distributions.Normal(), b))
    xa = x0 + z + _log_c(Distributions.cdf(Distributions.Normal(), a))

    qdivp = 4 / __PI * (_exp_c(xb) + _exp_c(xa))

    return 1 / (1 + qdivp)
end

# function rtinvchi2_dr(rng::Distributions.AbstractRNG, scale::Float64, trunc::Float64)
#     R = trunc / scale
#     E = tnorm_left_sp(rng, 1 / sqrt(R))
#     #E = Distributions.rand(rng, Distributions.TruncatedNormal(0.0, 1.0, 1 / sqrt(R), Inf))
#     X = scale / (E * E)
#     return X
# end

function rtigauss_dr(rng::Distributions.AbstractRNG, Z::Float64)
    #Z = abs(Z)
    t = __TRUNC
    X = t + 1.0
    if (__TRUNC_RECIP > Z)
        alpha = 0.0
        while (Random.rand(rng) > alpha)
            E1 = Random.randexp(rng)
            E2 = Random.randexp(rng)
            while (E1^2 > (2 * E2 / t))
                E1 = Random.randexp(rng)
                E2 = Random.randexp(rng)
            end
            X = t / (1 + E1 * t)^2
            alpha = _exp_c(-0.5 * Z^2 * X)
        end
    else
        mu = 1 / Z
        while (X > t)
            Y = (Random.randn(rng)^2)
            half_mu = mu / 2
            mu_Y = mu * Y
            X = mu + half_mu * mu_Y - half_mu * sqrt(4 * mu_Y + mu_Y^2)
            if (Random.rand(rng) > (mu / (mu + X)))
                X = mu^2 / X
            end
        end
    end
    return X
end

## Samples from PG(1, z)
function Distributions.rand(rng::Distributions.AbstractRNG, s::PolyaGammaDevRoye1Sampler)
    Z = abs(s.z) * 0.5
    ## PG(1,z) = 1/4 J*(1,Z/2)
    fz = (__PI^2 / 8) + (0.5 * Z^2)
    X = 0.0
    while true
        #num.trials = num.trials + 1
        if (Random.rand(rng) < _mass_texpon(Z))
            ## Truncated Exponential
            X = _TRUNC + Random.randexp(rng) / fz
        else
            ## Truncated Inverse Normal
            #X = Distributions.rand(rng, TruncatedInverseGaussian(1 / Z, 1.0, 0.0, _TRUNC))
            X = rtigauss_dr(rng, Z)
        end
        S = _a(0, X)
        Y = Random.rand(rng) * S
        n = 0
        cont = true
        while cont
            n += 1
            if (mod(n, 2) == 1)
                S = S - _a(n, X)
                if (Y <= S)
                    return 0.25 * X
                end
            else
                S = S + _a(n, X)
                if (Y > S)
                    cont = false
                end
            end
        end
        # if (Y <= S)
        #     return 0.25 * X
        # end
    end
end

## Sample from PG(h, z) using Devroye-like method.
## h is a natural number and z is a positive real.
##--------------------------------------------------------------
function Distributions.rand(rng::Distributions.AbstractRNG, s::PolyaGammaDevRoyeSampler)
    devroye1sampler = PolyaGammaDevRoye1Sampler(s.h, s.z)
    return mapreduce(x -> Distributions.rand(rng, devroye1sampler), +, 1:s.h)
end


text:
using Clang.wrap_c

context = wrap_c.init(output_file="../src/libfann.jl",
					  header_library=x->"libfann",
					  common_file="../src/libfann_common.jl",
					  clang_diagnostics=true)

context.options.wrap_structs = true

wrap_c.wrap_c_headers(context, ["../deps/usr/include/doublefann.h",
								"../deps/usr/include/fann.h",
								"../deps/usr/include/fann_activation.h",
								"../deps/usr/include/fann_cascade.h",
								"../deps/usr/include/fann_cpp.h",
								"../deps/usr/include/fann_data.h",
								"../deps/usr/include/fann_error.h",
								"../deps/usr/include/fann_internal.h",
								"../deps/usr/include/fann_io.h",
								"../deps/usr/include/fann_train.h"])


text:
"""
    Layer(transformation, data, positional::AbstractVector, named::AbstractDictionary)

Algebraic object encoding a single layer of a visualization. It is composed of a dataset,
positional and named arguments, as well as a transformation to be applied to those.
`Layer` objects can be multiplied, yielding a novel `Layer` object, or added,
yielding a [`AlgebraOfGraphics.Layers`](@ref) object.
"""
Base.@kwdef struct Layer
    transformation::Any=identity
    data::Any=nothing
    positional::Arguments=Arguments()
    named::NamedArguments=NamedArguments()
end

transformation(f) = Layer(transformation=f)

data(df) = Layer(data=columns(df))

mapping(args...; kwargs...) = Layer(positional=collect(Any, args), named=NamedArguments(kwargs))

⨟(f, g) = f === identity ? g : g === identity ? f : g ∘ f

function Base.:*(l1::Layer, l2::Layer)
    transformation = l1.transformation ⨟ l2.transformation
    data = isnothing(l2.data) ? l1.data : l2.data
    positional = vcat(l1.positional, l2.positional)
    named = merge(l1.named, l2.named)
    return Layer(; transformation, data, positional, named)
end

## Format for layer after processing

Base.@kwdef struct ProcessedLayer
    plottype::PlotFunc=Any
    primary::NamedArguments=NamedArguments()
    positional::Arguments=Arguments()
    named::NamedArguments=NamedArguments()
    labels::MixedArguments=MixedArguments()
    attributes::NamedArguments=NamedArguments()
end

function ProcessedLayer(processedlayer::ProcessedLayer; kwargs...)
    nt = (;
        processedlayer.plottype,
        processedlayer.primary,
        processedlayer.positional,
        processedlayer.named,
        processedlayer.labels,
        processedlayer.attributes
    )
    return ProcessedLayer(; merge(nt, values(kwargs))...)
end

"""
    ProcessedLayer(layer::Layer)

Convert `layer` to equivalent processed layer.
"""
function ProcessedLayer(layer::Layer)
    processedlayer = process_mappings(layer)
    grouped_entry = isnothing(layer.data) ? processedlayer : group(processedlayer)
    primary = map(vs -> map(getuniquevalue, vs), grouped_entry.primary)
    return layer.transformation(ProcessedLayer(grouped_entry; primary))
end

unnest(vs::AbstractArray, indices) = map(k -> [el[k] for el in vs], indices)

unnest_arrays(vs) = unnest(vs, keys(first(vs)))
unnest_dictionaries(vs) = unnest(vs, Indices(keys(first(vs))))

slice(v, c) = map(el -> getnewindex(el, c), v)

function slice(processedlayer::ProcessedLayer, c)
    labels = slice(processedlayer.labels, c)
    primary = slice(processedlayer.primary, c)
    positional = slice(processedlayer.positional, c)
    named = slice(processedlayer.named, c)
    return ProcessedLayer(processedlayer; labels, primary, positional, named)
end

function Base.map(f, processedlayer::ProcessedLayer)
    axs = shape(processedlayer)
    outputs = map(CartesianIndices(axs)) do c
        return f(slice(processedlayer.positional, c), slice(processedlayer.named, c))
    end
    positional, named = unnest_arrays(map(first, outputs)), unnest_dictionaries(map(last, outputs))
    return ProcessedLayer(processedlayer; positional, named)
end

## Get scales from a `ProcessedLayer`

uniquevalues(v::ArrayLike) = collect(uniquesorted(vec(v)))

to_label(label::AbstractString) = label
to_label(labels::ArrayLike) = reduce(mergelabels, labels)

function categoricalscales(processedlayer::ProcessedLayer, palettes)
    categoricals = MixedArguments()
    merge!(categoricals, processedlayer.primary)
    merge!(categoricals, Dictionary(filter(iscategoricalcontainer, processedlayer.positional)))
    return map(keys(categoricals), categoricals) do key, val
        palette = key isa Integer ? automatic : get(palettes, key, automatic)
        datavalues = key isa Integer ? mapreduce(uniquevalues, mergesorted, val) : uniquevalues(val)
        label = to_label(get(processedlayer.labels, key, ""))
        return CategoricalScale(datavalues, palette, label)
    end
end

function has_zcolor(pl::ProcessedLayer)
    for field in (:primary, :named, :attributes)
        haskey(getproperty(pl, field), :color) && return false
    end
    return pl.plottype <: Union{Heatmap, Contour, Contourf, Surface}
end

# This method works on a "sliced" `ProcessedLayer`
function continuousscales(processedlayer::ProcessedLayer)
    continuous = MixedArguments()
    merge!(continuous, filter(iscontinuous, processedlayer.named))
    merge!(continuous, Dictionary(filter(iscontinuous, processedlayer.positional)))

    continuousscales = map(keys(continuous), continuous) do key, val
        extrema = Makie.extrema_nan(val)
        label = to_label(get(processedlayer.labels, key, ""))
        return ContinuousScale(extrema, label)
    end
    # TODO: also encode colormap here
    if has_zcolor(processedlayer) && !haskey(continuousscales, :color)
        colorscale = get(continuousscales, 3, nothing)
        isnothing(colorscale) || insert!(continuousscales, :color, colorscale)
    end
    return continuousscales
end

## Machinery to convert a `ProcessedLayer` to a grid of slices of `ProcessedLayer`s

function compute_grid_positions(categoricalscales, primary=NamedArguments())
    return map((:row, :col), (first, last)) do sym, f
        scale = get(categoricalscales, sym, nothing)
        lscale = get(categoricalscales, :layout, nothing)
        return if !isnothing(scale)
            rg = Base.OneTo(maximum(plotvalues(scale)))
            haskey(primary, sym) ? fill(primary[sym]) : rg
        elseif !isnothing(lscale)
            rg = Base.OneTo(maximum(f, plotvalues(lscale)))
            haskey(primary, :layout) ? fill(f(primary[:layout])) : rg
        else
            Base.OneTo(1)
        end
    end
end

function rescale(p::ProcessedLayer, categoricalscales::MixedArguments)
    primary = map(keys(p.primary), p.primary) do key, values
        scale = get(categoricalscales, key, nothing)
        return rescale(values, scale)
    end
    positional = map(keys(p.positional), p.positional) do key, values
        scale = get(categoricalscales, key, nothing)
        return rescale.(values, Ref(scale))
    end

    # compute dodging information
    dodge = get(categoricalscales, :dodge, nothing)
    attributes = if isa(dodge, CategoricalScale)
        set(p.attributes, :n_dodge => maximum(plotvalues(dodge)))
    else
        p.attributes
    end

    return ProcessedLayer(p; primary, positional, attributes)
end

# Determine whether entries from a `ProcessedLayer` should be merged
function mergeable(processedlayer::ProcessedLayer)
    plottype, primary = processedlayer.plottype, processedlayer.primary
    # merge violins for correct renormalization
    plottype <: Violin && return true
    # merge stacked barplots
    plottype <: BarPlot && haskey(primary, :stack) && return true
    # do not merge by default
    return false
end

# This method works on a list of "sliced" `ProcessedLayer`s
function concatenate(pls::AbstractVector{ProcessedLayer})
    pl = first(pls)
    ns = [mapreduce(length, assert_equal, Iterators.flatten([pl.positional, pl.named])) for pl in pls]

    primary = map(key -> reduce(vcat, [fill(pl.primary[key], n) for (pl, n) in zip(pls, ns)]), keys(pl.primary))
    positional = map(key -> reduce(vcat, [pl.positional[key] for pl in pls]), keys(pl.positional))
    named = map(key -> reduce(vcat, [pl.named[key] for pl in pls]), keys(pl.named))

    return ProcessedLayer(pl; primary, positional, named)
end

function append_processedlayers!(pls_grid, processedlayer::ProcessedLayer, categoricalscales::MixedArguments)
    processedlayer = rescale(processedlayer, categoricalscales)
    tmp_pls_grid = map(_ -> ProcessedLayer[], pls_grid)
    for c in CartesianIndices(shape(processedlayer))
        pl = slice(processedlayer, c)
        rows, cols = compute_grid_positions(categoricalscales, pl.primary)
        for i in rows, j in cols
            push!(tmp_pls_grid[i, j], pl)
        end
    end

    ismergeable = mergeable(processedlayer)
    for (pls, tmp_pls) in zip(pls_grid, tmp_pls_grid)
        isempty(tmp_pls) && continue
        if ismergeable
            push!(pls, concatenate(tmp_pls))
        else
            append!(pls, tmp_pls)
        end
    end
    return pls_grid
end


text:
# get the transcript by Ensembl transcript id
function gencode_transcript(index, tid)
    for (chr, genelist) in index
        for gene in genelist
            for transcript in gene.transcripts
                if transcript.id == tid
                    return transcript
                end
            end
        end
    end
    return nothing
end

# get gene details by gene name, return a gene list
# for most cases, it will return a array with one record
function gencode_genes(index, genename)
    genes = []
    for (chr, genelist) in index
        for gene in genelist
            if gene.name == genename
                push!(genes, gene)
            end
        end
    end
    return genes
end

# find which gene and which exon chr::pos locates in
function gencode_locate(index, chr, pos)
    if !haskey(index, chr)
        return false
    end
    # largest gene found in nature is about 2.4M bp
    const LARGEST_GENE_SIZE = 2_600_000
    # we search all genes with its start_pos in [pos - LARGEST_GENE_SIZE, pos]
    genes = index[chr]
    len = length(genes)

    # binary search for pos - LARGEST_GENE_SIZE
    left = 1
    right = len
    while left < right
        cur = div(left+right, 2)
        gene = genes[cur]
        if gene.start_pos == pos - LARGEST_GENE_SIZE
            left = cur
            break
        elseif gene.start_pos > pos - LARGEST_GENE_SIZE
            right = cur - 1
        else
            left = cur + 1
        end
    end
    search_left = left

    # binary search for pos
    right = len
    while left < right
        cur = div(left+right, 2)
        gene = genes[cur]
        if gene.start_pos == pos
            right = cur
            break
        elseif gene.start_pos > pos
            right = cur - 1
        else
            left = cur + 1
        end
    end
    search_right = right

    matches = []
    for s in search_left:search_right
        gene = genes[s]
        if gene.start_pos <= pos && gene.end_pos >= pos
            result = search_in_gene(gene, pos)
            if result != false
                push!(matches, result)
            end
        end
    end
    return matches
end

function search_in_gene(gene, pos)
    for t in gene.transcripts
        if t.start_pos>pos || t.end_pos < pos
            continue
        end
        if !haskey(t.attributes, "tag") || !contains(t.attributes["tag"], "basic")
            continue
        end
        for exon in t.exons
            # strand is -, then exons are aligned big coord -> small coord
            if gene.strand == "-"
                if exon.start_pos<=pos && exon.end_pos>=pos
                    return Dict("gene"=>gene.name, "transcript"=>t.id, "type"=>"exon", "number"=>exon.number)
                elseif exon.end_pos<pos
                    return Dict("gene"=>gene.name, "transcript"=>t.id, "type"=>"intron", "number"=>exon.number - 1)
                end
            else
                if exon.start_pos<=pos && exon.end_pos>=pos
                    return Dict("gene"=>gene.name, "transcript"=>t.id, "type"=>"exon", "number"=>exon.number)
                elseif exon.start_pos>pos
                    return Dict("gene"=>gene.name, "transcript"=>t.id, "type"=>"intron", "number"=>exon.number - 1)
                end
            end
        end
    end
    return false
end


text:
using ReservoirManagement
using FactCheck
using SIUnits
using SIUnits.ShortUnits
using TimeSeries
using JuMP

include("_utils.jl")

include("reservoir_ds.jl")
include("reservoir_river.jl")
include("reservoir_reservoir.jl")
include("utils.jl")
include("scenariosgeneration.jl")
include("stats.jl")
include("solvers.jl")

FactCheck.exitstatus()


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/singular-poisson/demo_singular-poisson.py.html

module singular_poisson

using FenicsPy

@assert has_linear_algebra_backend("PETSc") "DOLFIN has not been configured with PETSc. "

parameters["linear_algebra_backend"] = "PETSc"

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)

a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Assemble system
A = assemble(a)
b = assemble(L)

# Solution Function
u = FeFunction(V)

# Create Krylov solver
solver = PETScKrylovSolver("cg")
solver.set_operator(A)

# Create vector that spans the null space and normalize
null_vec = u.vector()
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")

# Create null space basis object and attach to PETSc matrix
null_space = VectorSpaceBasis([null_vec])
as_backend_type(A).set_nullspace(null_space)

null_space.orthogonalize(b)

solver.solve(u.vector(), b)

# Output file
File("singular_poisson/output.pvd") << u

#plot(u)

end  # module singular_poisson


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule squashfs_tools_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("squashfs_tools")
JLLWrappers.@generate_main_file("squashfs_tools", UUID("eed32e3e-a7c5-5bf9-9121-5cf3ab653887"))
end  # module squashfs_tools_jll


text:
@everywhere function montecarlo!(L, T, conf, rng; nconfs=10_000)
    # set parameters & initialize
    β = 1.0 / T

    for j ∈ 1:L
        for k ∈ 1:L
                # Periodic boundary condition
            ip_1 = j + 1 > L ? j + 1 - L : j + 1
            im_1 = j - 1 < 1 ? j + L - 1 : j - 1
            jp_1 = k + 1 > L ? k + 1 - L : k + 1
            jm_1 = k - 1 < 1 ? k + L - 1 : k - 1

                # Change in energy
            spin_value = conf[j, k]
            spin_neighbor_sum = conf[im_1, k] + conf[ip_1, k] + conf[j, jp_1]
            spin_neighbor_sum += conf[j, jm_1]

            ΔE = 2.0 * spin_value * spin_neighbor_sum
                # Metropolis criteria
            if (ΔE <= 0) || (rand(rng) < exp(-β * ΔE))
                conf[j, k] *= -1 # flip sign
            end
        end
    end
end

@everywhere function ising_threaded(nconfs, L, sweeps)
    Ts = 1.2:0.1:3.4
    c = Matrix{Int32}(undef, L, L)

    # @threads for t in Ts
    rng = Xoroshiro128Plus()
    rand!(rng, c, [-1, 1])

    println("T = $t")
    @showprogress "Equilibrating..." for _ in 1:sweeps
        montecarlo!(L, t, c, rng; nconfs=nconfs)
    end

    println("Sampling...")
    flush(stdout)

    j = 0
    @inbounds for i in 1:sweeps
        montecarlo!(L, t, c, rng; nconfs=nconfs)
        if iszero(mod(i, nconfs))
            j += 1
            fileising = @sprintf "ising_%.2f_%d.jld2" t j
            @save joinpath("data", fileising) copy(c)
        end
    end
    println("Done.")
    flush(stdout)
    # end
end

@everywhere function ising_simple(nconfs, L, sweeps, t)
    # t = 1.2
    c = Matrix{Int32}(undef, L, L)

    # @threads for t in Ts
    rng = Xoroshiro128Plus()
    rand!(rng, c, [-1, 1])

    println("T = $t")
    @showprogress "Equilibrating..." for _ in 1:sweeps
        montecarlo!(L, t, c, rng; nconfs=nconfs)
    end

    println("Sampling...")
    flush(stdout)

    j = 0
    @inbounds for i in 1:sweeps
        montecarlo!(L, t, c, rng; nconfs=nconfs)
        if iszero(mod(i, nconfs))
            j += 1
            fileising = @sprintf "ising_%.2f_%d.jld2" t j
            @save joinpath("data", fileising) copy(c)
        end
    end
    println("Done.")
    flush(stdout)
end

# ising_threaded()
# ising_simple()

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
using Flux
using DiffEqSensitivity
using DifferentialEquations
using LinearAlgebra
using Statistics
using Distributions
using Random

include("./AssimilationProblem.jl")

function FourDVarSolve(prob::AssimilationProblem, nsteps::Int, opt)

    nobs = size(prob.w, 2)
    tobs_end = (nobs-1)*prob.Δt_obs
    Rinv = collect(inv(prob.R))

    function loss()
        res = Array(solve(ODEProblem(prob.f, prob.u₀, prob.tspan, prob.p), saveat=prob.Δt_obs:prob.Δt_obs:tobs_end))
        # ℓ = 0.5*(u0a-u0b)'*Binv*(u0a-u0b) + 0.5*sum((w[:,i]-h(res[:,i]))'*Rinv*(w[:,i]-h(res[:,i])) for i ∈ 1:size(res,2))

        ℓ = 0.5*sum((prob.w[:,i]-prob.h(res[:,i]))'*Rinv*(prob.w[:,i]-prob.h(res[:,i])) for i ∈ 1:size(res,2))
        return ℓ
    end

    data = Iterators.repeated((), nsteps)

    cb = function ()
        display(loss())
    end

    u0a = prob.u₀

    Flux.train!(loss, Flux.params(u0a), data, opt, cb=cb)

    return u0a
end


text:
"""
Cache for a single parameter/variable/dual.
Stores arrays chronologically by simulation timestamp.
"""
mutable struct OptimizationOutputCache
    key::OptimizationResultCacheKey
    "Contains both clean and dirty entries. Any key in data that is earlier than the first
    dirty timestamp must be clean."
    data::OrderedDict{Dates.DateTime, Array}
    "Oldest entry is first"
    dirty_timestamps::Deque{Dates.DateTime}
    stats::CacheStats
    size_per_entry::Int
    flush_rule::CacheFlushRule
end

function OptimizationOutputCache(key, flush_rule)
    return OptimizationOutputCache(
        key,
        OrderedDict{Dates.DateTime, Array}(),
        Deque{Dates.DateTime}(),
        CacheStats(),
        0,
        flush_rule,
    )
end

Base.length(x::OptimizationOutputCache) = length(x.data)
get_cache_hit_percentage(x::OptimizationOutputCache) = get_cache_hit_percentage(x.stats)
get_size(x::OptimizationOutputCache) = length(x) * x.size_per_entry
has_clean(x::OptimizationOutputCache) =
    !isempty(x.data) && !is_dirty(x, first(keys(x.data)))
has_dirty(x::OptimizationOutputCache) = !isempty(x.dirty_timestamps)
should_keep_in_cache(x::OptimizationOutputCache) = x.flush_rule.keep_in_cache

function get_dirty_size(cache::OptimizationOutputCache)
    return length(cache.dirty_timestamps) * cache.size_per_entry
end

function is_dirty(cache::OptimizationOutputCache, timestamp)
    isempty(cache.dirty_timestamps) && return false
    return timestamp >= first(cache.dirty_timestamps)
end

function Base.empty!(cache::OptimizationOutputCache)
    @assert isempty(cache.dirty_timestamps) "dirty cache was still present $(cache.key) $(cache.dirty_timestamps)"
    empty!(cache.data)
    cache.size_per_entry = 0
    return
end

"""
Add result to the cache.
"""
function add_result!(cache::OptimizationOutputCache, timestamp, array, system_cache_is_full)
    if cache.size_per_entry == 0
        cache.size_per_entry = length(array) * sizeof(first(array))
    end

    @debug "add_result!" cache.key timestamp get_size(cache)
    if haskey(cache.data, timestamp)
        throw(IS.InvalidValue("$timestamp is already stored in $(cache.key)"))
    end

    # Note that we buffer all writes in cache until we reach the flush size.
    # The entries using "should_keep_in_cache" can grow quite large for read caching.
    if system_cache_is_full && should_keep_in_cache(cache)
        if has_clean(cache)
            popfirst!(cache.data)
            @debug "replaced cache entry" LOG_GROUP_SIMULATION_STORE cache.key length(
                cache.data,
            )
        end
    end

    _add_result!(cache, timestamp, array)
    return cache.size_per_entry
end

function _add_result!(cache::OptimizationOutputCache, timestamp, data)
    cache.data[timestamp] = data
    push!(cache.dirty_timestamps, timestamp)
    return
end

function discard_results!(cache::OptimizationOutputCache, timestamps)
    for timestamp in timestamps
        pop!(cache.data, timestamp)
    end

    @debug "Removed $(first(timestamps)) - $(last(timestamps)) from cache" cache.key
    return
end

"""
Return all dirty data from the cache. Mark the timestamps as clean.
"""
function get_dirty_data_to_flush!(cache::OptimizationOutputCache)
    timestamps = [x for x in cache.dirty_timestamps]
    empty!(cache.dirty_timestamps)
    # Uncomment for performance testing of CacheFlush
    #TimerOutputs.@timeit RUN_SIMULATION_TIMER "Concatenate arrays for flush" begin
    arrays = (cache.data[x] for x in timestamps)
    arrays = cat(arrays..., dims=ndims(first(arrays)) + 1)
    #end

    return timestamps, arrays
end

function has_timestamp(cache::OptimizationOutputCache, timestamp)
    present = haskey(cache.data, timestamp)
    if present
        cache.stats.hits += 1
    else
        cache.stats.misses += 1
    end

    return present
end


text:
# A TweenManager updates all your tweens and timelines at once.
# Its main interest is that it handles the tween/timeline life-cycles for you,
# as well as the pooling constraints.
#
# Just give it a bunch of tweens or timelines and call [update] periodically,
# you don't need to care for anything else! Relax and enjoy your animations.

mutable struct TweenManager
    paused::Bool

    targets::Array{AbstractTween,1}

    function TweenManager(paused::Bool = true)
        o = new()
        o.paused = paused
        o.targets = Array{AbstractTween,1}()
        o
    end
end

# Adds a tween or timeline to the manager and starts or restarts it.
function add!(man::TweenManager, tween::AbstractTween)
    idx = findfirst(isequal(tween), man.targets)
    if isnothing(idx)
        println("TweenManager: adding target tween")
        push!(man.targets, tween)
    end

    start!(tween)
end

function contains(man::TweenManager, tween::AbstractTween)
    idx = findfirst(isequal(tween), man.targets)
    idx
end

function remove!(man::TweenManager, tween::AbstractTween)
    idx = contains(man, tween)
    if !isnothing(idx)
        println("Removing tween target")
        deleteat!(man.targets, idx)
    else
        println("Unable to find and remove tween target");
    end
end

function kill_all!(man::TweenManager)
    for target in man.targets
        kill!(target)
    end
end

function kill_target!(man::TweenManager, tween::AbstractTween)
    idx = findfirst(isequal(tween), man.targets)
    if !isnothing(idx)
        println("Killing idx:(", idx, ") tween target")
        kill!(man.targets, idx)
    else
        println("Unable to find and kill tween target")
    end
end

function pause!(man::TweenManager)
    man.paused = true
end

function resume!(man::TweenManager)
    man.paused = false
end

function start!(man::TweenManager, tween::AbstractTween)
    add!(man, tween)
end

function update!(man::TweenManager, delta::Float64)
    items_to_remove = Array{AbstractTween,1}()

    for target in man.targets
        if target.base.finished && target.base.auto_remove
            println("TweenManager: finished and auto removing")
            free!(target)
            push!(items_to_remove, target)
        end
    end

    for item in items_to_remove
        idx = findfirst(isequal(item), man.targets)
        deleteat!(man.targets, idx)
    end

    if !man.paused
        println("TweenManager: updating: ", delta)
        if delta >= 0.0
            for target in man.targets
                update!(target, delta)
            end
        else
            for i in length(man.targets):1
                update!(man.targets[i], delta)
            end
        end
    end
end

function count(man::TweenManager)
    length(man.targets)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function neighboringelements(t2t, elem)

t2te = t2t[:,elem];
nbelem = sort(unique(t2te[:]));
if nbelem[1] == 0
    nbelem = nbelem[2:end];
end

return nbelem

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using Base.Test, MuKanren, FactCheck
import Base.string
import MuKanren: conj, take, bind, mplus

facts("objects") do
  context("list functions") do
    @fact list(1, 2, 3) --> is_cons
  end
  context("cons construction") do
    p1 = list("a")
    @fact car(p1) --> "a"
    @fact cdr(p1) --> nil()
    p2 = list("a", "b")
    @fact car(p2) --> "a"
    @fact cdr(p2) --> list("b")
    p3 = list("a", "b")
    @fact car(p2) --> "a"
    @fact cdr(p2) --> list("b")
  end
end

facts("Goal construction") do
  context("ext-s") do
    @fact ext_s("a", 1, nil()) --> list(Pair("a", 1))
    @fact ext_s("a", 1, list(Pair("b", 2))) --> list(Pair("a", 1), Pair("b", 2))
  end

  context("assp") do
    res = assp(x->vars_equal(x, Var(1)), list(Pair(Var(1), "b"), Pair(Var(0), "a")))
    @fact res --> Pair(Var(1), "b")
    res = assp(x ->vars_equal(x, Var(3)), list(Pair(Var(4), 6), Pair(Var(3), 5), Pair(Var(2), 6), Pair(Var(1), 5)) )
    @fact res --> Pair(Var(3), 5)
  end

  empty_state=Pair(nil(), 0)
  call_empty_state(g) = g(empty_state)
  some_state = Pair(list(Pair(Var(0), 6)), 1)
  some_state2 = Pair(list(Pair(Var(0), 5)), 1)
  five = x -> equals(x, 5)
  fives = x -> disj(equals(x, 5), s_c -> () -> fives(x)(s_c))
  a_and_b = conj(
    call_fresh(a-> equals(a, 7)),
    call_fresh(b-> disj(
        equals(b, 5),
        equals(b, 6)
      ))
    )

  context("call/fresh") do
    @fact call_fresh(five)(empty_state) --> list(Pair(list(Pair(Var(0), 5)), 1))
    @fact call_fresh(five)(some_state) --> list(Pair(list(Pair(Var(1), 5), Pair(Var(0), 6)), 2))
    @fact call_fresh(five)(some_state2) --> list(Pair(list(Pair(Var(1), 5), Pair(Var(0), 5)), 2))
    #println(take(1, call_fresh(q-> fives(q))(empty_state)))
    @fact take(1, call_fresh(q-> fives(q))(empty_state)) --> [Pair(list(Pair(Var(0), 5)), 1)]
    @fact cdr(take(2, call_fresh(x -> fives(x))(empty_state))) --> [Pair(list(Pair(Var(0), 5)), 1)]
    #fives = x -> disj(equals(x, 5), s_c ->() -> ((fives(x))(s_c)))
    #println((call_fresh(fives))(empty_state))
    ccc = call_fresh(a-> equals(a, :([1, 2])))
    println(ccc(empty_state))
  end

  context("a-and-b") do
    a_and_b_result = a_and_b(empty_state)
    @fact car(a_and_b_result) --> Pair(list(Pair(Var(1), 5), Pair(Var(0), 7)), 2)
    @fact take(1, a_and_b_result) --> [Pair(list(Pair(Var(1), 5), Pair(Var(0), 7)), 2)]
    @fact car(cdr(a_and_b_result)) --> Pair(list(Pair(Var(1), 6), Pair(Var(0), 7)), 2)
    @fact cdr(cdr(a_and_b_result)) --> nil()

    @fact take(2, a_and_b(empty_state)) --> [
      Pair(list(Pair(Var(1), 5), Pair(Var(0), 7)), 2),
      Pair(list(Pair(Var(1), 6), Pair(Var(0), 7)), 2)]

    @fact take_all(a_and_b(empty_state)) --> [
      Pair(list(Pair(Var(1), 5), Pair(Var(0), 7)), 2),
      Pair(list(Pair(Var(1), 6), Pair(Var(0), 7)), 2)]
  end

  context("bind and mplus") do
    @fact bind(nil(), 1) --> mzero
    #@pending miniKanren.MicroKanren.mplus()
  end

  context("conj and disj") do
    exp1 = call_fresh(b->call_fresh(a->conj(equals(a, 3), equals(b, 4))))
    @fact exp1(empty_state) --> list(Pair(list(Pair(Var(0), 4), Pair(Var(1), 3)), 2))
    expression = conj(call_fresh(a->equals(a, 7)), call_fresh(b->disj(equals(b, 5), equals(b, 6))))
    @fact expression(empty_state) --> list(
      Pair(list(Pair(Var(1), 5), Pair(Var(0), 7)), 2),
      Pair(list(Pair(Var(1), 6), Pair(Var(0), 7)), 2))
    exp3 = call_fresh(a->disj(equals(a, 3), equals(a, 4)))
    println(exp3(empty_state))
    @fact exp3(empty_state) --> list(Pair(list(Pair(Var(0), 3)), 1), Pair(list(Pair(Var(0), 4)), 1))#TODO: fix
  end

end

facts("Function calls") do
  empty_state = Pair(nil(), 0)
  appendo = (l, s, out) -> disj(
    conj(equals(nil(), l), equals(s, out)),
    call_fresh(a-> call_fresh(d -> conj(
          equals(cons(a, d), l),
          call_fresh(res -> conj(
            equals(cons(a, res), out),
            s_c -> () -> appendo(d, s, res)(s_c)
            )
          )
        )
      )
    )
  )

  appendo2 = (l, s, out) -> disj(
      conj(equals(nil(), l), equals(s, out)),
      call_fresh(a ->
        call_fresh(d ->
          conj(
            equals(cons(a, d), l),
            call_fresh(res-> conj(
              s_c -> () -> appendo2(d, s, res)(s_c),
              equals(cons(a,res), out)
              ))
          )
        )
      )
  )

  ground_appendo = appendo(list("a"), list("b"), list("a", "b"))
  ground_appendo2 = appendo2(list("a"), list("b"), list("a", "b"))

  call_appendo = begin
    call_fresh(q->
      call_fresh(l ->
        call_fresh(s ->
          call_fresh(out -> conj(
              appendo(l, s, out),
              equals(cons(l, cons(s, cons(out, nil()))), q))
          )
        )
      )
    )
  end

  call_appendo2 = begin
    call_fresh(q ->
      call_fresh(l ->
        call_fresh(s ->
          call_fresh(out -> conj(
              appendo2(l, s, out),
              equals(cons(l, cons(s, cons(out, nil()))), q))
          )
        )
      )
    )
  end

  context("ground_appendo") do
    @fact car(ground_appendo(empty_state)()) --> Pair(list(list(Var(2), "b"), Pair(Var(1), nil()), Pair(Var(0), "a")), 3)
  end

  context("ground_appendo2") do
    @fact car(ground_appendo2(empty_state)()) -->
      Pair(list(list(Var(2), "b"), list(Var(1)), Pair(Var(0), "a")), 3)
  end

  context("appendo") do
    @fact take(2, call_appendo(empty_state)) --> [
        Pair(list(list(Var(0), Var(1), Var(2), Var(3)), Pair(Var(2), Var(3)), list(Var(1))), 4),
        Pair(list(list(Var(0), Var(1), Var(2), Var(3)), Pair(Var(2), Var(6)), list(Var(5)),
          cons(Var(3), cons(Var(4), Var(6))), cons(Var(1), cons(Var(4), Var(5)))), 7)
    ]
  end

  context("appendo2") do
    @fact take(2, call_appendo2(empty_state)) --> [
      Pair(list(list(Var(0), Var(1), Var(2), Var(3)), Pair(Var(2), Var(3)), list(Var(1))), 4),
      Pair(list(list(Var(0), Var(1), Var(2), Var(3)), cons(Var(3), cons(Var(4), Var(6))),
        Pair(Var(2), Var(6)), list(Var(5)), cons(Var(1), cons(Var(4), Var(5)))), 7)
    ]
  end

  context("reify-1st across appendo") do
    res = take(2, call_appendo(empty_state))
  println(res)
    @fact string(map(reify_state_1st_var, res)) --> string([
      list(nil(), "_.0", "_.0"),
      list(list("_.0"), "_.1", cons("_.0", "_.1"))])
    @fact string(map(reify_state_1st_var, res)) --> string([
      list(nil(), "_.0", "_.0"),
      list(list("_.0"), "_.1", cons("_.0", "_.1"))])
    @fact string(mk_reify(res)) --> string([
      list(nil(), "_.0", "_.0"),
      list(list("_.0"), "_.1", cons("_.0", "_.1"))])
  end

  context("many_non_ans") do
    relo = x -> begin
      call_fresh(x1 ->
        call_fresh(x2 ->
          conj(
            equals(x, cons(x1, x2)),
            disj(
              equals(x1, x2),
              s_c -> () -> relo(x)(s_c)
            )
          )
        )
      )
    end

    many_non_ans = call_fresh(x -> disj(relo(cons(5, 6)), equals(x, 3)))
    @fact take(1, many_non_ans(empty_state)) --> [Pair(list(Pair(Var(0), 3)), 1)]
  end

end


text:
using ITensors, Test
using Combinatorics: permutations

i = Index(2,"i")
j = Index(2,"j")
k = Index(2,"k")
l = Index(2,"l")

A = randomITensor(i, j, k, l)

@testset "Two index combiner" begin
    for inds_ij ∈ permutations([i,j])
        C = combiner(inds_ij...)
        B = A*C
        @test hasinds(B, l, k)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_il ∈ permutations([i,l])
        C = combiner(inds_il...)
        B = A*C
        @test hasinds(B, j, k)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_ik ∈ permutations([i,k])
        C = combiner(inds_ik...)
        B = A*C
        @test hasinds(B, j, l)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_jk ∈ permutations([j,k])
        C = combiner(inds_jk...)
        B = A*C
        @test hasinds(B, i, l)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_jl ∈ permutations([j,l])
        C = combiner(inds_jl...)
        B = A*C
        @test hasinds(B, i, k)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_kl ∈ permutations([k,l])
        C = combiner(inds_kl...)
        B = A*C
        @test hasinds(B, i, j)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
end

@testset "Three index combiner" begin
    for inds_ijl ∈ permutations([i,j,l])
        C = combiner(inds_ijl...)
        B = A*C
        @test hasindex(B, k)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_ijk ∈ permutations([i,j,k])
        C = combiner(inds_ijk...)
        B = A*C
        @test hasindex(B, l)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
    for inds_jkl ∈ permutations([j,k,l])
        C = combiner(inds_jkl...)
        B = A*C
        @test hasindex(B, i)
        @test combinedindex(C) == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test Array(permute(D, i, j, k, l)) == Array(A)
    end
end


text:
import Base: maximum

const tiny = 1e-12
const huge = 1e+12

function softmax(v)
    ret = exp.(clamp.(v .- maximum(v), -100.0, 0.0))
    ret ./ sum(ret)
end


text:
using PermPlain
using Test
using PermPlain: permcycles, cycstoperm
using Random

Random.seed!(1234)

makec() = PCycles([[1, 6, 14, 7, 3, 11], [2, 8, 12, 9, 10, 5, 15, 4]], 15)
makep() = [3, 7, 2, 4, 8, 10, 1, 6, 9, 5]

@testset "isperm, cyclelengths, numcycles" begin
    ntrials = 50
    for n in (3, 4, 5, 10, 15, 30)
        for trial in 1:ntrials
            p = randperm(n)
            sp = permsparse(p)
            @test permlist(permcycles(permmatrix(sp))) == p
            cl = cyclelengths(p)
            ncl = numcycles(p)
            for prep in (p, permcycles(p), permmatrix(p), sp)
                @test cl == cyclelengths(prep)
                @test ncl == numcycles(prep)
                @test isperm(prep)
            end
            @test permsgn(permcompose(p, p)) == 1
            @test permcommute(p, p)
        end
    end
end

@testset "permpower canoncycles" begin
    np = 3
    c = makec()
    p = permlist(c)
    pp = permpower(p, np);
    pp1 = permpower2(p, np);
    c1 = permpower(c,np);
    p2 = permlist(c1);
    @test pp == p2
    @test pp1 == p2
    @test permcycles(pp) == canoncycles(c1)
    @test cyc_pow_perm(c, np) == pp
end

@testset "permsparse" begin
    p = makep()
    sp = permsparse(permmatrix(p))
    @test permlist(permcycles(sp)) == p

    sp = permsparse(permcycles(p))
    @test permlist(permmatrix(sp)) == p
end

@testset "permcompose" begin
    ntrials = 100
    for n in (3, 4, 10)
        for trial in 1:ntrials
            p = randperm(n);
            q = randperm(n);
            pm = permmatrix(p);
            qm = permmatrix(q);
            @test permcompose(p, q) == permlist(qm * pm)
        end
    end
end

@testset "permkron permcommute" begin
    ntrials = 100
    nbad = 0
    ngood = 0
#    for n in (3, 5, 10)
    for n in (3, 4)
        for trial in 1:ntrials
            p = randperm(n);
            q = randperm(n);
            pq = permkron(p, q);
            pm = permmatrix(p);
            qm = permmatrix(q);
            pqm = kron(pm, qm);
            @test pq == permlist(pqm)
            @test isperm(pq)
            if permcommute(p, q) != permcommute(pm, qm)
                nbad += 1
                # println(p)
                # println(q)
                # println()
            else
                ngood += 1
            end
        end
    end
    @info "$nbad / $(nbad + ngood) permcommute tests failed"
end


text:
module Button

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export btn, btngroup

Genie.Renderer.Html.register_normal_element("q__btn", context = @__MODULE__)
Genie.Renderer.Html.register_normal_element("q__btn__group", context = @__MODULE__)

function btn( label::String = "",
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              content::Union{String,Vector,Function} = "",
              kwargs...)

  wrap() do
    q__btn([isa(content, Function) ? content() : join(content)],
            args...;
            attributes([:label => label, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function btn( content::Function,
              label::String = "",
              args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  btn(label, args...; wrap = wrap, content = content, kwargs...)
end

function btn( args...;
              label::String = "",
              content::Union{String,Vector,Function} = "",
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  btn(label, args...; wrap = wrap, content = content, kwargs...)
end

const button = btn

function btngroup(args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)
  wrap() do
    q__btn__group(args...; kwargs...)
  end
end

const buttongroup = btngroup

end


text:
function WrapDist(x,y)
    # Computes the wrap-around distance on the torus [0,1).
    #
    # Inputs:
    # x: first frequency
    # y: second frequency
    #
    # Output:
    # wrap around distance between x and y

    return min(abs(x-y),abs(x-y+1))
end


text:
"""
Provides logging functionality for Genie apps.
"""
module Loggers

using Dates
using Genie
using Logging

import Base.log
export log, initialize_logging

global _isloggingconfigured = false

function _default_logcolor(level)
    level < Logging.Info  ? Base.debug_color() :
    level < Logging.Warn  ? Base.info_color()  :
    level < Logging.Error ? Base.warn_color()  :
                    Base.error_color()
end
"""
  Customized logging format function.
"""
function metafmt(level, _module, group, id, file, line)
  # Set the prefix and suffix for the logging message here.
  color = _default_logcolor(level)
  prefix = "[$(string(now()))] $(level == Logging.Warn ? "Warning" : string(level)):"
  suffix = ""
  Logging.Info <= level < Logging.Warn && return color, prefix, suffix
  _module !== nothing && (suffix *= "$(_module)")
  if file !== nothing
    _module !== nothing && (suffix *= " ")
    suffix *= Base.contractuser(file)
    if line !== nothing
      suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
    end
  end
  !isempty(suffix) && (suffix = "@ " * suffix)
  return color, prefix, suffix
end

function _logstringtologgingenum(level::Union{String,Symbol})::Logging.LogLevel
  level = string(level)
  level == "debug" && return Logging.Debug
  level == "info" && return Logging.Info
  level == "warn" && return Logging.Warn
  level == "error" && return Logging.Error
  # Default
  return Logging.Info
end

"""
Initialize the logging.
"""
function initialize_logging()::Nothing
  global _isloggingconfigured
  if !_isloggingconfigured
    @info "Configuring a global logger to log level: '$(Genie.config.log_level)'. Message after this will be filtered, you can change this in your config file for your environment."
    # Default to stderr.
    #TODO: Can set to write to a stream here with the first parameter if file logging should be done.
    # Personally I just keep everything working on console so that console loggers like Splunk and Kibana can just capture them
    consoleLogger = ConsoleLogger(stderr, _logstringtologgingenum(Genie.config.log_level), meta_formatter=metafmt)
    global_logger(consoleLogger)
    _isloggingconfigured = true
  end
  return nothing
end

"""
    log(message, level = "info"; showst::Bool = true) :: Nothing
    log(message::Any, level::Any = "info"; showst::Bool = false) :: Nothing
    log(message::String, level::Symbol) :: Nothing

Logs `message` to all configured logs (STDOUT, FILE, etc) by delegating to `Lumberjack`.
Supported values for `level` are "info", "warn", "debug", "error".
If `level` is `error` or `critical` it will also dump the stacktrace onto STDOUT.

# Examples
```julia
```
"""
function log(message::Union{String,Symbol,Number,Exception}, level::Union{String,Symbol} = "info") :: Nothing
  message = string(message)
  level = string(level)

  # Just a check in case it wasn't initialized.
  global _isloggingconfigured
  if !_isloggingconfigured
    initializeLogging()
  end

  try
    level == "debug" && @debug message
    level == "info" && @info message
    level == "warn" && @warn message
    (level == "error" || level == "err" || level == "critical") && @error message
  catch ex
    @info string(ex)
  end

  nothing
end

"""
    truncate_logged_output(output::String) :: String

Truncates (shortens) output based on `output_length` settings and appends "..." -- to be used for limiting the output length when logging.

# Examples
```julia
julia> Genie.config.output_length
100

julia> Genie.config.output_length = 10
10

julia> Loggers.truncate_logged_output("abc " ^ 10)
"abc abc ab..."
```
"""
function truncate_logged_output(output::String) :: String
  length(output) > Genie.config.output_length && output[1:Genie.config.output_length] * "..."
end

function log_path(path = Genie.LOG_PATH) :: String
  try
    "$(joinpath(path, haskey(ENV, "GENIE_ENV") ? ENV["GENIE_ENV"] : "dev")).log"
  catch ex
    string(ex) |> println
    # println("...")
    "$(joinpath(path, "dev")).log"
  end
end
function log_path!(path = Genie.LOG_PATH) :: String
  if ! isfile(log_path(path))
    mkpath(path)
    touch(log_path(path))
  end

  log_path(path)
end

"""
    empty_log_queue() :: Vector{Tuple{String,Symbol}}

The Genie log queue is used to push log messages in the early phases of framework bootstrap,
when the logger itself is not available. Once the logger is ready, the queue is emptied and the
messages are logged.
Automatically invoked.
"""
function empty_log_queue() :: Nothing
  for log_message in Genie.GENIE_LOG_QUEUE
    log(log_message...)
  end

  empty!(Genie.GENIE_LOG_QUEUE)

  nothing
end

"""
    macro location()

Provides a macro that injects the FILE and the LINE where the logger was invoked.
"""
macro location()
  :(log(" in $(@__FILE__):$(@__LINE__)", :err))
end

# setup_loggers()
# empty_log_queue()

end


text:
module Aiyagari

using QuantEcon
using Interpolations
#using Optim
using Parameters
using LinearAlgebra: norm
using SparseArrays: sparse
using Arpack: eigs
using StructArrays
using ProgressMeter
using Reexport

include("Expectations.jl")
@reexport using .Expectations

include("ExogenousStates.jl")
@reexport using .ExogenousStates

include("expectations-new.jl")

abstract type HouseState end
struct IsState <: HouseState end
struct NoState <: HouseState end

abstract type Household end
abstract type CoupledHouseholds end

@with_kw struct Consumer{T} <: Household
  𝔼::T = Unconditional()
end

@with_kw struct Owner{T1,T2<:HouseState} <: Household
  𝔼::T1 = Unconditional()
  state::T2 = NoState()
end

@with_kw struct Renter{T} <: Household
  𝔼::T = Unconditional()
end

𝔼(hh::Household) = hh.𝔼

@with_kw struct OwnOrRent{O<:Owner,R<:Renter} <: CoupledHouseholds
  owner::O = Owner()
  renter::R = Renter()
end

households(chh::OwnOrRent) = [chh.owner, chh.renter]
Base.length(chh::CoupledHouseholds) = length(households(chh))

include("bellman.jl")
include("stationary-distribution.jl")
include("aggregate-state.jl")

export solve_bellman
export controlled_markov_chain, stationary_distribution
export AggregateState
export Household, Owner, Renter, Consumer, OwnOrRent
export HouseState, IsState, NoState

export MarkovChain

examples = ["huggett", "housing", "morestates"]

end # module