text:
using RocksDB

const UNINITIALIZED = 0
const FOLLOWER = 1
const CANDIDATE = 2
const LEADER = 3

current_state = UNINITIALIZED

struct RavanaException <: Exception
    msg::AbstractString
end

mutable struct raftInfo
    currentTerm # latest term server has seen
    votedFor    # candidateId that received vote in current term
    commitIndex # index of highest log entry known to be committed
    lastApplied # index of highest log entry applied to state machine
    lastIndex   # index of last entry in log
end

struct raftNode
    name           # IP Address of name of node
    port           # Port where Raft protocol server is active
    nodeId         # Node id of the Raft node
    clusterId      # Cluster that the node belongs to
end

log       = nothing # Handle to the persistent log
nodes     = nothing # Array of nodes in cluster
nodeId    = nothing # This node id
clusterId = nothing # Cluster id
ipAddress = nothing # Protocol server IP
ipPort    = nothing # Protocol server port
rInfo     = nothing # Raft node state

leaderAddress = nothing
leaderPort    = nothing

const TASK_KILL = 1
mutable struct followerState
    nextIndex::Int # Index of next entry to replicate
    t::Task        # Task syncying up follower
    flag::Int      # To communicate with task
end

followers = Dict{UInt128, followerState}() # List of followers
clients = Dict{Int, Condition}()  # List of clients waiting on log replication
repCond = Condition()        # Condition to wake up replicator tasks
tasks = Dict{String, Task}() # Dict of first level tasks

const CANDIDATE_SLEEP_GRANULARITY = 0.001
const HEARTBEAT_TIMEOUT = 0.01
# Election timeout must be
#   1. > heartbeat timeout
#   2. > worst one way communication time to another node
const ELECTION_TIMEOUT  = 0.1

raft_increment_term() = (rInfo.currentTerm = rInfo.currentTerm + 1)
set_current_term(t) = (rInfo.currentTerm = t)
set_last_index(i) = (rInfo.lastIndex = i)
set_commit_index(index) = (rInfo.commitIndex = index)

function raft_uninitialized()
    return
end

function raft_follower()
    @info ("Moved to FOLLOWER state")
    start_election_timer()
end

function raft_candidate()
    @info ("Moved to CANDIDATE state")
    votes = 0
    for n in nodes
        @async begin  # Parallel requests for votes
            @debug "Calling request vote"
            (remoteTerm, status) = remote_request_vote(n)
            status && (votes += 1)
        end
    end
    while !is_majority(votes) && (raft_get_state() == CANDIDATE)
        sleep(HEARTBEAT_TIMEOUT)
    end
    raft_get_state() != CANDIDATE && return # State changed asynchronously
    is_majority(votes) && raft_set_state(LEADER)
    !is_majority(votes) && raft_set_state(FOLLOWER)
end

function raft_leader()
    @info "Moved to LEADER state"
    global leaderAddress = ipAddress
    global leaderPort = ipPort
    @info "leader: $(ipAddress):$(leaderPort)"
    start_heartbeat()
end

const raft_states = Dict(UNINITIALIZED   =>  raft_uninitialized,
                         FOLLOWER        =>  raft_follower,
                         CANDIDATE       =>  raft_candidate,
                         LEADER          =>  raft_leader)

function state_transition()
    f = raft_states[current_state]
    f()
end

raft_get_state() = (return current_state)

function raft_set_state(s)
    cur = raft_get_state()
    cur == s && return
    if (cur == LEADER) && (s != LEADER)
        leaderAddress = nothing
        leaderPort = nothing
    end
    global current_state = s
    state_transition()
end

function init_raft()
    current_state = UNINITIALIZED
    read_state()
    raft_server()
    ravana_server()
    raft_init_protocol()
end

function read_state()
    init_raft_log()
    global nodes = get_cluster_nodes()
    if nodes == nothing
        println("Cluster not initialized.")
        println("1. Call init_cluster() to start a new cluster with this machine as the first node")
        println("2. Call add_node_to_cluster(clusterId) to add this machine to an existing cluster")
        return
    end

    init_raft_info()
    init_cluster_state()
end

function get_cluster_nodes()
    db_get(log, "RavanaCluster")
end

function init_cluster_state()
    (nId, cluId)  = raft_get_node()
    global clusterId = cluId
    global nodeId = nId
    for i in nodes
        if i.nodeId == nodeId
            @assert cluId == i.clusterId
            global ipAddress = i.name
            global ipPort    = i.port
        end
    end

    raft_set_state(FOLLOWER)
    @info ("Node is part of cluster: ", clusterId)
end

function init_raft_info()
    global rInfo = raft_get_info()
end

function raft_init_protocol()
    state_machine_watcher()
    # state_watcher()
end

# A task to watch for changes to state of this node and
# trigger actions
const STATE_CHANGE_TIMEOUT = 0.001
function state_watcher()
    @async begin
        while true
            state_transition()
            sleep(STATE_CHANGE_TIMEOUT)
        end
    end
end

const STATE_MACHINE_TIMEOUT = 1
# A task to apply log entries to state machine
function state_machine_watcher()
    t = @task begin
        while true
            if raft_get_state() == UNINITIALIZED
                sleep(STATE_MACHINE_TIMEOUT)
                continue
            end
            if (rInfo.commitIndex > rInfo.lastApplied)
                execute_entry(rInfo.lastApplied + 1)
                rInfo.lastApplied += 1
            end
            raft_save_info(rInfo)
            sleep(STATE_MACHINE_TIMEOUT)
        end
    end
    tasks["state_machine_watcher"] = t
    schedule(t)
end

# Applies log entry at index to state machine
function execute_entry(index)
    (term, op, argv) = get_log_entry(index)
    raft_execute(op, op_table(op), argv)
end

# get_heartbeat(), set_heartbeat() and start_election_timer() called
# in the context of follower
receivedHeartbeat = false
get_heartbeat() = return receivedHeartbeat
set_heartbeat!(val) = (global receivedHeartbeat = val)

function start_election_timer()
    t = @task while true
        sleep_random()
        (raft_get_state() != FOLLOWER) && continue
        if !get_heartbeat()
            @debug ("No heartbeat received")
            raft_increment_term()
            raft_set_state(CANDIDATE)
        end
        set_heartbeat!(false)
    end
    tasks["election_timer"] = t
    schedule(t)
end

"""
    start_heartbeat()
Start heartbeat tasks for each follower
"""
function start_heartbeat()
    for n in nodes # Init
        (n.nodeId != nodeId) && init_replicator_task(n)
    end
    t = @task while true # Wake up replicators periodically
        sleep(HEARTBEAT_TIMEOUT)
        (raft_get_state() != LEADER) && continue
        notify(repCond, all=true)
    end
    tasks["log_replicator"] = t
    schedule(t)
end

function init_replicator_task(node)
    haskey(followers, node.nodeId) && return
    t = @task replicate(node)
    followers[node.nodeId] = fState = followerState(rInfo.lastIndex, t, 0)
    schedule(t)
end

"""
    replicate(node)

For the given node
1. If all local log entries are replicated, send NOOP entry as heartbeat
2. Else start sending entries that are not replicated
"""
function replicate(node)
    term = rInfo.currentTerm
    while true
        wait(repCond)  # Wait for work
        (raft_get_state() != LEADER) && continue # Only leader replicates
        if fState.nextIndex == (rInfo.lastIndex + 1)
            r = remote_append_entries(rInfo.currentTerm, OP_NOOP, nothing, rInfo.lastIndex, rInfo.currentTerm, rInfo.commitIndex, node)
        else
            replicate_entries(node)
        end
    end
end

"""
    replicate_entries(node)

For given node replicate log entries that are missing from leader (this node)
"""
function replicate_entries(node)
    fState = followers[node.nodeId]
    (prevTerm, pevOp, prevArgv) = get_log_entry(fState.nextIndex - 1)
    while (fState.nextIndex <= rInfo.lastIndex) && (raft_get_state() == LEADER)
        (term, op, argv) = get_log_entry(fState.nextIndex)
        (fTerm, fLastIndex, status) = remote_append_entries(term, op, argv, fState.nextIndex - 1, prevTerm, rInfo.commitIndex, node)
        if !status
            if fTerm > rInfo.currentTerm  # Some other leader superseded
                rInfo.currentTerm = term
                leaderAddress = node.name
                leaderPort = node.port
                @async raft_set_state(FOLLOWER)
                break
            end
            fState.nextIndex = fLastIndex
        else
            fState.nextIndex += 1
        end
        prevTerm = term
    end
end

"""
    raft_run_command(op, func, args)

Run the command from a client. First replicate to a majority,
then execute the command and return the result. This is the only
time commitIndex will be incremented - see the safety argument in
the paper.
"""
function raft_run_command(op, func, args)
    (raft_get_state() != LEADER) && (return nothing)
    r = get_log_entry(rInfo.lastIndex)
    @assert(r != nothing)
    (prevTerm, op, argv) = r
    payload = (rInfo.currentTerm, nodeId, rInfo.lastIndex, prevTerm, rInfo.commitIndex, op, argv)
    setup_condition(rInfo.lastIndex + 1) # Setup condition variable on index
    @debug "Append local log"
    local_append_entries(payload)
    @debug "Wake up tasks to append to other nodes' logs"
    notify(repCond, all=true)
    wait_on_majority(rInfo.lastIndex + 1)
    @debug "Replicated to majority! Execute command on this node"
    set_commit_index(rInfo.lastIndex + 1)
    ret = raft_execute(op, func, args)
    ret
end

"""
    setup_condition(index)

Setup condition variable for the index in the log. The condition
is released when a majority of nodes have replicated the entry at *index*.
Clients can wait on these to return a success to the caller of a command.
"""
function setup_condition(index)
    @assert(haskey(clients, index) == false)
    (length(nodes) == 1) && return
    clients[index] = Condition()
    return
end

"""
    wait_on_majority(index)

Waits for the log entry at index on this node to be replicated to a majority
of nodes in the cluster.
"""
function wait_on_majority(index)
    !haskey(clients, index) && return
    cond = clients[index]
    while !majority_replicated(index)
        wait(cond)
    end
    delete!(clients, index)
    return
end

"""
    majority_replicated(index)

Returns true if the log entry at *index* has been replicated to
a majority of nodes.
"""
function majority_replicated(index)
    votes = 1
    for n in nodes
        fState = followers[node.nodeId]
        (fState.nextIndex > index) && (votes += 1)
    end
    is_majority(votes) && return true
    false
end

"""
    remote_append_entries()
op         : Operation to replicate
argv       : Arguments to OP
prevIndex  : Index in log of entry just before current
prevTerm   : Term of prevIndex
commitIndex: Index of log that is committed on the leader
node       : Raft node to which this is sent

Called by leader to send log entries to followers.
"""
function remote_append_entries(term, op, argv, prevIndex, prevTerm, commitIndex, node)
    @info ("append_entries() term=", term, " op=", op, " prevIndex=", prevIndex, " term=", prevTerm, "node=", node)
    payload = (term, rInfo.nodeId, prevIndex, prevTerm, commitIndex, op, argv)
    r = raft_client(node.name, node.port, RAFT_APPEND_ENTRIES, payload)
end

"""
    local_append_entries()
Called on the followers to add log entries from leader to their logs.
"""
function local_append_entries(argv)
    (currentTerm, leaderNodeId, prevIndex, prevTerm, commitIndex, op, args) = argv

    # Add entry without question if this node is leader
    if leaderNodeId == nodeId
        @debug "Adding entry to leader/local log"
        (op != OP_NOOP) && add_log_entry(prevIndex + 1, currentTerm, op, args)
        rInfo.lastIndex = prevIndex + 1
        return (rInfo.currentTerm, prevIndex, true)
    end

    # If this node is superseded by another, become a follower
    if currentTerm >= rInfo.currentTerm
        leader = get_node_info(rInfo.nodeId)
        @debug "Found a another leader $(leader.name) superseding $(leaderAddress)"
        leaderAddress = leader.name
        leaderPort = leader.port
        @async raft_set_state(FOLLOWER)
        return (term, prevIndex, false)
    end

    entry = get_log_entry(prevIndex)

    # Don't apply entry in the following cases
    # 1. Entry preceding the current one is not written to log
    entry == nothing && return (rInfo.currentTerm, rInfo.lastIndex, false)
    (term, op, args) = entry
    # 2. Preceding entry's term does not match with leader
    if prevTerm != term
        delete_log_to_end(prevIndex)  # Delete all entries from prevIndex
        rInfo.lastIndex = prevIndex - 1
        return (term, prevIndex-1, false)
    end

    # Success case
    (op != OP_NOOP) && add_log_entry(prevIndex + 1, currentTerm, op, args)
    set_commit_index(min(commitIndex, prevIndex+1))
    set_current_term(currentTerm)
    set_last_index(prevIndex + 1)
    global received_heartbeat = true
    return (rInfo.currentTerm, prevIndex, true)
end

"""
    remote_request_vote()

Called by candidates in order to request votes to become the leader.
Messages are sent to all nodes in the Raft cluster.
"""
function remote_request_vote(node)
    (term, op, args) = get_log_entry(rInfo.commitIndex)
    payload = (rInfo.currentTerm, nodeId, rInfo.commitIndex, term)
    r = raft_client(node.name, node.port, RAFT_REQUEST_VOTE, payload)
    r
end

votedFor  = nothing
votedTerm = 1
"""
    local_request_vote()

Receives a request_vote request from a candidate and processes it.
"""
function local_request_vote(argv)
    (newTerm, leaderId, lastIndex, lastTerm) = argv

    # Grant vote to self unconditionally
    leaderId == nodeId && return (rInfo.currentTerm, true)

    # Deny vote for these conditions
    # 1. Candidate's term is < this node's term
    newTerm < rInfo.currentTerm && return (rInfo.currentTerm, false)
    # 2. This node voted for a node who's term is >= candidate's term
    votedTerm >= newTerm && return (rInfo.currentTerm, false)
    # 3. This node's log has more entries than candidate's log
    rInfo.lastIndex > lastIndex && (rInfo.currentTerm, false)

    # Success case follows
    votedTerm = newTerm
    # If this node is superseded by another, become a follower
    if newTerm >= rInfo.currentTerm
        leader = get_node_info(leaderId)
        leaderAddress = leader.name
        leaderPort = leader.port
        @async raft_set_state(FOLLOWER)
    end

    return (rInfo.currentTerm, true)
end

#=
   ## Raft log utils ##
=#
function init_raft_log()
    global log = RocksDB.open_db(basedir() * "raft_log", true)
end

function get_log_entry(index)
    r = db_get(log, index)
    r
end

function add_log_entry(index, term, op, args)
    @assert(rInfo.lastIndex == index-1)
    db_put(log, index, (term, op, args))
end

"""
    find_last_log_index(hint)
Find the index of the last valid entry in the log.
"""
function find_last_log_index(hint)
    @assert(get_log_entry(rInfo.lastApplied) != nothing)
    search_log(rInfo.lastApplied, hint)
end

# Binary search log for last entry
function search_log(first, last)
    mid = floor(Int, (first + last) / 2)
    if first == mid
        if get_log_entry(last) != nothing
            return last
        else
            return first
        end
    end
    if get_log_entry(mid + 1) != nothing
        search_log(mid, last)
    else
        search_log(first, mid)
    end
end

# Delete all entries in log beginning with index
function delete_log_to_end(index)
    while get_log_entry(index) != nothing
        db_delete(log, index)
        index += 1
    end
end

# Save some persistent variables
function raft_save_info(st::raftInfo)
    db_put(log, "RaftInfo", st)
end

# Get persistent variables
function raft_get_info()
    db_get(log, "RaftInfo")
end


text:
"""
This module implements a new type of `AbstractString`, a `StringView`,
that provides a string representation of any underlying array of bytes
(any `AbstractVector{UInt8}`), interpreted as UTF-8 encoded Unicode data.

Unlike Julia's built-in `String` type (which also wraps UTF-8 data), the
`StringView` type is a copy-free wrap of *any* `AbstractVector{UInt8}`
instance, and does not take "ownership" of or modify the array.   Otherwise,
a `StringView` is intended to be usable in any context where you might
have otherwise used `String`.
"""
module StringViews
export StringView

"""
    StringView{T<:AbstractVector{UInt8}} <: AbstractString

`StringView(array)` creates an `AbstractString` representation of
any `array` of `UInt8` data, interpreted as UTF-8 encoded Unicode.
It does *not* make a copy of or modify `array`.
"""
struct StringView{T<:AbstractVector{UInt8}} <: AbstractString
    data::T
end

const DenseStringView = StringView{<:Union{DenseVector{UInt8},<:Base.FastContiguousSubArray{UInt8,1,<:DenseVector{UInt8}}}}
const StringAndSub = Union{String,SubString{String}}
const StringViewAndSub = Union{StringView,SubString{<:StringView}}
const DenseStringViewAndSub = Union{DenseStringView,SubString{<:DenseStringView}}

Base.Vector{UInt8}(s::StringView{Vector{UInt8}}) = s.data
Base.Vector{UInt8}(s::StringViewAndSub) = Vector{UInt8}(codeunits(s))
Base.Array{UInt8}(s::StringViewAndSub) = Vector{UInt8}(s)
Base.String(s::StringViewAndSub) = String(copyto!(Base.StringVector(ncodeunits(s)), codeunits(s)))
StringView(s::StringView) = s
StringView(s::String) = StringView(codeunits(s))

Base.Symbol(s::DenseStringViewAndSub) =
    return ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), s, ncodeunits(s))

Base.pointer(s::DenseStringView) = pointer(s.data)
Base.pointer(s::DenseStringView, i::Integer) = pointer(s.data, i)
Base.pointer(x::SubString{<:DenseStringView}) = pointer(x.string) + x.offset
Base.pointer(x::SubString{<:DenseStringView}, i::Integer) = pointer(x.string) + x.offset + (i-1)
Base.unsafe_convert(::Type{Ptr{UInt8}}, s::DenseStringViewAndSub) = pointer(s)
Base.unsafe_convert(::Type{Ptr{Int8}}, s::DenseStringViewAndSub) = convert(Ptr{Int8}, pointer(s))
Base.cconvert(::Type{Ptr{UInt8}}, s::SubString{<:DenseStringView}) = s
Base.cconvert(::Type{Ptr{Int8}}, s::SubString{<:DenseStringView}) = s

Base.sizeof(s::StringView) = length(s.data)
Base.ncodeunits(s::StringView) = length(s.data)
Base.codeunit(s::StringView) = UInt8
Base.@propagate_inbounds Base.codeunit(s::StringView, i::Integer) = s.data[i]
Base.codeunits(s::StringView) = s.data
Base.codeunits(s::SubString{<:StringView}) = @view s.string.data[1+s.offset:s.offset+s.ncodeunits]

_memcmp(a, b, len) =
    ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), a, b, len % Csize_t) % Int
function _cmp(a, b)
    al, bl = sizeof(a), sizeof(b)
    c = _memcmp(a, b, min(al,bl))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(al,bl)
end
Base.cmp(a::DenseStringViewAndSub, b::DenseStringViewAndSub) = _cmp(a, b)
Base.cmp(a::DenseStringViewAndSub, b::StringAndSub) = _cmp(a, b)
Base.cmp(a::StringAndSub, b::DenseStringViewAndSub) = _cmp(a, b)
Base.:(==)(s1::StringViewAndSub, s2::StringViewAndSub) = codeunits(s1) == codeunits(s2)
Base.:(==)(s1::StringAndSub, s2::StringViewAndSub) = codeunits(s1) == codeunits(s2)
function Base.:(==)(a::StringAndSub, b::DenseStringViewAndSub)
    al = sizeof(a)
    return al == sizeof(b) && 0 == _memcmp(a, b, al)
end
Base.:(==)(s1::StringViewAndSub, s2::StringAndSub) = s2 == s1

Base.typemin(::Type{StringView{Vector{UInt8}}}) = StringView(Vector{UInt8}(undef,0))
Base.typemin(::T) where {T<:StringView} = typemin(T)

Base.isvalid(s::DenseStringViewAndSub) = ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), s, sizeof(s)) ≠ 0
Base.isvalid(s::StringViewAndSub) = all(isvalid, s)
Base.isvalid(::Type{String}, s::StringViewAndSub) = isvalid(s)

function Base.isascii(s::StringViewAndSub)
    @inbounds for i = 1:ncodeunits(s)
        codeunit(s, i) >= 0x80 && return false
    end
    return true
end

Base.write(io::IO, s::StringViewAndSub) = write(io, codeunits(s))
Base.print(io::IO, s::StringViewAndSub) = (write(io, s); nothing)

Base.@propagate_inbounds Base.thisind(s::StringViewAndSub, i::Int) = Base._thisind_str(s, i)
Base.@propagate_inbounds Base.nextind(s::StringViewAndSub, i::Int) = Base._nextind_str(s, i)
Base.isvalid(s::StringViewAndSub, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i

function Base.hash(s::DenseStringViewAndSub, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s, ncodeunits(s), h % UInt32) + h
end

# each string type must implement its own reverse because it is generally
# encoding-dependent
function Base.reverse(s::StringViewAndSub)::String
    # Read characters forwards from `s` and write backwards to `out`
    out = Base._string_n(sizeof(s))
    offs = sizeof(s) + 1
    for c in s
        offs -= ncodeunits(c)
        Base.__unsafe_string!(out, c, offs)
    end
    return out
end

include("decoding.jl")
include("regex.jl")
include("parse.jl")
include("util.jl")
include("search.jl")

end # module

ext: jl


text:
#function Base.checkbounds() end

function Base.checkindex(::Type{Bool}, idx::AbstractIndex, i)
    return _checkindex(index_by(idx, i), idx, i)
end
function Base.checkindex(::Type{Bool}, idx::AbstractIndex, i::AbstractVector)
    return _checkindex(index_by(idx, i), idx, i)
end
function Base.checkindex(::Type{Bool}, idx::AbstractIndex, i::AbstractRange)
    return _checkindex(index_by(idx, i), idx, i)
end
function Base.checkindex(::Type{Bool}, idx::AbstractIndex, i::Real)
    return _checkindex(index_by(idx, i), idx, i)
end

_checkindex(::ByKeyTrait, idx, i) = !isnothing(find_first(==(i), keys(idx)))
function _checkindex(::ByKeyTrait, idx, i::AbstractVector)
    for ii in i
        ii in keys(idx) || return false
    end
    return true
end
_checkindex(::ByValueTrait, idx, i) = checkindex(Bool, values(idx), i)
#=
function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::K) where {name,K}
    return !isnothing(find_first(==(i), keys(indx)))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,Int}, i::Int) where {name}
    return !isnothing(find_first(==(i), keys(indx)))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,CartesianIndex{1}}, i::CartesianIndex{1}) where {name}
    return !isnothing(find_first(==(i), keys(indx)))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::Int) where {name,K}
    return !isnothing(find_first(==(i), keys(values(indx))))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::CartesianIndex{1}) where {name,K}
    return !isnothing(find_first(==(i), keys(values(indx))))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::AbstractVector{K}) where {name,K}
    return !isnothing(find_first(==(i), keys(indx)))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,CartesianIndex{1}}, i::AbstractVector{CartesianIndex{1}}) where {name}
    return !isnothing(find_first(==(i), keys(indx)))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::AbstractVector{Int}) where {name,K}
    return !isnothing(find_first(==(i), keys(values(indx))))
end

function Base.checkindex(::Type{Bool}, indx::AbstractIndex{name,K}, i::AbstractVector{CartesianIndex{1}}) where {name,K}
    return !isnothing(find_first(==(i), keys(values(indx))))
end
=#


text:
module TestExamples

import Modia3D

# Test Modia3D - Examples

# test/dynamics
dynamicsPath = joinpath(Modia3D.path, "test", "dynamics")
include(joinpath(dynamicsPath, "Simulate_ControllerDamper.jl"))
include(joinpath(dynamicsPath, "Simulate_DamperMacro.jl"))
include(joinpath(dynamicsPath, "Simulate_FourBar.jl"))
include(joinpath(dynamicsPath, "Simulate_pathPlanning.jl"))
include(joinpath(dynamicsPath, "Simulate_PendulumWithFixedJoint.jl"))

# test/inertiaTensor
inertiaPath = joinpath(Modia3D.path, "test", "inertiaTensor")
include(joinpath(inertiaPath, "Test_2Rev_ZylZ_BarX.jl"))
include(joinpath(inertiaPath, "Test_3Rev_ZylZ_BarX_BarY.jl"))
include(joinpath(inertiaPath, "Test_InertiaTensor.jl"))

# test/revoluteJoints
testJointsPath = joinpath(Modia3D.path, "test", "joints")
include(joinpath(testJointsPath, "Test_KinematicRevoluteJoints.jl"))
include(joinpath(testJointsPath, "Test_Dynamic_Pendulum_xAxis.jl"))
include(joinpath(testJointsPath, "Test_Dynamic_Pendulum_yAxis.jl"))
include(joinpath(testJointsPath, "Test_Dynamic_Pendulum_zAxis.jl"))
include(joinpath(testJointsPath, "Test_Prismatic_xAxis.jl"))
include(joinpath(testJointsPath, "Test_Prismatic_yAxis.jl"))
include(joinpath(testJointsPath, "Test_Prismatic_zAxis.jl"))
include(joinpath(testJointsPath, "Move_FourBar_xAxis.jl"))
include(joinpath(testJointsPath, "Move_FourBar_yAxis.jl"))
include(joinpath(testJointsPath, "Move_FourBar_zAxis.jl"))

# test/kinematics
kinematicsPath = joinpath(Modia3D.path, "test", "kinematics")
include(joinpath(kinematicsPath, "Move_FourBar_noMacros.jl"))

# test/massComputation
massCompPath = joinpath(Modia3D.path, "test", "massComputation")
include(joinpath(massCompPath, "test_massComputation.jl"))

# test/signalToFlange
signalToFlangePath = joinpath(Modia3D.path, "test", "signalToFlange")
include(joinpath(signalToFlangePath, "Test_Signal1Assembly.jl"))
include(joinpath(signalToFlangePath, "Test_Signal4Assemblies.jl"))

# test/visual
visuPath = joinpath(Modia3D.path, "test", "visual")
include(joinpath(visuPath, "visual_3Bars.jl"))

# test/volumeComputation
volCompPath = joinpath(Modia3D.path, "test", "volumeComputation")
include(joinpath(volCompPath, "volume_computation3D_obj.jl"))

# test/withoutAssembly/kinematics
wAKinematicsPath = joinpath(Modia3D.path, "test", "withoutAssembly","kinematics")
include(joinpath(wAKinematicsPath, "Move_Pendulum.jl"))

# examples/withoutAssembly/visual
wAVisualPath = joinpath(Modia3D.path, "test", "withoutAssembly","visual")
include(joinpath(wAVisualPath, "Visualize_Beam.jl"))

println("\n ...test_Examples finished!")
end


text:
"""
    Constraint{T}

    Abstract type for graph edge object
"""
abstract type Constraint{T,N} <: Node{T} end
Base.length(::Constraint{T,N}) where {T,N} = N


text:
using Test

@testset "GenerateJuliaPrograms" begin
  include("GenerateJuliaPrograms.jl")
end

@testset "ParseJuliaPrograms" begin
  include("ParseJuliaPrograms.jl")
end

@testset "ParseRelations" begin
  include("ParseRelations.jl")
end


text:
module TestResilience
    using Test, EcologicalNetworks

    A = rand([0, 1, 0.5], 10, 10)
    N = UnipartiteQuantitativeNetwork(A)

    @assert s(N) ≈ sum(A) / 10
    @assert s(N, dims=1) ≈ sum(A, dims=2)
    @assert s(N, dims=2) ≈ sum(A, dims=1)'

    @assert resilience(N) ≈ s(N) + symmetry(N) * heterogeneity(N)

    # species only interact with themselves
    Nid = UnipartiteNetwork([true false false false;
                              false true false false;
                              false false true false;
                              false false false true])

    @assert heterogeneity(Nid) ≈ 0.0
    # @assert symmetry(Nid)   NAN

    # symetric network
    Nsym = UnipartiteNetwork([false true true true;
                               true false false false;
                               true false false false;
                               true false false false])
    @assert symmetry(Nsym) ≈ 1.0  # in- and out degree is same
    @assert heterogeneity(Nsym) ≈ 0.5

    # asymetric network
    Nasym = UnipartiteNetwork([false true true true;
                               false false false false;
                               false false false false;
                               false false false false])
    @assert symmetry(Nasym) ≈ -1.0  # in- and out degree is same
    @assert heterogeneity(Nasym) ≈ 0.75

    v = ones(10) / 10
    Nhomo = UnipartiteProbabilisticNetwork(v * v')
    @assert isapprox(heterogeneity(Nhomo), 0.0, atol=1e-10)
end


text:
get_free_memory_size(model::TierraModel) = model.memory_size - model.used_memory
normalize_address(x) = mod(x, 2^16)

function read_memory(model::TierraModel, address)
    model.memory[mod(address, model.memory_size)+one(UInt16)]
end
function read_memory(model::TierraModel, address, length)
    return [
        read_memory(model, curr_addr) for
        curr_addr = address:(address+length-one(UInt16))
    ]
end
function write_memory(
    model::TierraModel,
    organism::TierrianOrganism,
    address,
    instruction::UInt8,
)
    if organism_has_access(organism, address, 1)
        if MUTATE && rand() < FLAW_PROBABILITY
            instruction_to_write = get_mutated_instruction(instruction)
        else
            instruction_to_write = instruction
        end

        model.memory[mod(address, model.memory_size)+one(UInt16)] = instruction_to_write
    else # organism has no write access to this address
        organism.error_flag = true
    end
end
function write_memory(
    model::TierraModel,
    organism::TierrianOrganism,
    start_address,
    instructions::Vector{UInt8},
)
    for (i, instrucion) in enumerate(instructions)
        write_memory(model, organism, start_address + i - 1, instrucion)
    end
end

function organism_has_access(organism::TierrianOrganism, address, length)
    if _is_within_range(
        address,
        length,
        organism.start_address,
        organism.length,
    )
        true
    elseif organism.has_daughter && _is_within_range(
        address,
        length,
        organism.daughter_address,
        organism.daughter_length,
    )
        true
    else
        false
    end
end

function clear_memory!(model::TierraModel, start_address, length)
    length = min(length, model.memory_size - start_address)

    # first, we test if there are existing overlapping blocks

    for i = 1:Base.length(model.free_blocks)
        free_block = model.free_blocks[i]

        overlap = _determine_overlap(
            start_address,
            length,
            free_block.start_address,
            free_block.length,
        )

        if overlap in [ExactMatch(), Within(), WithinFromStart(), WithinToEnd()]
            return 0

        elseif overlap == LeftOverlap()
            return clear_memory!(
                model,
                start_address,
                free_block.start_address - start_address,
            )

        elseif overlap == RightOverlap()
            return clear_memory!(
                model,
                free_block.start_address + free_block.length,
                start_address + length - free_block.start_address -
                free_block.length,
            )

        elseif overlap == CompleteOverlap()
            return sum([
                clear_memory!(
                    model,
                    start_address,
                    free_block.start_address - start_address,
                ),
                clear_memory!(
                    model,
                    free_block.start_address + free_block.length,
                    start_address + length - free_block.start_address -
                    free_block.length,
                ),
            ])
        end
    end

    # if there is no overlapping block, we serach for potential neighboring free blocks to merge

    index_left_neighbor = -1
    index_right_neighbor = -1

    for i = 1:Base.length(model.free_blocks)
        free_block = model.free_blocks[i]

        overlap = _determine_overlap(
            start_address,
            length,
            free_block.start_address,
            free_block.length,
        )

        if overlap == LeftNeighbor()
            @assert index_left_neighbor == -1
            index_left_neighbor = i
        elseif overlap == RightNeighbor()
            @assert index_right_neighbor == -1
            index_right_neighbor = i
        end
    end

    if index_left_neighbor != -1 && index_right_neighbor == -1
        free_block = model.free_blocks[index_left_neighbor]

        model.free_blocks[index_left_neighbor] =
            FreeMemoryBlock(start_address, free_block.length + length)
        model.used_memory -= length

        return length

    elseif index_left_neighbor == -1 && index_right_neighbor != -1
        free_block = model.free_blocks[index_right_neighbor]

        model.free_blocks[index_right_neighbor] = FreeMemoryBlock(
            free_block.start_address,
            free_block.length + length,
        )
        model.used_memory -= length

        return length

    elseif index_left_neighbor != -1 && index_right_neighbor != -1
        free_block_right = model.free_blocks[index_left_neighbor] # the free block where the block to clear is the left neighbor is its right neighbor!!
        free_block_left = model.free_blocks[index_right_neighbor] # TODO

        model.free_blocks[index_left_neighbor] = FreeMemoryBlock(
            free_block_left.start_address,
            free_block_left.length + free_block_right.length + length,
        )
        deleteat!(model.free_blocks, index_right_neighbor)

        model.used_memory -= length

        return length
    end

    # no overlapping or neighboring free memory blocks => create a new one

    push!(model.free_blocks, FreeMemoryBlock(start_address, length))
    model.used_memory -= length
    return length
end

function allocate_free_memory!(model::TierraModel, length)
    if get_free_memory_size(model) < length # no memory left
        return -1
    end

    for _ = 1:100
        block_index = rand(1:Base.length(model.free_blocks))
        block = model.free_blocks[block_index]

        if length <= block.length
            alloc_start, new_free_block =
                _allocate_in_block!(model, block, length)

            if new_free_block === nothing || new_free_block.length == 0
                deleteat!(model.free_blocks, block_index)
            else
                model.free_blocks[block_index] = new_free_block
            end

            model.used_memory += length

            return alloc_start
        end
    end

    return -1
end

function _allocate_in_block!(model::TierraModel, block::FreeMemoryBlock, length)
    if block.length == length
        return block.start_address, nothing
    end

    block_start = block.start_address

    if rand(["beginning", "end"]) == "beginning"
        alloc_start = block_start
        alloc_end = block_start + length
        new_free_block = FreeMemoryBlock(alloc_end, block.length - length)
    else
        alloc_start = block_start + block.length - length
        new_free_block = FreeMemoryBlock(block_start, block.length - length)
    end

    return UInt16(alloc_start), new_free_block
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: SymbolicContinuousCallback, SymbolicContinuousCallbacks, NULL_AFFECT, get_callback

@parameters t
@variables x(t)=0
D = Differential(t)

eqs = [D(x) ~ 1]
affect = [x ~ 0]

## Test SymbolicContinuousCallback
@testset "SymbolicContinuousCallback constructors" begin
    e = SymbolicContinuousCallback(eqs[])
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    e = SymbolicContinuousCallback(eqs)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    e = SymbolicContinuousCallback(eqs, NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    e = SymbolicContinuousCallback(eqs[], NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    e = SymbolicContinuousCallback(eqs => NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    e = SymbolicContinuousCallback(eqs[] => NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT

    ## With affect

    e = SymbolicContinuousCallback(eqs[], affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    e = SymbolicContinuousCallback(eqs, affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    e = SymbolicContinuousCallback(eqs, affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    e = SymbolicContinuousCallback(eqs[], affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    e = SymbolicContinuousCallback(eqs => affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    e = SymbolicContinuousCallback(eqs[] => affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect

    # test plural constructor

    e = SymbolicContinuousCallbacks(eqs[])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == NULL_AFFECT

    e = SymbolicContinuousCallbacks(eqs)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == NULL_AFFECT

    e = SymbolicContinuousCallbacks(eqs[] => affect)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks(eqs => affect)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks([eqs[] => affect])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks([eqs => affect])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks(SymbolicContinuousCallbacks([eqs => affect]))
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

end

##

@named sys = ODESystem(eqs, continuous_events = [x ~ 1])
@test getfield(sys, :continuous_events)[] == SymbolicContinuousCallback(Equation[x ~ 1], NULL_AFFECT)
@test isequal(equations(getfield(sys, :continuous_events))[], x ~ 1)
fsys = flatten(sys)
@test isequal(equations(getfield(fsys, :continuous_events))[], x ~ 1)

@named sys2 = ODESystem([D(x) ~ 1], continuous_events = [x ~ 2], systems=[sys])
@test getfield(sys2, :continuous_events)[] == SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT)
@test all(ModelingToolkit.continuous_events(sys2) .== [SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT), SymbolicContinuousCallback(Equation[sys.x ~ 1], NULL_AFFECT)])

@test isequal(equations(getfield(sys2, :continuous_events))[1], x ~ 2)
@test length(ModelingToolkit.continuous_events(sys2)) == 2
@test isequal(ModelingToolkit.continuous_events(sys2)[1].eqs[], x ~ 2)
@test isequal(ModelingToolkit.continuous_events(sys2)[2].eqs[], sys.x ~ 1)

# Functions should be generated for root-finding equations
prob = ODEProblem(sys, Pair[], (0.0, 2.0))
p0 = 0
t0 = 0
@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.ContinuousCallback
cb = ModelingToolkit.generate_rootfinding_callback(sys)
cond = cb.condition
out = [0.0]
cond.rf_ip(out, [0], p0, t0)
@test out[] ≈ -1 # signature is u,p,t
cond.rf_ip(out, [1], p0, t0)
@test out[] ≈ 0  # signature is u,p,t
cond.rf_ip(out, [2], p0, t0)
@test out[] ≈ 1  # signature is u,p,t

prob = ODEProblem(sys, Pair[], (0.0, 2.0))
sol = solve(prob, Tsit5())
@test minimum(t->abs(t-1), sol.t) < 1e-10 # test that the solver stepped at the root

# Test that a user provided callback is respected
test_callback = DiscreteCallback(x->x, x->x)
prob = ODEProblem(sys, Pair[], (0.0, 2.0), callback = test_callback)
cbs = get_callback(prob)
@test cbs isa CallbackSet
@test cbs.discrete_callbacks[1] == test_callback

prob = ODEProblem(sys2, Pair[], (0.0, 3.0))
cb = get_callback(prob)
@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback

cond = cb.condition
out = [0.0, 0.0]
# the root to find is 2
cond.rf_ip(out, [0,0], p0, t0)
@test out[1] ≈ -2 # signature is u,p,t
cond.rf_ip(out, [1,0], p0, t0)
@test out[1] ≈ -1  # signature is u,p,t
cond.rf_ip(out, [2,0], p0, t0) # this should return 0
@test out[1] ≈ 0  # signature is u,p,t

# the root to find is 1
out = [0.0, 0.0]
cond.rf_ip(out, [0,0], p0, t0)
@test out[2] ≈ -1 # signature is u,p,t
cond.rf_ip(out, [0,1], p0, t0) # this should return 0
@test out[2] ≈ 0  # signature is u,p,t
cond.rf_ip(out, [0,2], p0, t0)
@test out[2] ≈ 1  # signature is u,p,t

sol = solve(prob, Tsit5())
@test minimum(t->abs(t-1), sol.t) < 1e-10 # test that the solver stepped at the first root
@test minimum(t->abs(t-2), sol.t) < 1e-10 # test that the solver stepped at the second root

@named sys = ODESystem(eqs, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same state
prob = ODEProblem(sys, Pair[], (0.0, 3.0))
@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
sol = solve(prob, Tsit5())
@test minimum(t->abs(t-1), sol.t) < 1e-10 # test that the solver stepped at the first root
@test minimum(t->abs(t-2), sol.t) < 1e-10 # test that the solver stepped at the second root

## Test bouncing ball with equation affect
@variables t x(t)=1 v(t)=0
D = Differential(t)

root_eqs = [x ~ 0]
affect   = [v ~ -v]

@named ball = ODESystem([
    D(x) ~ v
    D(v) ~ -9.8
], t, continuous_events = root_eqs => affect)

@test getfield(ball, :continuous_events)[] == SymbolicContinuousCallback(Equation[x ~ 0], Equation[v ~ -v])
ball = structural_simplify(ball)

tspan = (0.0,5.0)
prob = ODEProblem(ball, Pair[], tspan)
sol = solve(prob,Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
# plot(sol)

## Test bouncing ball in 2D with walls
@variables t x(t)=1 y(t)=0 vx(t)=0 vy(t)=1
D = Differential(t)

continuous_events = [
    [x ~ 0] => [vx ~ -vx]
    [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]
]

@named ball = ODESystem([
    D(x) ~ vx
    D(y) ~ vy
    D(vx) ~ -9.8
    D(vy) ~ -0.01vy # there is some small air resistance
], t, continuous_events = continuous_events)

ball = structural_simplify(ball)

tspan = (0.0,5.0)
prob = ODEProblem(ball, Pair[], tspan)

cb = get_callback(prob)
@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
@test getfield(ball, :continuous_events)[1] == SymbolicContinuousCallback(Equation[x ~ 0], Equation[vx ~ -vx])
@test getfield(ball, :continuous_events)[2] == SymbolicContinuousCallback(Equation[y ~ -1.5, y ~ 1.5], Equation[vy ~ -vy])
cond = cb.condition
out = [0.0, 0.0, 0.0]
cond.rf_ip(out, [0,0,0,0], p0, t0)
@test out ≈ [0, 1.5, -1.5]

sol = solve(prob,Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
@test minimum(sol[y]) ≈ -1.5 # check wall conditions
@test maximum(sol[y]) ≈ 1.5  # check wall conditions

# tv = sort([LinRange(0, 5, 200); sol.t])
# plot(sol(tv)[y], sol(tv)[x], line_z=tv)
# vline!([-1.5, 1.5], l=(:black, 5), primary=false)
# hline!([0], l=(:black, 5), primary=false)

## Test multi-variable affect
# in this test, there are two variables affected by a single event.
continuous_events = [
    [x ~ 0] => [vx ~ -vx, vy ~ -vy]
]

@named ball = ODESystem([
    D(x) ~ vx
    D(y) ~ vy
    D(vx) ~ -1
    D(vy) ~ 0
], t, continuous_events = continuous_events)

ball = structural_simplify(ball)

tspan = (0.0,5.0)
prob = ODEProblem(ball, Pair[], tspan)
sol = solve(prob,Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
@test -minimum(sol[y]) ≈ maximum(sol[y]) ≈ sqrt(2)  # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number)

# tv = sort([LinRange(0, 5, 200); sol.t])
# plot(sol(tv)[y], sol(tv)[x], line_z=tv)
# vline!([-1.5, 1.5], l=(:black, 5), primary=false)
# hline!([0], l=(:black, 5), primary=false)


text:
module MATLAB

using Base.Sys: islinux, iswindows, isapple
using Libdl
using SparseArrays

import Base: eltype, close, size, copy, ndims, unsafe_convert

# mxarray
export MxArray, mxClassID, mxComplexity,
       mxclassid, data_ptr,
       classid, nrows, ncols, nelems, elsize

export is_double, is_single,
       is_int8, is_uint8, is_int16, is_uint16,
       is_int32, is_uint32, is_int64, is_uint64,
       is_numeric, is_complex, is_sparse, is_empty,
       is_logical, is_char, is_struct, is_cell

export mxarray, mxsparse, delete,
       mxcellarray, get_cell, set_cell,
       mxstruct, mxstructarray, mxnfields, get_fieldname, get_field, set_field,
       jvalue, jarray, jscalar, jvector, jmatrix, jsparse, jstring, jdict

# engine & matfile
export MSession, MatFile,
       get_default_msession, restart_default_msession, close_default_msession,
       eval_string, get_mvariable, get_variable, put_variable, put_variables,
       variable_names, read_matfile, write_matfile,
       mxcall,
       @mput, @mget, @mat_str

if iswindows()
    export show_msession, hide_msession, get_msession_visiblity
end

# exceptions
struct MEngineError <: Exception
    message::String
end

include("init.jl") # initialize Refs
include("mxbase.jl")
include("mxarray.jl")
include("matfile.jl")
include("engine.jl")
include("matstr.jl")

if iswindows()
    # workaround "primary message table for module 77" error
    # creates a dummy Engine session and keeps it open so the libraries used by all other
    # Engine clients are not loaded and unloaded repeatedly
    # see: https://www.mathworks.com/matlabcentral/answers/305877-what-is-the-primary-message-table-for-module-77

    # initialization is delayed untill first call to MSession
    const persistent_msession_ref = Ref{MSession}()
    const persistent_msession_assigned = Ref(false)

    function assign_persistent_msession()
        if persistent_msession_assigned[] == false
            persistent_msession_assigned[] = true
            persistent_msession_ref[] = MSession(0)
        end
    end
end

function __init__()

    # load libraries

    libmx[]  = Libdl.dlopen(joinpath(matlab_libpath(), "libmx"), Libdl.RTLD_GLOBAL)
    libmat[] = Libdl.dlopen(joinpath(matlab_libpath(), "libmat"), Libdl.RTLD_GLOBAL)
    libeng[] = Libdl.dlopen(joinpath(matlab_libpath(), "libeng"), Libdl.RTLD_GLOBAL)

    # engine functions

    eng_open[]          = engfunc(:engOpen)
    eng_close[]         = engfunc(:engClose)
    eng_set_visible[]   = engfunc(:engSetVisible)
    eng_get_visible[]   = engfunc(:engGetVisible)
    eng_output_buffer[] = engfunc(:engOutputBuffer)
    eng_eval_string[]   = engfunc(:engEvalString)
    eng_put_variable[]  = engfunc(:engPutVariable)
    eng_get_variable[]  = engfunc(:engGetVariable)

    # mxarray functions

    mx_destroy_array[]   = mxfunc(:mxDestroyArray)
    mx_duplicate_array[] = mxfunc(:mxDuplicateArray)

    # load functions to access mxarray

    mx_free[]         = mxfunc(:mxFree)

    mx_get_classid[]  = mxfunc(:mxGetClassID)
    mx_get_m[]        = mxfunc(:mxGetM)
    mx_get_n[]        = mxfunc(:mxGetN)
    mx_get_nelems[]   = mxfunc(:mxGetNumberOfElements)
    mx_get_ndims[]    = mxfunc(:mxGetNumberOfDimensions_730)
    mx_get_elemsize[] = mxfunc(:mxGetElementSize)
    mx_get_data[]     = mxfunc(:mxGetData)
    mx_get_dims[]     = mxfunc(:mxGetDimensions_730)
    mx_get_nfields[]  = mxfunc(:mxGetNumberOfFields)
    mx_get_pr[]       = mxfunc(:mxGetPr)
    mx_get_pi[]       = mxfunc(:mxGetPi)
    mx_get_ir[]       = mxfunc(:mxGetIr_730)
    mx_get_jc[]       = mxfunc(:mxGetJc_730)

    mx_is_double[]    = mxfunc(:mxIsDouble)
    mx_is_single[]    = mxfunc(:mxIsSingle)
    mx_is_int64[]     = mxfunc(:mxIsInt64)
    mx_is_uint64[]    = mxfunc(:mxIsUint64)
    mx_is_int32[]     = mxfunc(:mxIsInt32)
    mx_is_uint32[]    = mxfunc(:mxIsUint32)
    mx_is_int16[]     = mxfunc(:mxIsInt16)
    mx_is_uint16[]    = mxfunc(:mxIsUint16)
    mx_is_int8[]      = mxfunc(:mxIsInt8)
    mx_is_uint8[]     = mxfunc(:mxIsUint8)
    mx_is_char[]      = mxfunc(:mxIsChar)

    mx_is_numeric[]   = mxfunc(:mxIsNumeric)
    mx_is_logical[]   = mxfunc(:mxIsLogical)
    mx_is_complex[]   = mxfunc(:mxIsComplex)
    mx_is_sparse[]    = mxfunc(:mxIsSparse)
    mx_is_empty[]     = mxfunc(:mxIsEmpty)
    mx_is_struct[]    = mxfunc(:mxIsStruct)
    mx_is_cell[]      = mxfunc(:mxIsCell)

    # load functions to create & delete MATLAB array

    mx_create_numeric_matrix[]   = mxfunc(:mxCreateNumericMatrix_730)
    mx_create_numeric_array[]    = mxfunc(:mxCreateNumericArray_730)

    mx_create_double_scalar[]  = mxfunc(:mxCreateDoubleScalar)
    mx_create_logical_scalar[] = mxfunc(:mxCreateLogicalScalar)

    mx_create_sparse[]         = mxfunc(:mxCreateSparse_730)
    mx_create_sparse_logical[] = mxfunc(:mxCreateSparseLogicalMatrix_730)

    mx_create_string[]         = mxfunc(:mxCreateString)
    mx_create_char_array[]     = mxfunc(:mxCreateCharArray_730)

    mx_create_cell_array[]     = mxfunc(:mxCreateCellArray_730)

    mx_create_struct_matrix[]  = mxfunc(:mxCreateStructMatrix_730)
    mx_create_struct_array[]   = mxfunc(:mxCreateStructArray_730)

    mx_get_cell[]              = mxfunc(:mxGetCell_730)
    mx_set_cell[]              = mxfunc(:mxSetCell_730)

    mx_get_field[]             = mxfunc(:mxGetField_730)
    mx_set_field[]             = mxfunc(:mxSetField_730)
    mx_get_field_bynum[]       = mxfunc(:mxGetFieldByNumber_730)
    mx_get_fieldname[]         = mxfunc(:mxGetFieldNameByNumber)

    mx_get_string[]            = mxfunc(:mxGetString_730)

    # load I/O mat functions

    mat_open[]         = matfunc(:matOpen)
    mat_close[]        = matfunc(:matClose)
    mat_get_variable[] = matfunc(:matGetVariable)
    mat_put_variable[] = matfunc(:matPutVariable)
    mat_get_dir[]      = matfunc(:matGetDir)

end

###########################################################
#
#   deprecations
#
###########################################################

end


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: BoundingBox, get_overlapping_fields, infer_init
import Celeste: Log, SDSSIO

const usage_info =
"""
Usage:
  estimate-box-runtime.jl <ramin> <ramax> <decmin> <decmax>
"""

const stagedir = ENV["CELESTE_STAGE_DIR"]

const max_hardness = 10_000

function to_tasks(box)
    rcfs = get_overlapping_fields(box, stagedir)
    catalog = SDSSIO.read_photoobj_files(rcfs, stagedir,
                        duplicate_policy=:primary)
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_sources = find(entry_in_range, catalog)

    num_rcfs, num_targets = length(rcfs), length(target_sources)

    hardness = num_rcfs * num_targets
    if hardness <= max_hardness
        box_str = "$(box.ramin) $(box.ramax) $(box.decmin) $(box.decmax)"
        println("$hardness\t$num_rcfs\t$num_targets\t$box_str")
    else
        sl = (box.ramax - box.ramin) / 2
        to_tasks(BoundingBox(box.ramin, box.ramin + sl,
                             box.decmin, box.decmin + sl))
        to_tasks(BoundingBox(box.ramin + sl, box.ramax,
                             box.decmin, box.decmin + sl))
        to_tasks(BoundingBox(box.ramin, box.ramin + sl,
                             box.decmin + sl , box.decmax))
        to_tasks(BoundingBox(box.ramin + sl, box.ramax,
                             box.decmin + sl, box.decmax))
    end
end

if length(ARGS) != 4
    println(usage_info)
else
    to_tasks(BoundingBox(ARGS...))
end


text:
using Split, DataFrames
using Base.Test

# Array test
a = [1,2,3,4,5]

# DataArray test data
da = @data([1,2,3,4,5])
(size_da,) = size(da)
indices_da = [1:size_da;]

# DataFrames test data
df = DataFrame(A = 1:5, B = ["A", "B", "C", "D", "E"])

# Test data
fraction = 0.5
seed = 1234

# Test that the shuffled input is different than the output
@test Split.seeded_shuffle(indices_da, seed) != indices_da

# Test that the shuffled input is equal to itself when the seed is the same
@test Split.seeded_shuffle(indices_da, seed) == Split.seeded_shuffle(indices_da, seed)

# Test that the fractional split is correct
(a, b) = Split.sample_indices(size_da, fraction, seed)
@test size(a) == (2,) && size(b) == (3,)

# Test that the type for array is correct
@test typeof(Split.sample(a, fraction, seed)) <: Tuple{Array, Array}

# Test that the type for array is correct
@test typeof(Split.sample(da, fraction, seed)) <: Tuple{DataArray, DataArray}

# Test that the type for array is correct
@test typeof(Split.sample(df, fraction, seed)) <: Tuple{DataFrame, DataFrame}


text:
# This file was generated, do not modify it. # hide
y, X = unpack(df, ==(:Scaled_Sound));

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module CCallbacks

function parse_cargs(arg_spec,arg_sel,do_escape)
    ret_type = arg_spec.args[2]
    cargs = arg_spec.args[1].args[arg_sel]

    arg_names = []
    arg_types = []

    for carg in cargs
        carg.head == :(::) || error("Invalid argument specification $(carg)")
        push!(arg_names, do_escape ? :($(esc(carg.args[1]))) : :($(carg.args[1])))
        push!(arg_types, carg.args[2])
    end

    arg_names,arg_types,ret_type
end

# Stolen from https://gist.github.com/simonbyrne/c4146dc286fd5387385ca911e8318509
@eval macro $(Symbol("ccall"))(expr)
    expr.head == :(::) && expr.args[1].head == :call || error("Invalid use of @ccall")

    arg_names,arg_types,ret_type = parse_cargs(expr, 2:lastindex(expr.args[1].args), true)
    fname = expr.args[1].args[1]

    tupexpr = :(())
    ccexpr = :(ccall($(esc(fname)), $(esc(ret_type)), $(esc(tupexpr))))
    append!(ccexpr.args, arg_names)
    append!(tupexpr.args, arg_types)

    ccexpr
end

macro ccallback(arg_spec, body)
    arg_spec.head == :(::) && arg_spec.args[1].head == :tuple ||
        error("Invalid argument specification")
    body.head == :block || error("Callback block missing")

    arg_names,arg_types,ret_type = parse_cargs(arg_spec, :, false)

    # TODO: Need to figure out how to escape variable references in
    # the body block.
    local ccallback_expr = :(() -> $(body))
    ccallback_expr.args[1].args = arg_names
    # For convenience, such that we don't have to remember to put
    # nothing:s at the end of all callbacks declared Cvoid.
    ret_type == :Cvoid && push!(ccallback_expr.args[2].args[2].args, nothing)

    cfun = Expr(:cfunction, Ptr{Cvoid},
                QuoteNode(ccallback_expr),
                ret_type,
                Expr(:call, GlobalRef(Core, :svec), arg_types...),
                QuoteNode(:ccall))
    esc(cfun)
end

export @ccall, @ccallback

end # module


text:
# This file was generated, do not modify it. # hide
train, test = MLJ.partition(MLJ.eachindex(targets), 0.70, rng=52)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "CovidSim"
version = v"0.8.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mrc-ide/covid-sim.git", "51959a733cbb856e4019031994deba2b1a7a9e4a"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/covid-sim/

if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ../patches/0001-fix-lib-cases.patch
fi

mkdir build && cd build/
# Build with GCC also for FreeBSD and macOS to use libgomp
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake \
    -DCMAKE_BUILD_TYPE=Release \
    ../src
make -j${nprocs}
mkdir -p "${bindir}"
cp "CovidSim${exeext}" "${bindir}/"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("CovidSim", :CovidSim)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to handle the configuration objects of PrettyTables
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export pretty_table_with_conf, clear_pt_conf!, set_pt_conf, set_pt_conf!

################################################################################
#                                   Printing
################################################################################

"""
    pretty_table_with_conf(conf::PrettyTablesConf, args...; kwargs...)

Call `pretty_table` using the default configuration in `conf`. The `args...` and
`kwargs...` can be the same as those passed to `pretty_tables`. Notice that all
the configurations in `kwargs...` will overwrite the ones in `conf`.

The object `conf` can be created by the function `set_pt_conf` in which the
keyword parameters can be any one supported by the function `pretty_table` as
shown in the following.

"""
function pretty_table_with_conf(conf::PrettyTablesConf, args...; kwargs...)

    # Copy the configuration so that the user object is not modified.
    _local_conf = deepcopy(conf)

    # Apply the new configurations to it.
    set_pt_conf!(_local_conf; kwargs...)

    # Get the named tuple with the configurations.
    nt = _conf_to_nt(_local_conf)

    # Print the table.
    pretty_table(args...; nt...)
end

################################################################################
#                                Set and clear
################################################################################

"""
    clear_pt_conf!(conf::PrettyTablesConf)

Clear all configurations in `conf`.

"""
function clear_pt_conf!(conf::PrettyTablesConf)
    empty!(conf.confs)
    return nothing
end

"""
    set_pt_conf(;kwargs...)

Create a new configuration object based on the arguments in `kwargs`.

"""
@inline function set_pt_conf(;kwargs...)
    conf = PrettyTablesConf()
    set_pt_conf!(conf; kwargs...)
    return conf
end

"""
    set_pt_conf!(conf; kwargs...)

Apply the configurations in `kwargs` to the object `conf`.

"""
@inline function set_pt_conf!(conf::PrettyTablesConf; kwargs...)
    for kw in kwargs
        conf.confs[kw[1]] = kw[2]
    end
    return nothing
end

################################################################################
#                              Private functions
################################################################################

"""
    _conf_to_nt(conf::PrettyTablesConf)

Convert the configuration object `conf` to a named tuple so that it can be
passed to `pretty_table`.

"""
@inline function _conf_to_nt(conf::PrettyTablesConf)
    # Get the named tuple with the configurations.
    dictkeys = (collect(keys(conf.confs))...,)
    dictvals = (collect(values(conf.confs))...,)
    nt = NamedTuple{dictkeys}(dictvals)
end


text:
# same letter of different case cancel
#
# could do this in place but my hunch is that it's expensive to delete in middle
# of a vector so let's use an output vector we push to one at a time, while
# iterating over the input

input = read("05.input", String) |> strip |> collect

cancels(a::Char, b::Char) = abs(a-b) == 32

function cancel!(output, next)
    if !isempty(output) && cancels(next, last(output))
        pop!(output)
    else
        push!(output, next)
    end
    output
end

cancel(input) = reduce(cancel!, input, init=Char[])

cancel(collect("dabAcCaCBAcCcaDA")) |> String == "dabCBAcaDA"

# star 1
canceled = cancel(input)
length(canceled)

# star 2:
minimum((length(cancel(Iterators.filter(!in((t, uppercase(t))), input))), t)
         for t
         in unique(lowercase(c) for c in input))

# in-place version: don't delete, just move (and resize)
function canceli!(input)
    # the last index in the working area
    j = 0
    for i in eachindex(input)
        if j > 0 && cancels(input[j], input[i])
            j -= 1
        else
            j += 1
            input[j] = input[i]
        end
    end
    resize!(input, j)
    input
end

canceli!(collect("dabAcCaCBAcCcaDA")) |> String == "dabCBAcaDA"

using BenchmarkTools

input = read("05.input", String) |> strip |> collect

@btime canceli!(copy($input));
#  280.754 μs (2 allocations: 195.39 KiB)
@btime cancel($input);
#  572.529 μs (15 allocations: 128.69 KiB)


text:
# Form clique-reduction graph
function ExtractTriMat(A::SparseMatrixCSC{Float64,Int64})
    @assert(issymmetric(A))
    n = size(A,1)
    Am = zeros(n,n)
    for i = 1:n-2
        for j = i+1:n-1
            Aji = A[j,i]
            if Aji == 0
                continue
            end

            for k = j+1:n

                if A[k,j]+A[k,i] == 2
                    Am[j,i]+=1/2
                    Am[k,i]+=1/2
                    Am[k,j]+=1/2
                end
            end
        end
    end
    Am = sparse(Am + Am')
    return Am
end

# Form clique-reduction graph
function ExtractAllTriMat(A::SparseMatrixCSC{Float64,Int64})
    n = size(A,1)
    Am = zeros(n,n)
    for i = 1:n-2
        for j = i+1:n-1
            Aji = A[j,i]
            Aij = A[i,j]
            if Aji == 0 || Aij == 0
                continue
            end

            for k = j+1:n

                if A[k,j]+A[k,i]+A[j,k]+A[i,k] == 4
                    Am[j,i]+=1/2
                    Am[k,i]+=1/2
                    Am[k,j]+=1/2
                end
            end
        end
    end
    Am = sparse(Am + Am')
    return Am
end

# Extract triangles as a hyperedge list
function ExtractTriangle_List(A::SparseMatrixCSC{Float64,Int64})
    n = size(A,1)
    H1 = Vector{Int64}()
    H2 = Vector{Int64}()
    H3 = Vector{Int64}()
    for i = 1:n-2
        for j = i+1:n-1
            Aji = A[j,i]
            Aij = A[i,j]
            if Aji == 0 || Aij == 0
                continue
            end

            for k = j+1:n

                if A[k,j]+A[k,i]+A[j,k]+A[i,k] == 4
                    push!(H1,i)
                    push!(H2,j)
                    push!(H3,k)
                end
            end
        end
    end
    H = [H1 H2 H3]
    return H
end

function ExtractFanMat(A::SparseMatrixCSC{Float64,Int64})
    n = size(A,1)
    Am = zeros(n,n)
    for i = 1:n
        for j = i+1:n
            for k = j+1:n
                for l = k+1:n
                    # Yuck, quad loop.
                    ps= [i j k l; i k j l; i l j k; j k i l; j l i k; k l i j]
                    for w = 1:size(ps,1)
                        a = ps[w,1]; b = ps[w,2]; c = ps[w,3]; d = ps[w,4]
                        if A[a,c] > 0 && A[b,c] > 0 && A[a,d] > 0 && A[b,d] > 0 && A[c,a] ==0 && A[d,a] ==0 && A[c,b] ==0 && A[d,b] ==0
                            Am[a,b] += 1/2
                            Am[b,a] += 1/2
                            Am[c,d] += 1/2
                            Am[d,c] += 1/2
                        end
                    end
                end
            end
        end
    end
    Am = sparse(Am+Am')
    return Am
end

function ExtractFanList(A::SparseMatrixCSC{Float64,Int64})
    n = size(A,1)
    H1 = Vector{Int64}()
    H2 = Vector{Int64}()
    H3 = Vector{Int64}()
    H4 = Vector{Int64}()
    for i = 1:n
        for j = i+1:n
            for k = j+1:n
                for l = k+1:n
                    # Yuck, quad loop.
                    ps= [i j k l; i k j l; i l j k; j k i l; j l i k; k l i j]
                    for w = 1:size(ps,1)
                        a = ps[w,1]; b = ps[w,2]; c = ps[w,3]; d = ps[w,4]
                        if A[a,c] > 0 && A[b,c] > 0 && A[a,d] > 0 && A[b,d] > 0 && A[c,a] ==0 && A[d,a] ==0 && A[c,b] ==0 && A[d,b] ==0
                            push!(H1,a)
                            push!(H2,b)
                            push!(H3,c)
                            push!(H4,d)
                        end
                    end
                end
            end
        end
    end
    H = [H1 H2 H3 H4]
    return H
end

function turntounweighted(Aorig::SparseMatrixCSC{Float64,Int64})
    A = copy(Aorig)
    n = size(A,1)
    for i = 1:n
        for j= 1:n
            if A[i,j] > 0
                A[i,j] = 1
                A[j,i] = 1
            end
        end
    end
    return A
end

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
module MountainCar
# defines move_car! for many different types

function move_car!(x::Vector, u::Vector, dt::Number)
    x[2] += 0.001*u[1] - 0.0025*cos(3*x[1])

    # Correct bounds
    if (x[2] < -0.07)
        x[2] = -0.07
    elseif (x[2]>0.07)
        x[2] = 0.07
    end

    # Update position
    if x[1] >= 0.5
        x[2] = 0
        return x
    end

    x[1] += x[2]

    #check if the car reached a boundary
    if (x[1]<=-1.2)
        x[1] = -1.2;
        x[2] = 0.0;
    end
    return x
end

#"""
#x[1] position
#x[2] velocity
#"""
function cost_car(x::Vector, u::Vector)
    terminal = 0

    #cost = -(-0.01+0.0025/3.0*( sin(3.0*x[1])-sin(3.0*0.5) )); #continous reward
    cost = -10*(x[1]+0.5)^2 - 1*x[2]^2

    if (x[1]>=0.5)
        cost     = cost-100
        terminal = 1
    end

    return cost, terminal
end

rk4_func  = move_car!  # set standard Environment function
cost_func = cost_car # set standard Environment cost-function
s_init = [-pi/6 , 0.0]  #set initial position for environment
u_range = [1.0]

export rk4_func, cost_func, s_init, u_range

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
(f::NLopt.Algorithm)() = f

function __map_optimizer_args(prob::OptimizationProblem, opt::NLopt.Opt;
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    local_method::Union{NLopt.Algorithm, NLopt.Opt, Nothing} = nothing,
    local_maxiters::Union{Number, Nothing} = nothing,
    local_maxtime::Union{Number, Nothing} = nothing,
    local_options::Union{NamedTuple,Nothing} = nothing,
    kwargs...)

    if local_method !== nothing
        if isa(local_method,NLopt.Opt)
            if ndims(local_method) != length(prob.u0)
                error("Passed local NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
            end
            local_meth = local_method
        else
            local_meth = NLopt.Opt(local_method, length(prob.u0))
        end

        if !isnothing(local_options)
            for j in Dict(pairs(local_options))
                eval(Meta.parse("NLopt."*string(j.first)*"!"))(local_meth, j.second)
            end
        end

        if !(isnothing(local_maxiters))
            NLopt.maxeval!(local_meth, local_maxiters)
        end

        if !(isnothing(local_maxtime))
            NLopt.maxtime!(local_meth, local_maxtime)
        end

        NLopt.local_optimizer!(opt, local_meth)
    end

    # add optimiser options from kwargs
    for j in kwargs
        eval(Meta.parse("NLopt."*string(j.first)*"!"))(opt, j.second)
    end

    if prob.ub !== nothing
        NLopt.upper_bounds!(opt, prob.ub)
    end

    if prob.lb !== nothing
        NLopt.lower_bounds!(opt, prob.lb)
    end

    if !(isnothing(maxiters))
        NLopt.maxeval!(opt, maxiters)
    end

    if !(isnothing(maxtime))
        NLopt.maxtime!(opt, maxtime)
    end

    if !isnothing(abstol)
        NLopt.ftol_abs!(opt, abstol)
    end
    if !isnothing(reltol)
        NLopt.ftol_rel!(opt, reltol)
    end

    return nothing
end

function __solve(prob::OptimizationProblem, opt::Union{NLopt.Algorithm, NLopt.Opt};
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 local_method::Union{NLopt.Algorithm, NLopt.Opt, Nothing} = nothing,
                 local_maxiters::Union{Number, Nothing} = nothing,
                 local_maxtime::Union{Number, Nothing} = nothing,
                 local_options::Union{NamedTuple,Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false,
                 cb = (args...) -> (false),
                 kwargs...)
    local x

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)
    local_maxiters = _check_and_convert_maxiters(local_maxiters)
    local_maxtime = _check_and_convert_maxtime(local_maxtime)

    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p)

    _loss = function(θ)
        x = f.f(θ, prob.p)
        cb(θ,x...)
        return x[1]
    end

    fg! = function (θ,G)
        if length(G) > 0
            f.grad(G, θ)
        end

        return _loss(θ)
    end

    if isa(opt,NLopt.Opt)
        if ndims(opt) != length(prob.u0)
            error("Passed NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
        end
        opt_setup = opt
    else
        opt_setup = NLopt.Opt(opt, length(prob.u0))
    end

    prob.sense === MaxSense ? NLopt.max_objective!(opt_setup, fg!) : NLopt.min_objective!(opt_setup, fg!)

    _map_optimizer_args(prob,opt_setup, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol, local_method=local_method, local_maxiters =local_maxiters, local_options=local_options; kwargs...)

    t0 = time()
    (minf,minx,ret) = NLopt.optimize(opt_setup, prob.u0)
    t1 = time()

    SciMLBase.build_solution(prob, opt, minx, minf; original=opt_setup, retcode=ret)
end


text:
using GitHub

function init_metadata(packagespec, url)
    return Dict(
        "name" => packagespec.name,
        "uuid" => packagespec.uuid,
        "slug" => Base.package_slug(packagespec.uuid, 5),
        "url" => url,
        "version" => packagespec.version,
        "installable" => missing,
        "doctype" => missing,
        "success" => false,
        "documenter_errored" => missing
    )
end

function package_metadata(packagespec, url)
    meta = init_metadata(packagespec, url)

    rmatch = match(r".*/(.*)/(.*?(?:\.jl)?)(?:.git)?$", url)
    if rmatch !== nothing
        meta["owner"] = rmatch[1]
        meta["name"] = rmatch[2]
    end

    merge!(meta, update_metadata(packagespec, url, rmatch[1], rmatch[2]))

    return meta
end

function update_metadata(packagespec, url, repo_owner, repo_name)
    meta = Dict{String, Any}()

    authpath = if haskey(ENV, "DOCGEN_GITHUB_AUTH_FILE")
        ENV["DOCGEN_GITHUB_AUTH_FILE"]
    else
        joinpath(@__DIR__, "gh_auth.txt")
    end

    if !isfile(authpath)
        @warn("No GitHub token found. Skipping metadata retrieval.")
        return meta
    end

    if !occursin("github.com", url)
        @warn("Can't retrieve metadata (not hosted on GitHub).")
        return meta
    end

    @info("Querying metadata.")
    try
        gh_auth = authenticate(readchomp(authpath))
        repo_info = repo(repo_owner * "/" * repo_name, auth = gh_auth)
        meta["description"] = something(repo_info.description, "")
        meta["stargazers_count"]  = something(repo_info.stargazers_count, 0)
        meta["homepage"]  = something(repo_info.homepage, 0)
        topics_dict, page = topics(repo_info, auth = gh_auth)
        meta["tags"] = something(topics_dict["names"], [])
        meta["contributors"] = contributor_user.(contributors(repo_info, auth = gh_auth)[1])
    catch err
        @error(string("Couldn't get info for ", url), error = err)
    end
    @info("Done querying metadata.")

    return meta
end

function contributor_user(dict)
    Dict(
        "name" => dict["contributor"].login,
        "type" => dict["contributor"].typ,
        "url" => string(dict["contributor"].html_url),
        "contributions" => dict["contributions"]
    )
end

function topics(repo, api = GitHub.DEFAULT_API; options...)
    results, page_data = GitHub.gh_get_paged_json(api, "/repos/$(GitHub.name(repo))/topics";
                                                  headers = Dict("Accept" => "application/vnd.github.mercy-preview+json"), options...)
    return results, page_data
end


text:
using Flux: logitbinarycrossentropy
using Zygote: @adjoint
using ForwardDiff: derivative
using CuArrays
CuArrays.allowscalar(false)

prevw = rand(150000);
target = rand([0, 1.], length(prevw));

@adjoint logitbinarycrossentropy(w, t) = logitbinarycrossentropy(w, t), delta -> (delta*(1/(1+exp(-w)) - t), delta)

g(prevw, target) = Zygote.gradient(prevw->logitbinarycrossentropy(prevw, target), prevw)[1]
g2(prevw, target) = Zygote.gradient(logitbinarycrossentropy, prevw, target)[1]
g3(prevw, target) = ForwardDiff.gradient(x->logitbinarycrossentropy(x[1], target), [prevw])[1]
g4(prevw, target) = ForwardDiff.derivative(prevw->logitbinarycrossentropy(prevw, target), prevw)
g5(prevw, target) = Zygote.gradient(prevw) do prevw
	Zygote.forwarddiff(prew->logitbinarycrossentropy(prew, target), prevw)
end[1]
g6(prevw, target) = deriv(LogitProbLoss(), target, 1 / (1 + exp(-prevw)))
easy(prevw, target) = 1 / (1 + exp(-prevw)) - target

gprev = gpu(prevw)
gtarget = gpu(target)

g6.(gprev, gtarget)

res = logitbinarycrossentropy.(prevw, target)

value(LogitProbLoss(), target, 1 ./ (1 .+ exp.(-prevw)))

deriv(LogitProbLoss(), target, prevw)

deriv2(LogitProbLoss(), target, prevw)

using LossFunctions

value(LogitProbLoss(), target, prevw)

@benchmark g.($prevw, $target)
@benchmark g2.($prevw, $target)
@benchmark g3.($prevw, $target)
@benchmark g4.($prevw, $target)
@benchmark g5.($prevw, $target)
@benchmark easy.($prevw, $target)
@benchmark g6.($prevw, $target)


text:
using Dataton2019
using Test
using DataDeps, DataFrames, LightGraphs, CSV, ZipFile, JSON

@info "Hola Dataton2019!"
@testset "Datos" begin
    #global t1 = readdir(datadep"icijparadise2018")[1]
    #@test  "csv_paradise_papers.2018-02-14.zip" == t1
    #global t2 = readdir(datadep"icijpanamapapers2018")[1]
    #@test "csv_panama_papers.2018-02-14.zip" == t2
    global t3 = readdir(datadep"icijoffshore2018")[1]
    @test  "csv_offshore_leaks.2018-02-14.zip" == t3
    #global t4 = readdir(datadep"icijbahamas2017")[1]
    #@test  "csv_bahamas_leaks.2017-12-19.zip" == t4

    # Esta es una version mini de el archivo "contrataciones_arr.json"
    # Encontrado en https://github.com/pdnmx/datosdataton2019
    # Minimizamos el archivo con el comando 'split -n 100 contartaciones_arr.json'
    # y cambiamos el nombre de 'xaa' a 'mini.json'
    #
    # Encontramos la ultima entrada que se ha recortado y le pegamos un "]"
    #
    x = readlines(open("mini.json"))[1];
    #
    # julia> findlast("{\"_id\"",x)
    # 22869056:22869061
    #
    s6 = JSON.parse(x[1:22869054] * "]");

    @test length(s6) == 6462

end

# Hacemos el unzip del archivo de bahamas2017

#pars2018 = ZipFile.Reader(datadep"icijparadise2018" * "/" * t1)
#panama2018 = ZipFile.Reader(datadep"icijpanamapapers2018" * "/" * t2)
offshore2018 = ZipFile.Reader(datadep"icijoffshore2018" * "/" * t3)
#bahamas2017 = ZipFile.Reader(datadep"icijbahamas2017" * "/" * t4)

# Este es el formato de los archivos a investigar

# uncompressedsize method  mtime            name
# ----------------------------------------------
#        234614743 Deflate 2018-02-13 17-44 paradise_papers.edges.csv
#         47379240 Deflate 2018-02-13 17-40 paradise_papers.nodes.address.csv
#         64159484 Deflate 2018-02-13 17-40 paradise_papers.nodes.entity.csv
#           239615 Deflate 2018-02-13 17-40 paradise_papers.nodes.intermediary.csv
#         52548340 Deflate 2018-02-13 17-40 paradise_papers.nodes.officer.csv
#           371926 Deflate 2018-02-13 17-40 paradise_papers.nodes.other.csv

# pars2018address = CSV.read(pars2018.files[2]; threaded = true)
# pars2018entity = CSV.read(pars2018.files[3]; threaded = true)
# pars2018intermediares = CSV.read(pars2018.files[4]; threaded = true)
# pars2018officer = CSV.read(pars2018.files[5]; threaded = true)
# pars2018other = CSV.read(pars2018.files[6]; threaded = true)

# panama2018address = CSV.read(panama2018.files[1]; threaded = true)
# panama2018entity = CSV.read(panama2018.files[2]; threaded = true)
# panama2018interm = CSV.read(panama2018.files[3]; threaded = true)
# panama2018officer = CSV.read(panama2018.files[4]; threaded = true)

# bahamas2017address = CSV.read(bahamas2017.files[1]; threaded = true)
# bahamas2017entity = CSV.read(bahamas2017.files[2]; threaded = true)
# bahamas2017interm = CSV.read(bahamas2017.files[3]; threaded = true)
# bahamas2017officer = CSV.read(bahamas2017.files[4]; threaded = true)

# Extraemos a CSVs

offshore2018address = CSV.read(offshore2018.files[1]; threaded = true)
offshore2018entity = CSV.read(offshore2018.files[2]; threaded = true)
offshore2018interm = CSV.read(offshore2018.files[3]; threaded = true)
offshore2018officer = CSV.read(offshore2018.files[4]; threaded = true)

occursin.(offshore2018entity[:address], "Mexico") |> sum
occursin.(offshore2018entity[:address], "México") |> sum
occursin.(offshore2018entity[:address], "MX") |> sum
occursin.(offshore2018entity[:address], "Mx") |> sum
occursin.(offshore2018entity[:address], "MEX") |> sum


text:
# Now that you're warmed up, it's time to play the real game.
#
# A second compartment opens, this time labeled Dirac dice. Out of it falls a
# single three-sided die.
#
# As you experiment with the die, you feel a little strange. An informational
# brochure in the compartment explains that this is a quantum die: when you
# roll it, the universe splits into multiple copies, one copy for each possible
# outcome of the die. In this case, rolling the die always splits the universe
# into three copies: one where the outcome of the roll was 1, one where it was
# 2, and one where it was 3.
#
# The game is played the same as before, although to prevent things from
# getting too far out of hand, the game now ends when either player's score
# reaches at least 21.
#
# Using the same starting positions as in the example above, player 1 wins in
# 444356092776315 universes, while player 2 merely wins in 341960390180808
# universes.
#
# Using your given starting positions, determine every possible outcome. Find
# the player that wins in more universes; in how many universes does that
# player win?

# read input
positions = map(l -> Int8(l[end] - '0'), eachline("input.txt"))

# map from total of three rolls to the number of ways to get that total
rolls = Dict{Int8,Int}()
for t in [Int8(x+y+z) for x in 1:3, y in 1:3, z in 1:3]
  rolls[t] = get(rolls, t, 0) + 1
end

# map from a state (player positions and scores) to the number of ways to get
# to that state
states = Dict((positions...,0,0) => 1)

# For each state, calculate the next states by playing all possible moves for
# the current player's turn.
wins = zeros(Int, length(positions))
turn = 0
while length(states) > 0
  nextstates = Dict{keytype(states),valtype(states)}()
  for (state,n) in states
    for (roll, nrolls) in rolls
      (p1,p2,s1,s2) = state
      nextn = n * nrolls
      if (turn & 1) == 0
        p1 = (p1 + roll - 1) % 10 + 1
        s1 += p1
        if s1 >= 21
          wins[1] += nextn
          continue
        end
      else
        p2 = (p2 + roll - 1) % 10 + 1
        s2 += p2
        if s2 >= 21
          wins[2] += nextn
          continue
        end
      end
      k = (p1,p2,s1,s2)
      nextstates[k] = get(nextstates, k, 0) + nextn
    end
  end
  global states = nextstates
  global turn = xor(turn, 1)
end

println("Result: ", maximum(wins))


text:
using GeometryBasics
using GeometryBasics: Mesh
using Colors
#
# Meshes in 3d space
#

"Compute a mesh from points (vertices) and faces"
function mesh3d(ps,fs; vertexColors = nothing )
    pps = [Point(ps[i,1],ps[i,2],ps[i,3]) for i in 1:size(ps)[1]]
    ffs = [TriangleFace(fs[i,1],fs[i,2],fs[i,3]) for i in 1:size(fs)[1]]
    if(vertexColors != nothing)
        return  meta(Mesh(pps, ffs),vertexColors=vertexColors)
    end
    Mesh(pps, ffs)
end

"Compute Mesh from the coordinates of triangles"
function mesh3d(faces_coords::Array{<:Real,2})
    nps = Point{3, Float64}[]
    nfs = NgonFace{3,Int}[]
    normals = Point{3,Float64}[]
    vcnt = 1
    for i in 1:size(faces_coords,1)
        p1 = faces_coords[i,1:3]
        p2 = faces_coords[i,4:6]
        p3 = faces_coords[i,7:9]
        push!(nps, Point(p1[1],p1[2],p1[3]))
        push!(nps, Point(p2[1],p2[2],p2[3]))
        push!(nps, Point(p3[1],p3[2],p3[3]))
        push!(nfs, TriangleFace(vcnt,vcnt+1,vcnt+2))
        vcnt = vcnt + 3
    end
    Mesh(nps,nfs)
end

"Compute line segments from vertices coordinates and edges"
function edges3d(ps,es; color=RGBA(1,1,1,1), linewidth=1)
    xs = ps[:,1]
    ys = ps[:,2]
    zs = ps[:,3]
    xes = transpose(xs[es])
    yes = transpose(ys[es])
    zes = transpose(zs[es])
    xes = vec(xes)
    yes = vec(yes)
    zes = vec(zes)
    MeshCat.LineSegments(Point.(xes,yes,zes),LineBasicMaterial(color=color, linewidth=linewidth))
end

"Compute line segments from the edges coordinates"
function edges3d(es_ps::Array{<:Real,2}; color=RGBA(1,1,1,1), linewidth=1)
    cnt = 2*size(es_ps,1)
    x = zeros(cnt)
    y = zeros(cnt)
    z = zeros(cnt)
    x[1:2:cnt] = es_ps[:,1]
    x[2:2:cnt] = es_ps[:,4]
    y[1:2:cnt] = es_ps[:,2]
    y[2:2:cnt] = es_ps[:,5]
    z[1:2:cnt] = es_ps[:,3]
    z[2:2:cnt] = es_ps[:,6]
    mat = LineBasicMaterial(color=color, linewidth=linewidth)
    LineSegments(Point.(x,y,z),mat)
end


text:
module NonNegInts

export NNInt32, NNInt64

import Base: +, -, *, ÷, div, fld, cld, mod, rem

struct NNInt{T<:Union{UInt32, UInt64}} <: Unsigned
    value::T
end

const NNInt32 = NNInt{UInt32}
const NNInt64 = NNInt{UInt64}

function NNInt32(x::I) where {I<:Signed}
    (signbit(x) || x>typemax(UInt32)) && throw(DomainError("$x"))
    z = x%UInt32
    return NNInt{UInt32}(z)
end

function NNInt32(x::I) where {I<:Unsigned}
    x>typemax(UInt32) && throw(DomainError("$x"))
    z = x%UInt32
    return NNInt{UInt32}(z)
end

function NNInt64(x::I) where {I<:Signed}
    (signbit(x) || x>typemax(UInt64)) && throw(DomainError("$x"))
    z = x%UInt64
    return NNInt{UInt64}(z)
end

function NNInt64(x::I) where {I<:Unsigned}
    x>typemax(UInt64) && throw(DomainError("$x"))
    z = x%UInt64
    return NNInt{UInt64}(z)
end

+(x::NNInt32, y::NNInt32) = NNInt32(x.value + y.value)
-(x::NNInt32, y::NNInt32) = NNInt32(x.value - y.value)
*(x::NNInt32, y::NNInt32) = NNInt32(x.value * y.value)
÷(x::NNInt32, y::NNInt32) = NNInt32(x.value ÷ y.value)

fld(x::NNInt32, y::NNInt32) = NNInt32(fld(x.value, y.value))
cld(x::NNInt32, y::NNInt32) = NNInt32(cld(x.value, y.value))
rem(x::NNInt32, y::NNInt32) = NNInt32(rem(x.value, y.value))
mod(x::NNInt32, y::NNInt32) = NNInt32(mod(x.value, y.value))

end # NonNegInts


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

using Markdown
using InteractiveUtils

# ╔═╡ a343426c-952d-11ec-0693-3dbb707e1e76
begin
    import Pkg
    Pkg.activate("..")

    import GMT

    using DelimitedFiles
    using STMMRunner
    using STMMRunner.MSTM.V4
end

# ╔═╡ 1b80ded1-6068-4370-9aa7-3f04ab013b11
begin
    spheres = [
		Sphere(
			r = 0.7 + 0.5 * rand(),
			x = i,
			y = j,
			z = k,
			m = 1.5 + 1.0im * rand()
		) for i = -6.0:3.0:6.0 for j = -6.0:3.0:6.0 for k = -6.0:3.0:6.0
	]

   param = STMMConfig(
		number_processors = Sys.CPU_THREADS ÷ 2,
		sm_number_processors = Sys.CPU_THREADS ÷ 2,
		run_print_file = "",
		working_directory = "",
		orientation = FixedOrientation,
		α = 70.0,
		calculate_near_field = true,
		near_field_x_min = 0.0,
		near_field_y_min = -10.0,
	    near_field_z_min = -10.0,
	    near_field_x_max = 0.0,
	    near_field_y_max = 10.0,
	    near_field_z_max = 10.0,
	    near_field_step_size = 0.1,
		spheres = spheres,
	)
end

# ╔═╡ ffcb61ce-843f-4108-9942-bc00ae3ca36a
results = run_mstm(param; keep = false)

# ╔═╡ 7ef49694-b211-4f66-9c29-6b9704c775fc
let
    nf = results.near_field.field
    y = nf.y |> Set |> collect |> sort
    z = nf.z |> Set |> collect |> sort
    sp = results.near_field.spheres

    grd = GMT.xyz2grd(
        [nf.y nf.z real.(nf.Ex₌)],
        region = (first(y), last(y), first(z), last(z)),
        inc = "$(y[2]-y[1])/$(z[2]-z[1])",
    )
    GMT.grdimage(grd; xaxis = "af+l@[y@[", yaxis = "af+l@[z@[")
    GMT.plot!([sp.y sp.z sp.r * 14 * 2 / (last(nf.y) - first(nf.y))]; symbol = "c")
    GMT.colorbar!(;
        xaxis = "af",
        yaxis = "+l@[\\mathbf{Re}E_{{\\parallel}x}@[",
        nolines = true,
        show = true,
        savefig = "E∥x_re_v4.png",
    )
end

# ╔═╡ 224530e6-1fee-44b0-a7d8-be879224b361
let
    spheres = param.spheres
    x = map(s -> s.x, spheres)
    y = map(s -> s.y, spheres)
    z = map(s -> s.z, spheres)
    r = map(s -> s.r, spheres)
    GMT.plot3d(
        [x y z fill(-20, length(x)) r 2r];
        symbol = "e",
        show = true,
        savefig = "setting_v4.png",
    )
end

# ╔═╡ Cell order:
# ╠═a343426c-952d-11ec-0693-3dbb707e1e76
# ╠═1b80ded1-6068-4370-9aa7-3f04ab013b11
# ╠═ffcb61ce-843f-4108-9942-bc00ae3ca36a
# ╠═7ef49694-b211-4f66-9c29-6b9704c775fc
# ╠═224530e6-1fee-44b0-a7d8-be879224b361


text:
export UKSettlementCalendar, UKLSECalendar, UKLMECalendar

abstract UnitedKingdomCalendar <: WesternCalendar

type UKSettlementCalendar <: UnitedKingdomCalendar; end
type UKLSECalendar <: UnitedKingdomCalendar; end
type UKLMECalendar <: UnitedKingdomCalendar; end

string(::UKSettlementCalendar) = "UK Generic Settlement"
string(::UKLSECalendar) = "London Stock Exchange"
string(::UKLMECalendar) = "London Metal Exchange"

function isWorkingDay(c::UKSettlementCalendar, dt::CalendarTime)
	d=day(dt)
	dd=dayofyear(dt)
	m=month(dt)
	w=dayofweek(dt)
	y=year(dt)
	em=easterMonday(c,y)

	!(isWeekend(c,w)
            # New Year's Day (possibly moved to Monday)
            || ((d == 1 || ((d == 2 || d == 3) && w == Monday)) && m == January)
            # Good Friday
            || (dd == em - 3)
            # Easter MONDAY
            || (dd == em)
            # first MONDAY of May (Early May Bank Holiday)
            || (d <= 7 && w == Monday && m == May)
            # last MONDAY of MAY (Spring Bank Holiday)
            || (d >= 25 && w == Monday && m == May && y != 2002)
            # last MONDAY of August (Summer Bank Holiday)
            || (d >= 25 && w == Monday && m == August)
            # Christmas (possibly moved to MONDAY or Tuesday)
            || ((d == 25 || (d == 27 && (w == Monday || w == Tuesday))) && m == December)
            # Boxing Day (possibly moved to MONDAY or TUESDAY)
            || ((d == 26 || (d == 28 && (w == Monday || w == Tuesday))) && m == December)
            # June 3rd, 2002 only (Golden Jubilee Bank Holiday)
            # June 4rd, 2002 only (special Spring Bank Holiday)
            || ((d == 3 || d == 4) && m == June && y == 2002)
            # June, 5th, 2012 only (Queens Diamond Jubilee)
            || (d == 5 && m == June && y == 2012)
            # DECEMBER 31st, 1999 only
            || (d == 31 && m == December && y == 1999))
end

function isWorkingDay(c::UKLSECalendar, dt::CalendarTime)
	d=day(dt);
	dd=dayofyear(dt)
	m=month(dt);
	w=dayofweek(dt);
	y=year(dt);
	em=easterMonday(c,y)

	!(isWeekend(c,w)
            # New Year's Day (possibly moved to Monday)
            || ((d == 1 || ((d == 2 || d == 3) && w == Monday)) && m == January)
            # Good Friday
            || (dd == em-3)
            # Easter Monday
            || (dd == em)
            # first Monday of May (Early May Bank Holiday)
            || (d <= 7 && w == Monday && m == May)
            # last Monday of May (Spring Bank Holiday)
            || (d >= 25 && w == Monday && m == May && y != 2002)
            # last Monday of August (Summer Bank Holiday)
            || (d >= 25 && w == Monday && m == August)
            # Christmas (possibly moved to Monday or Tuesday)
            || ((d == 25 || (d == 27 && (w == Monday || w == Tuesday))) && m == December)
            # Boxing Day (possibly moved to Monday or Tuesday)
            || ((d == 26 || (d == 28 && (w == Monday || w == Tuesday)))
                && m == December)
            # June 3rd, 2002 only (Golden Jubilee Bank Holiday)
            # June 4rd, 2002 only (special Spring Bank Holiday)
            || ((d == 3 || d == 4) && m == June && y == 2002)
            # June, 5th, 2012 only (Queens Diamond Jubilee)
            || (d == 5 && m == June && y == 2012)
            # December 31st, 1999 only
            || (d == 31 && m == December && y == 1999)
            )
end

function isWorkingDay(c::UKLMECalendar, dt::CalendarTime)
	d=day(dt);
	dd=dayofyear(dt)
	m=month(dt);
	w=dayofweek(dt);
	y=year(dt);
	em=easterMonday(c,y)

	!(isWeekend(w)
            # New Year's Day (possibly moved to Monday)
            || ((d == 1 || ((d == 2 || d == 3) && w == Monday)) &&  m == January)
            # Good Friday
            || (dd == em-3)
            # Easter Monday
            || (dd == em)
            # first Monday of May (Early May Bank Holiday)
            || (d <= 7 && w == Monday && m == May)
            # last Monday of May (Spring Bank Holiday)
            || (d >= 25 && w == Monday && m == May && y != 2002)
            # last Monday of August (Summer Bank Holiday)
            || (d >= 25 && w == Monday && m == August)
            # Christmas (possibly moved to Monday or Tuesday)
            || ((d == 25 || (d == 27 && (w == Monday || w == Tuesday)))
                && m == December)
            # Boxing Day (possibly moved to Monday or Tuesday)
            || ((d == 26 || (d == 28 && (w == Monday || w == Tuesday))) && m == December)
            # June 3rd, 2002 only (Golden Jubilee Bank Holiday)
            # June 4rd, 2002 only (special Spring Bank Holiday)
            || ((d == 3 || d == 4) && m == June && y == 2002)
            # June, 5th, 2012 only (Queens Diamond Jubilee)
            || (d == 5 && m == June && y == 2012)
            # December 31st, 1999 only
            || (d == 31 && m == December && y == 1999)
            )
end


text:
using ModelParameters

function c_dam_level(volume)
    return 156.8 + 0.9463 * volume^0.2922
end

function c_dam_area(volume)
    return 0.0021 * volume^0.762
end

function c_dam_discharge(volume, max_storage)
    discharge = 0.0
    if volume > max_storage
        discharge = 0.001492 * (volume - max_storage)^1.5280
    end

    return max(0.0, discharge)
end

function c_dam_outflow(discharge, irrigation_extraction)
    return discharge + irrigation_extraction
end

Base.@kwdef mutable struct DamNode{P} <: NetworkNode
    @network_node

    max_storage::Float64
    storage_coef::P = Param(0.5, bounds=(0.00001, 10.0))

    # Function to use to calculate dam level from storage volume
    calc_dam_level::Function = c_dam_level

    # Function to use to calculate dam surface area from storage volume
    calc_dam_area::Function = c_dam_area

    # Function to use to calculate dam discharge from storage volume
    calc_dam_discharge::Function = c_dam_discharge

    # Function to calculate outflow from dam
    calc_dam_outflow::Function = c_dam_outflow

    storage::Array{Float64} = [0.0]

    effective_rainfall::Array{Float64} = []
    et::Array{Float64} = []
    inflow::Array{Float64} = []
    dam_area::Array{Float64} = []

    level::Array{Float64} = []
    discharge::Array{Float64} = []
    outflow::Array{Float64} = []

end

function DamNode(
    name::String,
    area::Float64,
    max_storage::Float64,
    storage_coef::Float64,
    initial_storage::Float64,
    calc_dam_level::Function,
    calc_dam_area::Function,
    calc_dam_discharge::Function,
    calc_dam_outflow::Function)
    return DamNode(name, area, max_storage, storage_coef,
                   calc_dam_level, calc_dam_area, calc_dam_discharge, calc_dam_outflow,
                   [initial_storage], [], [], [], [], [], [], [])
end

"""
    DamNode(name::String, spec::Dict)

Create DamNode from a given specification.
"""
function DamNode(name::String, spec::Dict)
    n = DamNode{Param}(; name=name, area=spec["area"],
                       max_storage=spec["max_storage"])

    node_params = spec["parameters"]
    for (k, p) in node_params
        s = Symbol(k)
        if p isa String
            p = eval(Meta.parse(p))
        end

        try
            if k == "initial_storage"
                setfield!(n, :storage, [p])
            else
                f = getfield(n, s)
                setfield!(n, s, Param(p, bounds=f.bounds))
            end
        catch err
            setfield!(n, s, p)
        end
    end

    return n
end

function level(node::DamNode)
    return last(node.level)
end

function area(node::DamNode)
    return last(node.dam_area)
end

function discharge(node::DamNode)
    return last(node.discharge)
end

function storage(node::DamNode)
    return last(node.storage)
end

function update_state(node::DamNode, storage, rainfall, et, area, discharge, outflow)
    push!(node.storage, storage)
    push!(node.effective_rainfall, rainfall)
    push!(node.et, et)
    push!(node.level, node.calc_dam_level(storage))
    push!(node.dam_area, area)
    push!(node.discharge, discharge)
    push!(node.outflow, outflow)

    return nothing
end

"""
    update_volume(volume, node_inflow, gamma, rain, evap, area, extractions, discharge, max_store)::Float64

Update dam volume for timestep.

# Arguments
- volume : current water volume in ML
- node_inflow : inflow from previous node in ML
- gamma : groundwater exchange (positive is gain from gw flow, negative is loss to infiltration)
- rain : rainfall input
- evap : evaporation loss
- infiltration : infiltration loss
- area : dam surface area in square kilometers
- extractions : water extraction from dam in ML
- discharge : discharge from dam in ML
- max_store : maximum dam storage in ML

# Returns
volume of water stored in dam
"""
function update_volume(volume, node_inflow, gamma, rain, evap, area, extractions, discharge, max_store)::Float64

    vol = volume + (node_inflow + gamma) + (rain - evap) * area - extractions - discharge

    return max(0.0, min(max_store, vol))
end

function run_node!(node::DamNode, climate::Climate;
                   inflow=nothing, extraction=nothing, exchange=nothing)
    timesteps = sim_length(climate)
    for ts in 1:timesteps
        run_node!(node, climate, ts;
                  inflow=inflow, extraction=extraction, exchange=exchange)
    end
end

"""
    run_node!(node::DamNode, climate::Climate, timestep::Int;
              inflow=nothing, extraction=nothing, exchange=nothing)

Run a specific node for a specified time step.

# Arguments
- `node::DamNode` :
- `climate::Climate` :
- `timestep::Int` : current time step
- `inflow::DataFrame` : Time series of inflows from any upstream node.
- `extraction::DataFrame` : Time series of water orders (expects column of `_releases`)
- `exchange::DataFrame` : Time series of groundwater flux
"""
function run_node!(node::DamNode, climate::Climate, timestep::Int;
                   inflow=nothing, extraction=nothing, exchange=nothing)
    ts = timestep
    if checkbounds(Bool, node.outflow, ts)
        if node.outflow[ts] != undef
            # already ran for this time step so no need to run
            return node.outflow[ts], node.level[ts]
        end
    end

    node_name = node.name
    rain, et = climate_values(node, climate, ts)
    wo = timestep_value(ts, node_name, "releases", extraction)
    ex = timestep_value(ts, node_name, "exchange", exchange)
    in_flow = timestep_value(ts, node_name, "inflow", inflow)
    vol = node.storage[ts]

    return run_node!(node, rain, et, vol, in_flow, wo, ex)
end

"""
Calculate outflow for the dam node for a single time step.

# Parameters
- node : DamNode
- rain : rainfall in mm
- et : evapotranspiration data in mm
- irrig_ext : irrigation extractions
- extractions : extraction data in ML
- gw_flux : groundwater interaction

# Returns
- outflow from dam
"""
function run_node!(node::DamNode,
                   rain::Float64,
                   et::Float64,
                   vol::Float64,
                   inflow::Float64,
                   extractions::Float64,
                   gw_flux::Float64)
    volume = vol
    dam_area = node.calc_dam_area(volume)
    discharge = node.calc_dam_discharge(volume, node.max_storage)

    updated_store = update_volume(volume, inflow, gw_flux, rain, et,
                                  dam_area, extractions, discharge, node.max_storage)
    outflow = node.calc_dam_outflow(discharge, extractions)

    update_state(node, updated_store, rain, et, dam_area, discharge, outflow)

    return outflow, level(node)
end

function reset!(node::DamNode)::Nothing
    node.storage = [node.storage[1]]
    node.effective_rainfall = []
    node.et = []
    node.level = []
    node.dam_area = []

    node.discharge = []
    node.outflow = []

    return nothing
end

"""
    update_params!(node::DamNode, storage_coef::Float64)::Nothing
"""
function update_params!(node::DamNode, storage_coef::Float64)::Nothing
    node.storage_coef = Param(storage_coef, bounds=node.storage_coef.bounds)
    return nothing
end