text:
#
# Groups
#

"""
    Group(f::JLDFile, name::AbstractString)

Construct an empty group named `name` at the top level of `JLDFile` `f`.
"""
Group(f::JLDFile, name::AbstractString) = Group(f.root_group, name)

"""
    Group(g::Group, name::AbstractString)

Construct a group named `name` as a child of group `g`.
"""
Group(g::Group{T}, name::AbstractString) where {T} = (g[name] = Group{T}(g.f))

"""
    lookup_offset(g::Group, name::AbstractString) -> RelOffset

Lookup the offset of a dataset in a group. Returns `UNDEFINED_ADDRESS` if the dataset is
not present. Does not inspect `unwritten_child_groups`.
"""
function lookup_offset(g::Group, name::AbstractString)
    if g.last_chunk_start_offset != -1
        # Has been saved to file, so written_links exists
        roffset = get(g.written_links, name, UNDEFINED_ADDRESS)
        roffset != UNDEFINED_ADDRESS && return roffset
    end
    roffset = get(g.unwritten_links, name, UNDEFINED_ADDRESS)
end

"""
    pathize(g::Group, name::AbstractString, create::Bool) -> Tuple{Group,String}

Converts a path to a group and name object. If `create` is true, any intermediate groups
will be created, and the dataset name will be checked for uniqueness with existing names.
"""
function pathize(g::Group, name::AbstractString, create::Bool)
    if '/' in name
        f = g.f
        dirs = split(name, '/')

        # Handles the absolute path case where the name starts with /
        if isempty(first(dirs))
            g = g.f.root_group
            start = 2
        else
            start = 1
        end

        for i = start:length(dirs)-1
            dir = dirs[i]
            isempty(dir) && continue

            # See if a group already exists
            offset = lookup_offset(g, dir)
            if offset == UNDEFINED_ADDRESS
                if haskey(g.unwritten_child_groups, dir)
                    # It's possible that lookup_offset fails because the group has not yet
                    # been written to the file
                    g = g.unwritten_child_groups[dir]
                elseif create
                    # No group exists, so create a new group
                    newg = Group{typeof(f)}(f)
                    g.unwritten_child_groups[dir] = newg
                    g = newg
                else
                    throw(KeyError(join(dirs[1:i], '/')))
                end
            elseif haskey(f.loaded_groups, offset)
                g = f.loaded_groups[offset]
            elseif !isgroup(f, offset)
                throw(ArgumentError("path $(join(dirs[1:i], '/')) points to a dataset, not a group"))
            else
                g = f.loaded_groups[offset] = load_group(f, offset)
            end
        end

        name = String(dirs[end])
    end

    if create && haskey(g, name)
        throw(ArgumentError("a group or dataset named $name is already present within this group"))
    end

    return (g, name)
end

function Base.getindex(g::Group, name::AbstractString)
    f = g.f
    f.n_times_opened == 0 && throw(ArgumentError("file is closed"))

    (g, name) = pathize(g, name, false)

    roffset = lookup_offset(g, name)
    if roffset == UNDEFINED_ADDRESS
        haskey(g.unwritten_child_groups, name) && return g.unwritten_child_groups[name]
        throw(KeyError(name))
    end

    if isgroup(f, roffset)
        let loaded_groups = f.loaded_groups
            get!(()->load_group(f, roffset), loaded_groups, roffset)
        end
    else
        load_dataset(f, roffset)
    end
end

function Base.write(g::Group, name::AbstractString, obj, wsession::JLDWriteSession=JLDWriteSession())
    if g.last_chunk_start_offset != -1 && g.continuation_message_goes_here == -1
        error("objects cannot be added to this group because it was created with a previous version of JLD2")
    end
    f = g.f
    prewrite(f)
    (g, name) = pathize(g, name, true)
    g[name] = write_dataset(f, obj, wsession)
    nothing
end

function Base.setindex!(g::Group, obj, name::AbstractString)
    write(g, name, obj)
    g
end

function Base.setindex!(g::Group, child::Group, name::AbstractString)
    if child.last_chunk_start_offset != -1
        throw(ArgumentError("cannot re-link a group that has already been written"))
    end
    prewrite(g.f)
    (g, name) = pathize(g, name, true)
    g.unwritten_child_groups[name] = child
    g
end

# Internal use only
function Base.setindex!(g::Group, offset::RelOffset, name::AbstractString)
    if g.last_chunk_start_offset != -1 && g.continuation_message_goes_here == -1
        error("objects cannot be added to this group because it was created with a previous version of JLD2")
    end
    g.unwritten_links[name] = offset
    g
end

function Base.haskey(g::Group, name::AbstractString)
    (g, name) = pathize(g, name, false)
    (g.last_chunk_start_offset != -1 && haskey(g.written_links, name)) ||
        haskey(g.unwritten_links, name) || haskey(g.unwritten_child_groups, name)
end

Base.isempty(g::Group) =
    (g.last_chunk_start_offset == -1 || isempty(g.written_links)) &&
    isempty(g.unwritten_links) && isempty(g.unwritten_child_groups)

function Base.keys(g::Group)
    ks = String[]
    if g.last_chunk_start_offset != -1
        append!(ks, keys(g.written_links))
    end
    append!(ks, keys(g.unwritten_links))
    append!(ks, keys(g.unwritten_child_groups))
    ks
end

struct LinkInfo
    version::UInt8
    flags::UInt8
    fractal_heap_address::RelOffset
    name_index_btree::RelOffset
end
define_packed(LinkInfo)

LinkInfo() = LinkInfo(0, 0, UNDEFINED_ADDRESS, UNDEFINED_ADDRESS)

@enum(CharacterSet,
      CSET_ASCII,
      CSET_UTF8)

const LM_CREATION_ORDER_PRESENT = UInt8(2^2)
const LM_LINK_TYPE_FIELD_PRESENT = UInt8(2^3)
const LM_LINK_NAME_CHARACTER_SET_FIELD_PRESENT = UInt8(2^4)

function read_link(io::IO)
    # Version
    version = read(io, UInt8)
    version == 1 || throw(UnsupportedVersionException())

    # Flags
    flags = read(io, UInt8)

    if (flags & LM_LINK_TYPE_FIELD_PRESENT) != 0
        read(io, UInt8) == 0 || throw(UnsupportedFeatureException())
    end

    if (flags & LM_CREATION_ORDER_PRESENT) != 0
        skip(io, 8)
    end

    # Link name character set
    cset = CSET_ASCII
    if (flags & LM_LINK_NAME_CHARACTER_SET_FIELD_PRESENT) != 0
        cset_byte = read(io, UInt8)
        cset = CharacterSet(cset_byte)
    end

    sz = read_size(io, flags)  # Size
    name = read(io, UInt8, sz) # Link name
    target = read(io, RelOffset)  # Link information

    (String(name), target)
end

function load_group(f::JLDFile, roffset::RelOffset)
    io = f.io
    chunk_start_offset::Int64 = fileoffset(f, roffset)
    seek(io, chunk_start_offset)

    cio = begin_checksum_read(io)
    sz = read_obj_start(cio)
    chunk_checksum_offset::Int64 = position(cio) + sz

    # Messages
    continuation_message_goes_here::Int64 = -1
    links = OrderedDict{String,RelOffset}()

    continuation_offset::Int64 = -1
    continuation_length::Length = 0

    while true
        if continuation_offset != -1
            seek(io, continuation_offset)
            chunk_checksum_offset = continuation_offset + continuation_length - 4
            continuation_offset = -1

            cio = begin_checksum_read(io)
            read(cio, UInt32) == OBJECT_HEADER_CONTINUATION_SIGNATURE || throw(InvalidDataException())
        end

        while (curpos = position(cio)) <= chunk_checksum_offset - 4
            msg = read(cio, HeaderMessage)
            endpos = curpos + sizeof(HeaderMessage) + msg.size
            if msg.msg_type == HM_NIL
                if continuation_message_goes_here == -1 &&
                   chunk_checksum_offset - curpos >= sizeof(HeaderMessage) + sizeof(RelOffset) + sizeof(Length)
                    continuation_message_goes_here = curpos
                end
            else
                continuation_message_goes_here = -1
                if msg.msg_type == HM_LINK_INFO
                    link_info = read(cio, LinkInfo)
                    link_info.fractal_heap_address == UNDEFINED_ADDRESS || throw(UnsupportedFeatureException())
                elseif msg.msg_type == HM_GROUP_INFO
                    # This message doesn't help us much, so we ignore it for now
                elseif msg.msg_type == HM_LINK_MESSAGE
                    name, loffset = read_link(cio)
                    links[name] = loffset
                elseif msg.msg_type == HM_OBJECT_HEADER_CONTINUATION
                    continuation_offset = chunk_start_offset = fileoffset(f, read(cio, RelOffset))
                    continuation_length = read(cio, Length)
                elseif (msg.flags & 2^3) != 0
                    throw(UnsupportedFeatureException())
                end
            end
            seek(cio, endpos)
        end

        # Checksum
        seek(cio, chunk_checksum_offset)
        end_checksum(cio) == read(io, UInt32) || throw(InvalidDataException())

        continuation_offset == -1 && break
    end

    Group{typeof(f)}(f, chunk_start_offset, continuation_message_goes_here,
                     chunk_checksum_offset, OrderedDict{String,RelOffset}(),
                     OrderedDict{String,Group}(), links)
end

"""
    link_size(name::String)

Returns the size of a link message, including message header.
"""
link_size(name::String) =
    3 + size_size(sizeof(name)) + sizeof(name) + sizeof(RelOffset)

"""
    links_size(pairs)

Returns the size of several link messages. `pairs` is an iterator of
`String => RelOffset` pairs.
"""
function links_size(pairs)
    sz = 0
    for (name::String,) in pairs
        sz += link_size(name) + sizeof(HeaderMessage)
    end
    sz
end

"""
    group_payload_size(pairs)

Returns the size of a group payload, including link info, group info, and link messages,
but not the object header. `pairs` is an iterator of `String => RelOffset` pairs. Provides
space after the last object message for a continuation message.
"""
group_payload_size(pairs) =
    sizeof(HeaderMessage) + sizeof(LinkInfo) + sizeof(HeaderMessage) + 2 + links_size(pairs) +
    sizeof(HeaderMessage) + sizeof(RelOffset) + sizeof(Length)

group_continuation_size(pairs) =
    sizeof(OBJECT_HEADER_CONTINUATION_SIGNATURE) + links_size(pairs) +
    sizeof(HeaderMessage) + sizeof(RelOffset) + sizeof(Length) + 4 # Checksum is included

"""
    save_group(g::Group) -> RelOffset

Stores a group to a file, updating it if it has already been saved. Returns
UNDEFINED_ADDRESS if the group was already stored, or the offset of the new group
otherwise.
"""
function save_group(g::Group)
    # Save unwritten child groups
    if !isempty(g.unwritten_child_groups)
        for (name, child_group) in g.unwritten_child_groups
            retval = save_group(child_group)
            g.unwritten_links[name] = retval
        end
        empty!(g.unwritten_child_groups)
    end

    # Save links
    f = g.f
    io = f.io
    if g.last_chunk_start_offset == -1
        psz = group_payload_size(g.unwritten_links)
        sz = sizeof(ObjectStart) + size_size(psz) + psz

        g.last_chunk_start_offset = f.end_of_data
        retval = h5offset(f, f.end_of_data)
        seek(io, f.end_of_data)
        cio = begin_checksum_write(io, sz)

        # Object header
        write(cio, ObjectStart(size_flag(psz)))
        write_size(cio, psz)
        x = position(cio)

        # Link info message
        write(cio, HeaderMessage(HM_LINK_INFO, sizeof(LinkInfo), 0))
        write(cio, LinkInfo())

        # Group info message
        write(cio, HeaderMessage(HM_GROUP_INFO, 2, 0))
        write(cio, UInt16(0))
    else
        # If no changes, no need to save
        isempty(g.unwritten_links) && return UNDEFINED_ADDRESS

        retval = UNDEFINED_ADDRESS
        continuation_start = f.end_of_data
        csz = group_continuation_size(g.unwritten_links)

        # Object continuation message
        seek(io, g.continuation_message_goes_here)
        write(io, HeaderMessage(HM_OBJECT_HEADER_CONTINUATION, sizeof(RelOffset) + sizeof(Length), 0))
        write(io, h5offset(f, continuation_start))
        write(io, Length(csz))

        # Re-calculate checksum
        seek(io, g.last_chunk_start_offset)
        cio = begin_checksum_read(io)
        seek(cio, g.last_chunk_checksum_offset)
        seek(io, g.last_chunk_checksum_offset)
        write(io, end_checksum(cio))

        # Object continuation
        seek(io, continuation_start)
        g.last_chunk_start_offset = continuation_start
        cio = begin_checksum_write(io, csz - 4)
        write(cio, OBJECT_HEADER_CONTINUATION_SIGNATURE)
    end

    # Links
    for (name, offset) in g.unwritten_links
        write(cio, HeaderMessage(HM_LINK_MESSAGE, link_size(name), 0))
        write(cio, UInt8(1))             # Version

        # Flags
        flags = size_flag(sizeof(name)) | LM_LINK_NAME_CHARACTER_SET_FIELD_PRESENT
        write(cio, flags::UInt8)

        write(cio, UInt8(CSET_UTF8))     # Link name character set
        write_size(cio, sizeof(name))    # Length of link name
        write(cio, name)                 # Link name
        write(cio, offset)               # Link target
    end
    empty!(g.unwritten_links)

    # Extra space for object continuation
    g.continuation_message_goes_here = position(cio)
    write(cio, HeaderMessage(HM_NIL, sizeof(RelOffset)+sizeof(Length), 0))
    write(cio, RelOffset(0))
    write(cio, Length(0))

    # Checksum
    write(io, end_checksum(cio))
    f.end_of_data = position(io)
    g.last_chunk_checksum_offset = f.end_of_data - 4

    return retval
end

function show_group(io::IO, g::Group, prefix::String=" ", skiptypes::Bool=false)
    if g.f.n_times_opened == 0
        print(io, "  (closed)")
        return
    end

    ks = collect(keys(g))
    skiptypes && filter!(x -> x != "_types", ks)

    if isempty(ks) && prefix == " "
        print(io, "  (no datasets)")
        return
    end

    for i = 1:length(ks)
        k = ks[i]
        islast = i == length(ks)
        isagroup = haskey(g.unwritten_child_groups, k) || isgroup(g.f, lookup_offset(g, k))
        print(io, prefix, islast ? "└─" : "├─", isagroup ? "📂 " : "🔢 ", k)
        if isagroup
            newg = g[k]
            if !isempty(newg)
                print(io, '\n')
                show_group(io, newg, prefix*(islast ? "   " : "│  "))
            end
        end
        !islast && print(io, '\n')
    end
end

function Base.show(io::IO, g::Group)
    println(io, "JLD2.Group")
    show_group(io, g)
end


text:
function simulate_LIF_neuron(input_current;
                                simulation_time = 5,
                                dt = dt,
                                v_rest=-70, # * b2.mV
                                v_reset=-65, # * b2.mV,
                                firing_threshold=-50, # * b2.mV,
                                membrane_resistance=10., # * b2.Mohm,
                                membrane_time_scale=8., # * b2.ms,
                                abs_refractory_period=2.0, # * b2.ms,
                                doEuler = false)

    params = LIFParams(; tau = membrane_time_scale, reset = v_reset,
                            threshold = firing_threshold, delay = dt,
                            refractory_period = abs_refractory_period)
    network_init = FeedForwardNet([1], params;
                                recorder = AllSpikesRecorder(),
                                plasticityrule = NoPlasticity)

    doEuler && (network_init = EulerNet(network_init; Δt = dt))
    network_init.v[1] = v_rest
    network = deepcopy(network_init)

    input = input_current + network_init.parameters.reset

    b = @benchmark integratenet!(x[1], x[2], x[3]) setup = (x = (deepcopy($network_init), copy($input), copy($simulation_time))) evals = 1;
    profiling_time = BenchmarkTools.median(b).time

    integratenet!(network, input, simulation_time)
    rec = deepcopy(network.recorder)
    return network, rec, profiling_time / 1e9
end

function simulate_balanced_network(input_current;
                                n_of_neurons = 10^3,
                                J = 0.2,
                                g = 5.0,
                                f = 3/4,
                                simulation_time = 5,
                                dt = dt,
                                v_rest=-70, # * b2.mV
                                v_reset=-65, # * b2.mV,
                                firing_threshold=-50, # * b2.mV,
                                membrane_resistance=10., # * b2.Mohm,
                                membrane_time_scale=8., # * b2.ms,
                                abs_refractory_period=2.0, # * b2.ms,
                                doEuler = false)
    params = LIFParams(; tau = membrane_time_scale, reset = v_reset,
                            threshold = firing_threshold, delay = dt,
                            refractory_period = abs_refractory_period)
    network_init = BalancedNet(n_of_neurons, J, g; f = f, parameters = params)

    doEuler && (network_init = EulerNet(network_init; Δt = dt))
    network_init.v[:] .= v_rest
    network = deepcopy(network_init)

    input = input_current + network_init.parameters.reset

    b = @benchmark integratenet!(x[1], x[2], x[3]) setup = (x = (deepcopy($network_init), copy($input), copy($simulation_time))) evals = 1;
    profiling_time = BenchmarkTools.median(b).time

    integratenet!(network, input, simulation_time)
    rec = deepcopy(network.recorder)
    return network, rec, profiling_time / 1e9
end

ext: jl


text:
# SPDX-License-Identifier: X11
# 2020-10-18 (Timed)
# Serval vs Monster (100pt)

H, A = parse.(Int, split(readline()))
println(cld(H, A))

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
#=
Pr(C(x) = i) = w_i / Σ w_i

The derivation is easy, see [wrh] "Proof of Correctness".
=#

export WrhCoordinator, insert_node!, remove_node!, choose_nodes, num_nodes, node_weight, reseed!

# Make it mutable for reseeding
mutable struct WrhNode
    id::String
    seed::UInt64
    weight::Float64
end

"Weighted Rendezvous Hashing"
struct WrhCoordinator <: RvCoordinator
    nodes::Dict{String, WrhNode}

    WrhCoordinator() = new(Dict{String, WrhNode}())
end

function insert_node!(coordinator::WrhCoordinator, id::String, weight::Float64)
    node = WrhNode(id, hash(id), weight)
    coordinator.nodes[id] = node
end

function remove_node!(coordinator::WrhCoordinator, id::String)
    delete!(coordinator.nodes, id)
end

function choose_nodes(
    coordinator::WrhCoordinator,
    data::String,
    n::Union{Int, Nothing}=nothing,
)::Vector{String}
    nodes = [
        begin
            h = hash(data, node.seed)
            h_float = (h & MASK) / MAX_HASH
            s = -node.weight / log(h_float)
            s => id
        end
        for (id, node) in coordinator.nodes
    ]

    n::Int = isnothing(n) ? length(nodes) : n
    partialsort!(nodes, 1:n;
        rev = true  # descending
    ) .|> last
end

function num_nodes(coordinator::WrhCoordinator)::Int
    length(coordinator.nodes)
end

function node_weight(coordinator::WrhCoordinator, id::String)::Float64
    coordinator.nodes[id].weight
end

function reseed!(coordinator::WrhCoordinator)
    for node in values(coordinator.nodes)
        node.seed = rand(UInt64)
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Genie
import Genie.Router: route
import Genie.Renderer: json

Genie.config.run_as_server = true
Genie.config.log_level = :error
Genie.config.app_env = "prod"
Genie.config.server_host = "0.0.0.0"

route("/") do
  (:data => "Hello World!") |> json
end

Genie.startup(4500)

max_stars_repo_name: kislerdm/web-server-benchmark

max_issues_repo_name: kislerdm/web-server-benchmark

max_forks_repo_name: kislerdm/web-server-benchmark


text:
module distances

export distance

using StringDistances

const metrics = Dict(:jaro => Jaro, :dameraulevenshtein => DamerauLevenshtein, :levenshtein => Levenshtein, :ratcliffobershelp => RatcliffObershelp)

"""
Returns: Distance between val1 and val2 as measured by the metric.

Distances are scaled to be in [0, 1], with 0 indicating that the values are equal and 1 indicating that the values are totally different.
"""
function distance(metric::Symbol, val1, val2)
    ismissing(val1) && return 1.0
    ismissing(val2) && return 1.0
    1.0 - compare(val1, val2, metrics[metric]())
end

end


text:
"""
    pet_to_stp(pet_table, detector_SSD)

    Table, SolidStateDetectors.Simulation -> Table

Construct a table with "stepping info" based on
    position-energy-time information of simulated energy depositions
    given in <pet_table> and detector geometry given in <detector_SSD>.

Current steps include:
    - removing events outside of detector PV
    - clustering
    - grouping by detector
    - removing events reconstructed to be outisde of the detector
"""
function pet_to_stp(pet_table::Table, detector_SSD::SolidStateDetectors.Simulation)
    ## in case it's the Geant4 CSV file from LegendTestData, it's already grouped correctly
    hits_by_evtno = pet_table

    ## in case it's a g4simple HDF5 file

    # Save only events that occur in the detector PV
    # volID = 1 selects only events in the detector
    if hasproperty(pet_table, :volID)
        @info "...checking volID"
        pet_pv = group_by_column(pet_table, :volID)[1]

        # Need to turn into a normal Table before using internal SSD functions (group_by_evtno, cluster_detector_hits, etc)
        pet_flat = Table(
            evtno = pet_pv.evtno,
            detno = pet_pv.detno,
            thit = pet_pv.thit,
            edep = pet_pv.edep,
            pos = pet_pv.pos
        )

        # group hits by event number
        hits_by_evtno = RadiationDetectorSignals.group_by_evtno(pet_flat)
    end

    # cluster hits based on distance
    @info "...clustering"
    println("$(sum(length.(hits_by_evtno.edep))) hits before clustering")
    events_clustered = @time SolidStateDetectors.cluster_detector_hits(hits_by_evtno, 0.2u"mm")
    println("$(sum(length.(events_clustered.edep))) hits after clustering")

    hits = RadiationDetectorSignals.ungroup_by_evtno(events_clustered)

    # Waveform generation has to be per detector.
    # Let's reshuffle the detector hits, grouping by event number and detector
    # @info("...grouping by detector")
    events_per_det = RadiationDetectorSignals.group_by_evtno_and_detno(hits)

    # The hits are now grouped by event number, but separately for each detector, and sorted by detector number:
    # issorted(events_per_det.detno) # what does this do? It's not "!"

    ## ---------- This part is done for the case of multiple detectors
    ## For now we consider that the input file should have only one detector
    ## TODO
    # This makes it easy to group them by detector number ...
    # events_det1 = filter(evt -> evt.detno == 0, events_per_det)
    ## ---------- This part is done for the case of multiple detectors

    @info "...removing events outside of the detector"

    # We need to filter out the few events that, due to numerical effects, lie outside of the detector
    # (the proper solution is to shift them slightly, this feature will be added in the future)
    # stp_events = events_det1[findall(pts -> all(p -> T.(ustrip.(uconvert.(u"m", p))) in detector_SSD.detector, pts), events_det1.pos)]
    stp_events = events_per_det[findall(pts -> all(p -> T.(ustrip.(uconvert.(u"m", p))) in detector_SSD.detector, pts), events_per_det.pos)]

    stp_events

end

function pet_to_stp(pet_filename::AbstractString, det_metadata::AbstractString)
    # construct simplified config based on given info
    sim_config = PropDict(:detector_metadata => det_metadata, :input_file => pet_filename)

    pet_to_stp(sim_config)
end

function pet_to_stp(sim_config_file::AbstractString)
    # given config file already contains inputs as well
    pet_to_stp(propdict(sim_config_file))
end

function pet_to_stp(sim_config::PropDict)
    @info "---------------------- pet -> stp (stepping info)"
    # SSD detector for geometry purposes (not simulation)
    det_config_SSD = ssd_config(sim_config.detector_metadata, Environment())
    detector_SSD = Simulation(SolidStateDetector{T}(det_config_SSD))

    println("Processing file: $(sim_config.input_file) for detector $(sim_config.detector_metadata)")
    pet_table = read_pet(sim_config.input_file)

    # launch pet->stp
    pet_to_stp(pet_table, detector_SSD)
end

"""
AbstractString -> Table
"""
function read_pet(pet_filename::AbstractString)
    # check file extension and determine the type
    file_ext = splitext(pet_filename)[end]
    file_type = file_ext == ".csv" ? LegendTextIO.Geant4CSVInput : LegendHDF5IO.Geant4HDF5Input

    # read_pet(pet_filename, file_type)
    open(pet_filename, file_type) do io
        read(io)
    end
end

# function read_pet(pet_filename::AbstractString, file_type::Type{Geant4CSVInput})
#     open(pet_filename, file_type) do io
#         read(io)
#     end
# end

# function read_pet(pet_filename::AbstractString, ::Type{Geant4HDF5Input})
#     g4sfile = h5open(pet_filename, "r")
#     g4sntuple = g4sfile["default_ntuples"]["g4sntuple"]

#     evtno = read(g4sntuple["event"]["pages"])
#     detno = read(g4sntuple["iRep"]["pages"])
#     thit = read(g4sntuple["t"]["pages"]).*u"ns"
#     edep = read(g4sntuple["Edep"]["pages"]).*u"MeV"
#     volID = read(g4sntuple["volID"]["pages"])

#     x = read(g4sntuple["x"]["pages"])
#     y = read(g4sntuple["y"]["pages"])
#     z = read(g4sntuple["z"]["pages"])

#     # Construct array of positions for input to SSD
#     n_ind = length(evtno)
#     pos = [ SVector{3}(([ x[i], y[i], z[i] ] .* u"mm")...) for i in 1:n_ind ]

#     # Construct a Table with the arrays we just constructed from the g4sfile data
#     TypedTables.Table(
#             evtno = evtno,
#             detno = detno,
#             thit = thit,
#             edep = edep,
#             pos = pos,
#             volID = volID,
#     )
# end

"""
    group_by_column(table, colname)

TypedTables.Table, Symbol -> TypedTables.Table

Group table <table> by given column <colname>.

I think this function already exists somewhere.
Gotta ask Oliver and maybe use that package instead.
"""
function group_by_column(table::TypedTables.Table, colname::Symbol)
    sorting_idxs = sortperm(getproperty(table, colname))
    sorted = table[sorting_idxs]
    TypedTables.Table(consgroupedview(getproperty(sorted, colname), Tables.columns(sorted)))
end


text:
mutable struct SpecificInfo <: InfoDist
    joint::Matrix{Int}
    m1::Vector{Int}
    m2::Vector{Int}
    b1::Int
    b2::Int
    N::Int
    function SpecificInfo(b1::Int, b2::Int)
        if b1 < 1 || b2 < 1
            throw(ArgumentError("the support of each random variable must be at least 2"))
        end
        new(zeros(Int, b1, b2), zeros(Int, b1), zeros(Int, b2), b1, b2, 0)
    end
end

function SpecificInfo(xs::AbstractVector{Int}, ys::AbstractVector{Int})
    if isempty(xs) || isempty(ys)
        throw(ArgumentError("arguments must not be empty"))
    end
    xmin, xmax = extrema(xs)
    ymin, ymax = extrema(ys)
    if xmin < 1 || ymin < 1
        throw(ArgumentError("observations must be positive, nonzero"))
    end
    observe!(SpecificInfo(max(2, xmax), max(2, ymax)), xs, ys)
end

function observe!(dist::SpecificInfo, xs::AbstractVector{Int}, ys::AbstractVector{Int})
    if length(xs) != length(ys)
        throw(ArgumentError("arguments must have the same length"))
    end
    dist.N += length(xs)
    for i in eachindex(xs)
        dist.joint[xs[i], ys[i]] += 1
        dist.m1[xs[i]] += 1
        dist.m2[ys[i]] += 1
    end
    dist
end

function estimate(dist::SpecificInfo)
    si = zeros(dist.b1)
    N = dist.N
    for x1 in eachindex(dist.m1)
        n1 = dist.m1[x1]
        for x2 in eachindex(dist.m2)
            if !iszero(dist.joint[x1,x2])
                j, n2 = dist.joint[x1,x2], dist.m2[x2]
                si[x1] += j * log2((N * j) / (n1 * n2))
            end
        end
        si[x1] /= n1
    end
    si
end

@inline function clear!(dist::SpecificInfo)
    dist.joint[:] .= 0
    dist.m1[:] .= 0
    dist.m2[:] .= 0
    dist.N = 0
    dist
end

function specificinfo!(si::SpecificInfo, stimulus::AbstractVector{Int},
                       responses::AbstractMatrix{Int})
    specificinfo!(si, stimulus, box(responses))
end

function specificinfo!(si::SpecificInfo, stimulus::AbstractVector{Int},
                       responses::AbstractMatrix{Int}, subset::AbstractVector{Int})
    specificinfo!(si, stimulus, @view responses[subset, :])
end

function specificinfo!(si::SpecificInfo, stimulus::AbstractVector{Int},
                       responses::AbstractVector{Int})
    estimate(observe!(dist, stimulus, responses))
end

function specificinfo(stimulus::AbstractVector{Int}, responses::AbstractMatrix{Int})
    specificinfo(stimulus, box(responses))
end

function specificinfo(stimulus::AbstractVector{Int}, responses::AbstractMatrix{Int},
                      subset::AbstractVector{Int})
    specificinfo(stimulus, @view responses[subset, :])
end

function specificinfo(stimulus::AbstractVector{Int}, responses::AbstractVector{Int})
    estimate(SpecificInfo(stimulus, responses))
end


text:
F(n) = n < 1 ? one(n) : n - M(F(n - 1))
M(n) = n < 1 ? zero(n) : n - F(M(n - 1))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Gui_Test

using Intan, SpikeSorting, Gtk.ShortNames, MAT, JLD

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

myamp=RHD2164("PortA1")
d=Debug(string(dirname(Base.source_path()),"/data/qq.mat"),"qq")
myt=Task_NoTask()
myfpga=FPGA(1,myamp)
(myrhd,ss,myfpgas)=makeRHD([myfpga],debug=d)

handles = makegui(myrhd,ss,myt,myfpgas)

sleep(1.0)

    @test handles.sc.mi == (0.0,0.0)

#=
Callback Testing
=#

#Initialization
Intan.init_cb(handles.init.handle,(handles,myrhd,myt,myfpgas))

    @test myfpgas[1].numDataStreams == 2
    @test myfpgas[1].dataStreamEnabled == [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

#Run
setproperty!(handles.run,:active,true)
sleep(1.0)
Intan.run_cb(handles.run.handle,(handles,myrhd,ss,myt,myfpgas))

    myreads=myrhd.reads
    for i=1:5
        sleep(0.5)
        @test myrhd.reads > myreads
        myreads=myrhd.reads
    end

#Calibration
sleep(1.0)

    @test myrhd.cal == 3

#Slider 1
sleep(1.0)

    for i=1:4
        setproperty!(handles.adj,:value,i)
        Intan.update_c1(handles.adj.handle,(handles,))
        sleep(1.0)
        @test handles.num16 == i
        @test handles.sc.spike == 16*i-16+handles.num
    end

#Slider 2
sleep(1.0)

    for i=1:4
        setproperty!(handles.adj2,:value,i)
        Intan.update_c2(handles)
        sleep(1.0)
        @test handles.num == i
        @test handles.sc.spike == 16*handles.num16-16+i
    end

#=
Threshold Test
=#

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.thres_widgets.show),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.thres_widgets.show,"clicked",Bool,press)

sleep(5.0)

#=
Spike Button Tests
=#

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.spike_widgets.refresh),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.spike_widgets.refresh,"clicked",Bool,press)

sleep(1.0)

#=
Gain Button Tests
=#

#=
Right canvas callbacks
=#

#16 Channel Select

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),1.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

    @test handles.sc.spike == 49

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),200.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

    @test handles.sc.spike == 53

#32 Channel Select

#Change to 32 channel
press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb1[2]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.rb1[2],"clicked",Bool,press)
sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),1.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

    @test handles.sc.spike == 33

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),150.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

    @test handles.sc.spike == 39

#=
RadioButtons
=#
for i=2:5
    press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb1[i]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
    signal_emit(handles.rb1[i],"clicked",Bool,press)
    sleep(1.0)

	@test handles.c_right_top == i

end

for i=1:6
    press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb2[i]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
    signal_emit(handles.rb2[i],"clicked",Bool,press)
    sleep(1.0)

	@test handles.c_right_bottom == i

end

#=
Soft Scope
=#

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb2[4]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.rb2[4],"clicked",Bool,press)

sleep(1.0)

for i=0:4
    Intan.scope_popup_v_cb(handles.run.handle,(handles,i))
    sleep(1.0)
end

Intan.scope_popup_v_cb(handles.run.handle,(handles,2))
sleep(1.0)

for i=0:4
    Intan.scope_popup_t_cb(handles.run.handle,(handles,i))
    sleep(1.0)
end

Intan.scope_popup_t_cb(handles.run.handle,(handles,2))
sleep(1.0)

for i=0:1
    Intan.scope_popup_thres_cb(handles.run.handle,(handles,i))
    sleep(1.0)
end

#=
Popup Tests
=#

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb1[1]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.rb1[1],"clicked",Bool,press)
sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),1.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

Intan.popup_enable_cb(handles.run.handle,(handles,myrhd))

    @test handles.enabled[handles.sc.spike] == true

Intan.popup_disable_cb(handles.run.handle,(handles,myrhd))

    @test handles.enabled[handles.sc.spike] == false

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.rb1[2]),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.rb1[2],"clicked",Bool,press)
sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.c),Int8(0),UInt32(0),1.0,1.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.c,"button-press-event",Bool,press)
sleep(1.0)

Intan.popup_enable_cb(handles.run.handle,(handles,myrhd))

    @test handles.enabled[handles.sc.spike] == true

Intan.popup_disable_cb(handles.run.handle,(handles,myrhd))

    @test handles.enabled[handles.sc.spike] == false

#=
SAVE LOAD Test
=#

sleep(1.0);

setproperty!(handles.run,:active,false)
sleep(1.0);
#=
Intan.run_cb(handles.run.handle,(handles,myrhd,ss,myt,myfpgas))

sleep(1.0);

myv=parse_v(myrhd.save.v)

    @test size(myv,2) == 64

mys_m=save_ts_mat(myrhd.save.ts)
mys_j=save_ts_jld(myrhd.save.ts)

    @test length(mys_m) == length(mys_j)

=#
#=
Close
=#

destroy(handles.win)

end


text:
module MINLP

using JuMP, BARON, Test

m = Model(BARON.Optimizer)
ub = [2, 2, 1]
@variable(m, 0 ≤ x[i=1:3] ≤ ub[i])
@variable(m, y[1:3], Bin)

@NLconstraints(m, begin
    0.8log(x[2]+1) + 0.96log(x[1]-x[2]+1) - 0.8x[3]         ≥  0
       log(x[2]+1) + 1.20log(x[1]-x[2]+1)  -   x[3] - 2y[3] ≥ -2
       x[2] ≤  x[1]
       x[2] ≤ 2y[1]
       x[1] - x[2] ≤ 2y[2]
       y[1] + y[2] ≤ 1
end)

@NLobjective(m, Min, 5y[1] + 6y[2] + 8y[3] + 10x[1] - 7x[3] - 18log(x[2]+1) -
                        19.2log(x[1]-x[2]+1) + 10)

optimize!(m)

@test isapprox(value(x[1]), 1.300975890892825, rtol=1e-6)
@test isapprox(value(x[2]), 0.0, rtol=1e-6)
@test isapprox(value(x[3]), 1.0, rtol=1e-6)
@test isapprox(value(y[1]), 0.0, rtol=1e-6)
@test isapprox(value(y[2]), 1.0, rtol=1e-6)
@test isapprox(value(y[3]), 0.0, rtol=1e-6)
@test isapprox(objective_value(m), 6.00975890893, rtol=1e-6)

end # module


text:
module Parser

export parse_domain, parse_problem, parse_pddl, @pddl, @pddl_str
export load_domain, load_problem

using ParserCombinator, Julog
using ..PDDL: Signature, GenericDomain, GenericProblem, GenericAction
using ..PDDL: DEFAULT_REQUIREMENTS, IMPLIED_REQUIREMENTS

"PDDL keyword."
struct Keyword
    name::Symbol
end
Base.show(io::IO, kw::Keyword) = print(io, "KW:", kw.name)

"Parsers for top-level PDDL descirptions."
const top_level_parsers = Dict{Symbol,Function}()

"Header field parsers for top-level PDDL descriptions (domains, problems, etc.)."
const head_field_parsers = Dict{Symbol,Dict{Symbol,Function}}(
    :domain => Dict{Symbol,Function}(),
    :problem => Dict{Symbol,Function}()
)

"Body field parsers for top-level PDDL descriptions (domains, problems, etc.)."
const body_field_parsers = Dict{Symbol,Dict{Symbol,Function}}(
    :domain => Dict{Symbol,Function}(),
    :problem => Dict{Symbol,Function}()
)

# Utility functions
include("utils.jl")
# Lisp-style grammar for PDDL
include("grammar.jl")
# Parsers for PDDL formulas
include("formulas.jl")
# Parser for top-level PDDL descriptions
include("descriptions.jl")
# Parsers for PDDL domains
include("domain.jl")
# Parsers for PDDL problems
include("problem.jl")

"Parse to PDDL structure based on initial keyword."
function parse_pddl(expr::Vector)
    if isa(expr[1], Keyword)
        kw = expr[1].name
        for desc in keys(top_level_parsers)
            if kw in keys(head_field_parsers[desc])
                return head_field_parsers[desc][kw](expr)
            elseif kw in keys(body_field_parsers[desc])
                return body_field_parsers[desc][kw](expr)
            end
        end
        error("Keyword $kw not recognized.")
    elseif expr[1] == :define
        kw = expr[2][1]
        return top_level_parsers[kw](expr)
    else
        return parse_formula(expr)
    end
end
parse_pddl(sym::Symbol) = parse_formula(sym)
parse_pddl(str::AbstractString) = parse_pddl(parse_string(str))
parse_pddl(strs::AbstractString...) = [parse_pddl(parse_string(s)) for s in strs]

"Parse string(s) to PDDL construct."
macro pddl(str::AbstractString)
    return parse_pddl(str)
end

macro pddl(strs::AbstractString...)
    return collect(parse_pddl.(strs))
end

"Parse string to PDDL construct."
macro pddl_str(str::AbstractString)
    return parse_pddl(str)
end

"Load PDDL domain from specified path."
function load_domain(path::AbstractString)
    str = open(f->read(f, String), path)
    return parse_domain(str)
end

"Load PDDL problem from specified path."
function load_problem(path::AbstractString)
    str = open(f->read(f, String), path)
    return parse_problem(str)
end

end


text:
function device_mass_matrix_entries!(::AbstractArray, ::DynamicWrapper{T}) where {T}
    error("Mass Matrix not implemented for models $T")
end

function device_mass_matrix_entries!(
    mass_matrix::AbstractArray,
    dynamic_device::DynamicWrapper{DynG},
) where {DynG <: PSY.DynamicGenerator}
    global_index = get_global_index(dynamic_device)
    mass_matrix_tg_entries!(mass_matrix, PSY.get_prime_mover(dynamic_device), global_index)
    mass_matrix_pss_entries!(mass_matrix, PSY.get_pss(dynamic_device), global_index)
    mass_matrix_avr_entries!(mass_matrix, PSY.get_avr(dynamic_device), global_index)
    mass_matrix_machine_entries!(mass_matrix, PSY.get_machine(dynamic_device), global_index)
    mass_matrix_shaft_entries!(mass_matrix, PSY.get_shaft(dynamic_device), global_index)
    return
end

function device!(
    device_states::AbstractArray{T},
    output_ode::AbstractArray{T},
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    global_vars::AbstractArray{T},
    inner_vars::AbstractArray{T},
    dynamic_device::DynamicWrapper{DynG},
    t,
) where {DynG <: PSY.DynamicGenerator, T <: Real}
    if get_connection_status(dynamic_device) < 1.0
        output_ode .= zero(T)
        return
    end

    inner_vars[VR_gen_var] = voltage_r
    inner_vars[VI_gen_var] = voltage_i

    sys_ω = global_vars[GLOBAL_VAR_SYS_FREQ_INDEX]

    #Update Inner Vars
    _update_inner_vars!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)

    #Obtain ODEs and Mechanical Power for Turbine Governor
    mdl_tg_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)

    #Obtain ODEs for PSS
    mdl_pss_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)

    #Obtain ODEs for AVR
    mdl_avr_ode!(device_states, output_ode, inner_vars, dynamic_device)

    #Obtain ODEs for Machine
    mdl_machine_ode!(
        device_states,
        output_ode,
        inner_vars,
        current_r,
        current_i,
        dynamic_device,
    )

    #Obtain ODEs for PSY.Shaft
    mdl_shaft_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)
    return
end

function device!(
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    ::AbstractArray{T},
    ::AbstractArray{T},
    device::StaticWrapper{PSY.Source, U},
    t,
) where {T <: Real, U <: BusCategory}
    mdl_source!(voltage_r, voltage_i, current_r, current_i, device)
    return
end

function device!(
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    ::AbstractArray{T},
    ::AbstractArray{T},
    device::StaticWrapper{PSY.PowerLoad, U},
    t,
) where {T <: Real, U <: LoadCategory}
    # For now model all loads as constant impedance
    if get_connection_status(device) < 1.0
        current_r .= zero(T)
        current_i .= zero(T)
    else
        mdl_Zload!(voltage_r, voltage_i, current_r, current_i, device.device)
    end
    return
end

function device_mass_matrix_entries!(
    mass_matrix::AbstractArray{Float64},
    dynamic_device::DynamicWrapper{DynI},
) where {DynI <: PSY.DynamicInverter}
    global_index = get_global_index(dynamic_device)
    mass_matrix_DCside_entries!(
        mass_matrix,
        PSY.get_dc_source(dynamic_device),
        global_index,
    )
    mass_matrix_freq_estimator_entries!(
        mass_matrix,
        PSY.get_freq_estimator(dynamic_device),
        global_index,
    )
    mass_matrix_outer_entries!(
        mass_matrix,
        PSY.get_outer_control(dynamic_device),
        global_index,
    )
    mass_matrix_inner_entries!(
        mass_matrix,
        PSY.get_inner_control(dynamic_device),
        global_index,
    )
    mass_matrix_converter_entries!(
        mass_matrix,
        PSY.get_converter(dynamic_device),
        global_index,
    )
    mass_matrix_filter_entries!(mass_matrix, PSY.get_filter(dynamic_device), global_index)
    return
end

function device!(
    device_states::AbstractArray{T},
    output_ode::AbstractArray{T},
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    global_vars::AbstractArray{T},
    inner_vars::AbstractArray{T},
    dynamic_device::DynamicWrapper{DynI},
    t,
) where {DynI <: PSY.DynamicInverter, T <: Real}
    if get_connection_status(dynamic_device) < 1.0
        output_ode .= zero(T)
        return
    end

    #Obtain global vars
    sys_ω = global_vars[GLOBAL_VAR_SYS_FREQ_INDEX]

    #Update Voltage data
    inner_vars[Vr_inv_var] = voltage_r
    inner_vars[Vi_inv_var] = voltage_i

    #Update V_ref
    V_ref = get_V_ref(dynamic_device)
    inner_vars[V_oc_var] = V_ref

    #Update current inner_vars
    _update_inner_vars!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)

    #Obtain ODES for DC side
    mdl_DCside_ode!(device_states, output_ode, sys_ω, inner_vars, dynamic_device)

    #Obtain ODEs for PLL
    mdl_freq_estimator_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)

    #Obtain ODEs for OuterLoop
    mdl_outer_ode!(device_states, output_ode, inner_vars, sys_ω, dynamic_device)

    #Obtain inner controller ODEs and modulation commands
    mdl_inner_ode!(device_states, output_ode, inner_vars, dynamic_device)

    #Obtain converter relations
    mdl_converter_ode!(device_states, output_ode, inner_vars, dynamic_device)

    #Obtain ODEs for output filter
    mdl_filter_ode!(
        device_states,
        output_ode,
        current_r,
        current_i,
        inner_vars,
        sys_ω,
        dynamic_device,
    )

    return
end

function device_mass_matrix_entries!(
    mass_matrix::AbstractArray,
    dynamic_device::DynamicWrapper{PSY.PeriodicVariableSource},
)
    global_index = get_global_index(dynamic_device)
    mass_matrix_pvs_entries!(mass_matrix, dynamic_device, global_index)
    return
end

function mass_matrix_pvs_entries!(
    mass_matrix,
    pvs::DynamicWrapper{PSY.PeriodicVariableSource},
    global_index::ImmutableDict{Symbol, Int64},
)
    @debug "Using default mass matrix entries $pvs"
end

function device!(
    device_states::AbstractArray{T},
    output_ode::AbstractArray{T},
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    ::AbstractArray{T},
    ::AbstractArray{T},
    dynamic_device::DynamicWrapper{PSY.PeriodicVariableSource},
    t,
) where {T <: Real}
    ω_θ = PSY.get_internal_angle_frequencies(get_device(dynamic_device))
    ω_V = PSY.get_internal_angle_frequencies(get_device(dynamic_device))

    dV = 0
    for (ix, A) in
        enumerate(PSY.get_internal_voltage_coefficients(get_device(dynamic_device)))
        t <= 0 && continue
        dV += ω_V[ix] * (A[1] * cos(ω_V[ix] * t) - A[2] * sin(ω_V[ix] * t))
    end

    dθ = 0
    for (ix, A) in
        enumerate(PSY.get_internal_angle_coefficients(get_device(dynamic_device)))
        t <= 0 && continue
        dθ += ω_θ[ix] * (A[1] * cos(ω_θ[ix] * t) - A[2] * sin(ω_θ[ix] * t))
    end

    # Internal Voltage states
    V_R = device_states[1] * cos(device_states[2])
    V_I = device_states[1] * sin(device_states[2])
    output_ode[1] = dV
    output_ode[2] = dθ

    #update current
    R_th = PSY.get_R_th(get_device(dynamic_device))
    X_th = PSY.get_X_th(get_device(dynamic_device))
    Zmag = R_th^2 + X_th^2
    current_r[1] += R_th * (V_R - voltage_r[1]) / Zmag + X_th * (V_I - voltage_i[1]) / Zmag #in system pu flowing out
    current_i[1] += R_th * (V_I - voltage_i[1]) / Zmag - X_th * (V_R - voltage_r[1]) / Zmag #in system pu flowing out

    return
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, TG, P}},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #do nothing
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{PSY.DynamicGenerator{M, S, A, TG, P}},
) where {
    M <: Union{PSY.RoundRotorQuadratic, PSY.RoundRotorExponential},
    S <: PSY.Shaft,
    A <: PSY.AVR,
    TG <: PSY.TurbineGov,
    P <: PSY.PSS,
}
    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, M)

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ed_p = internal_states[2]
    ψ_kd = internal_states[3]
    ψ_kq = internal_states[4]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, M)
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #RI to dq transformation
    V_dq = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψq_pp = γ_q1 * ed_p + ψ_kq * (1 - γ_q1)
    ψd_pp = γ_d1 * eq_p + γ_d2 * (Xd_p - Xl) * ψ_kd
    ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)
    #Currents
    I_d =
        (1.0 / (R^2 + Xq_pp * Xd_pp)) *
        (-R * (V_dq[1] - ψq_pp) + Xq_pp * (-V_dq[2] + ψd_pp))
    Se = saturation_function(machine, ψ_pp)
    Xad_Ifd = eq_p + (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd + γ_d2 * eq_p) + Se * ψd_pp

    #Update Xad_Ifd
    inner_vars[Xad_Ifd_var] = Xad_Ifd
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.SalientPoleQuadratic, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, typeof(machine))

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ψ_kd = internal_states[2]
    ψq_pp = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, typeof(machine))
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #RI to dq transformation
    V_d, V_q = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd

    #Currents
    I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp) + Xd_pp * (ψd_pp - V_q))
    Se = saturation_function(machine, eq_p)
    Xad_Ifd =
        eq_p + Se * eq_p + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))

    #Update Xad_Ifd
    inner_vars[Xad_Ifd_var] = Xad_Ifd
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{
        PSY.DynamicGenerator{PSY.SalientPoleExponential, S, A, TG, P},
    },
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #Obtain indices for component w/r to device
    machine = PSY.get_machine(dynamic_device)
    local_ix = get_local_state_ix(dynamic_device, typeof(machine))

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    eq_p = internal_states[1]
    ψ_kd = internal_states[2]
    ψq_pp = internal_states[3]

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, typeof(machine))
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = inner_vars[VR_gen_var]
    V_tI = inner_vars[VI_gen_var]

    #Get parameters
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    #RI to dq transformation
    V_d, V_q = ri_dq(δ) * [V_tR; V_tI]

    #Additional Fluxes
    ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd
    ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)

    #Currents
    I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp) + Xq_pp * (-V_q + ψd_pp))
    Se = saturation_function(machine, ψ_pp)
    Xad_Ifd =
        eq_p + Se * ψd_pp + (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))

    #Update Xad_Ifd
    inner_vars[Xad_Ifd_var] = Xad_Ifd
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{PSY.DynamicInverter{C, O, IC, DC, P, F}},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}
    #do nothing
end

function _update_inner_vars!(
    device_states,
    output_ode,
    ω_sys,
    inner_vars,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{PSY.RenewableEnergyConverterTypeA, O, IC, DC, P, PSY.RLFilter},
    },
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}
    V_R = inner_vars[Vr_inv_var]
    V_I = inner_vars[Vi_inv_var]
    V_t = sqrt(V_R^2 + V_I^2)
    θ = atan(V_I / V_R)

    #Get Converter parameters
    converter = PSY.get_converter(dynamic_device)
    Brkpt = PSY.get_Brkpt(converter)
    Zerox = PSY.get_Zerox(converter)
    Lvpl1 = PSY.get_Lvpl1(converter)
    Vo_lim = PSY.get_Vo_lim(converter)
    Lv_pnt0, Lv_pnt1 = PSY.get_Lv_pnts(converter)
    K_hv = PSY.get_K_hv(converter)
    Lvpl_sw = PSY.get_Lvpl_sw(converter)
    R_source = PSY.get_R_source(converter)
    X_source = PSY.get_X_source(converter)

    #Define internal states for Converter
    converter_ix = get_local_state_ix(dynamic_device, PSY.RenewableEnergyConverterTypeA)
    converter_states = @view device_states[converter_ix]
    Ip = converter_states[1]
    Iq = converter_states[2]
    Vmeas = converter_states[3]

    #Saturate Ip if LVPL is active
    Ip_sat = Ip
    if Lvpl_sw == 1
        LVPL = get_LVPL_gain(Vmeas, Zerox, Brkpt, Lvpl1)
        Ip_sat = Ip <= LVPL ? Ip : LVPL
    end
    G_lv = get_LV_current_gain(V_t, Lv_pnt0, Lv_pnt1)
    Iq_extra = max(K_hv * (V_t - Vo_lim), 0.0)
    Id_cnv = G_lv * Ip_sat
    Iq_cnv = -Iq - Iq_extra
    #Reference Transformation
    Ir_cnv = Id_cnv * cos(θ) - Iq_cnv * sin(θ)
    Ii_cnv = Id_cnv * sin(θ) + Iq_cnv * cos(θ)

    #Obtain parameters
    filt = PSY.get_filter(dynamic_device)
    rf = PSY.get_rf(filt)
    lf = PSY.get_lf(filt)

    function V_cnv_calc(Ir_cnv, Ii_cnv, Vr_inv, Vi_inv)
        if lf != 0.0 || rf != 0.0
            Z_source_mag_sq = R_source^2 + X_source^2
            Zf_mag_sq = rf^2 + lf^2
            r_total_ratio = rf / Zf_mag_sq + R_source / Z_source_mag_sq
            l_total_ratio = lf / Zf_mag_sq + X_source / Z_source_mag_sq
            denom = 1.0 / (r_total_ratio^2 + l_total_ratio^2)
            rf_ratio = rf / Zf_mag_sq
            denom * (
                (Ir_cnv + Vi_inv * lf_ratio + Vr_inv * rf_ratio) * r_total_ratio -
                (Ii_cnv + Vi_inv * rf_ratio - Vr_inv * lf_ratio) * l_total_ratio
            )
            Vi_cnv =
                denom * (
                    (Ir_cnv + Vi_inv * lf_ratio + Vr_inv * rf_ratio) * l_total_ratio +
                    (Ii_cnv + Vi_inv * rf_ratio - Vr_inv * lf_ratio) * r_total_ratio
                )
        else
            Vr_cnv = Vr_inv
            Vi_cnv = Vi_inv
        end
        return Vr_cnv, Vi_cnv
    end

    Vr_cnv, Vi_cnv = V_cnv_calc(Ir_cnv, Ii_cnv, V_R, V_I)

    #Compute output currents
    if lf != 0.0 || rf != 0.0
        Vr_inv = V_R
        Vi_inv = V_I
        Zmag_squared = rf^2 + lf^2
        Ir_filt = (1.0 / Zmag_squared) * ((Vr_cnv - Vr_inv) * rf + (Vi_cnv - Vi_inv) * lf)
        Ii_filt = (1.0 / Zmag_squared) * ((Vi_cnv - Vi_inv) * rf - (Vr_cnv - Vr_inv) * lf)
    else
        Ir_filt = Ir_cnv
        Ii_filt = Ii_cnv
    end

    #Update inner_vars
    inner_vars[Ir_cnv_var] = Ir_cnv
    inner_vars[Ii_cnv_var] = Ii_cnv
    inner_vars[Vr_cnv_var] = Vr_cnv
    inner_vars[Vi_cnv_var] = Vi_cnv
    inner_vars[Ir_inv_var] = Ir_filt
    inner_vars[Ii_inv_var] = Ii_filt
end


text:
module Hierarchy

export Tree, parent_node

include("ancestors.jl")

# - A node is a dot-separated string representing hierarchy.
# - E.g. `""`, `a`, `a.b`, `a.bb.ccc`.
# - `a.b..c` & `a.b.c.` & `.a.b` are not valid.
# A parent node is the nearest ancestor node in `Tree`, fallback to `""` of not found.

struct Tree
    parent_node::Dict{String, String}
    # node => descendants
    # @assert has_node(t, descendant)
    placeholders::Dict{String, Vector{String}}

    Tree() = new(Dict{String, String}(), Dict{String, Vector{String}}())
end

has_node(t::Tree, node::String) = haskey(t.parent_node, node)
has_placeholder(t::Tree, node::String) = haskey(t.placeholders, node)
"""
For external usage only.
"""
parent_node(t::Tree, node::String) = t.parent_node[node]

"""
- Register `node` in `t`.
- Update parent nodes.
"""
function Base.push!(t::Tree, node::String)
    if has_node(t, node)
        return
    end

    if has_placeholder(t, node)
        # Fulfill a placeholder.
        descendants = pop!(t.placeholders, node)
        fix_descendants!(t, node, descendants)
    end

    parent_node = prepare_parent!(t, node)
    t.parent_node[node] = parent_node
end

"""
- `@assert has_node(t, node)`.
- Return the parent node of `node`.
- Register `node` as descendant of all its ancestor placeholders up to its parent node,
- since the further ancestor placeholders do not affect the parent node of `node` if they are fulfilled.
"""
function prepare_parent!(t::Tree, node::String)::String
    # Default.
    parent_node = ""

    for ancestor in Ancestors(node)
        if has_node(t, ancestor)
            # Hierarchy of `ancestor` is complete.
            # Found `parent_node`.
            parent_node = ancestor
            return parent_node
        end

        if has_placeholder(t, ancestor)
            # Register as a descendant.
            push!(t.placeholders[ancestor], node)
        else
            # Create a `Placeholder`.
            t.placeholders[ancestor] = String[node]
        end
    end

    parent_node
end

"""
- `@assert has_placeholder(t, placeholder)`.
- Replace the parent nodes of the descendants with `node` if `is_ancestor_or_self(parent_node, node)`
"""
function fix_descendants!(t::Tree, node::String, descendants::Vector{String})
    for descendant in descendants
        parent_node = t.parent_node[descendant]
        # Here we don't need strict ancestor check as it is slower.
        if is_ancestor_or_self(parent_node, node)
            t.parent_node[descendant] = node
        end
    end
end

end

max_stars_repo_name: colinfang/MiniLogging

max_issues_repo_name: colinfang/MiniLogging

max_forks_repo_name: colinfang/MiniLogging


text:
export LabeledTensor

# abstractions
abstract type AbstractTensor{T, N} end

"""
    LabeledTensor{T,N,AT<:AbstractArray{T,N}, LT, MT} <: AbstractTensor{T, N}
    LabeledTensor(array::AbstractArray, labels::AbstractVector[, meta])

Tensor with labeled legs.
When multiplying two tensors, legs with same labels will be treated as inner degree of freedom and get contracted.

The optional argument `meta` is the meta information, default is `nothing`. One can use it to store additional information like plotting layout.
"""
struct LabeledTensor{T,N,AT<:AbstractArray{T,N}, LT, MT} <: AbstractTensor{T, N}
    array::AT
    labels::Vector{LT}
    meta::MT
end
Base.ndims(t::LabeledTensor) = ndims(t.array)
LinearAlgebra.norm(t::LabeledTensor, p::Real=2) = norm(t.array, p)

function LabeledTensor(tensor::AbstractArray, labels::AbstractVector)
    @assert ndims(tensor) == length(labels) "dimension of tensor $(ndims(tensor)) != number of labels $(length(labels))"
    LabeledTensor(tensor, labels, nothing)
end

function Base.:(*)(A::LabeledTensor, B::LabeledTensor)
    labels_AB = setdiff(A.labels, B.labels) ∪ setdiff(B.labels, A.labels)
    LabeledTensor(tensorcontract(A.labels, A.array, B.labels, B.array, labels_AB), labels_AB, merge_meta(A.meta, B.meta))
end

function Base.isapprox(a::LabeledTensor, b::LabeledTensor; kwargs...)
    isapprox(a.array, b.array; kwargs...) && a.labels == b.labels
end
Base.size(t::LabeledTensor) = Base.size(t.array)
Base.copy(t::LabeledTensor) = LabeledTensor(copy(t.array), t.labels)
Base.similar(::Type{<:LabeledTensor}, arr::AbstractArray, labels::AbstractVector, meta=nothing) = LabeledTensor(arr, labels, meta)
Base.similar(::LabeledTensor, arr::AbstractArray, labels::AbstractVector, meta=nothing) = LabeledTensor(arr, labels, meta)
LinearAlgebra.rmul!(t::LabeledTensor, factor) = (rmul!(t.array, factor); t)

function Base.show(io::IO, lt::LabeledTensor)
    print(io, "$(typeof(lt).name.name){$(eltype(lt.array))}($(join(lt.labels, ", ")))")
end

function Base.show(io::IO, ::MIME"plain/text", lt::LabeledTensor)
    Base.show(io, lt)
end

function mul_dim(t::LabeledTensor, m::AbstractMatrix; dim::Int)
    data = t.array
    iA = ntuple(i->i, ndims(data))
    iB = (dim, -dim)
    iC = ntuple(i->i==dim ? -dim : i, ndims(data))
    LabeledTensor(tensorcontract(iA, data, iB, m, iC), t.labels)
end

struct PlotMeta
    loc::Tuple{Float64, Float64}
    name::String
end
merge_meta(m1::PlotMeta, m2::PlotMeta) = PlotMeta((m1.loc .+ m2.loc) ./ 2, m1.name*m2.name)
merge_meta(m1::Nothing, m2::Nothing) = nothing
dispmeta(m::PlotMeta) = m.name


text:
# interface to C++ library "figtree" from
#   Vlad I. Morariu, Balaji Vasan Srinivasan, Vikas C. Raykar,
#   Ramani Duraiswami, and Larry S. Davis. Automatic online tuning for
#   fast Gaussian summation. Advances in Neural Information Processing
#   Systems (NIPS), 2008.
# and website that hosts code is
#   http://www.umiacs.umd.edu/~morariu/figtree/

"""
    fgt!(v,d,M,N,h,ε,x,y,q,W)

Calls Fast Gauss Transform code compiled into shared library. Automatically chooses
   'speedup' parameters

Notes:

INPUTS:
v: [Vector{Float64}, M] interpolated values
d: [Int] dimension of data
M: [Int] number of targets
N: [Int] number of sources
h: [Float64] bandwidth, h=sqrt(2)*stddev if standard deviation is known
ε: [Float64] maximum absolute error
x: [Vector{Float64}, N×d] list of sources
y: [Vector{Float64}, M×d] list of targets
q: [Vector{Float64}, N] weights on sources
W: [Int] number of weight functions (e.g. if W=2, then length(q) = 2N)
"""
function fgt!(v::Array{Float64},d::Int,M::Int,N::Int,h::Float64,ε::Float64,
              x::Array{Float64},y::Array{Float64},q::Array{Float64},W::Int)
    # zero-out v output array
    fill!(v,zero(Float64))

    # call c function
    ccall((:figtree,libfigtree),
            Int,
            (Int, Int, Int, Int, Ptr{Float64}, Float64, Ptr{Float64}, Ptr{Float64}, Float64, Ptr{Float64}, Int, Int, Int, Int),
            d, N, M, W, x, h, q, y, ε, v, 4, 1, 0, 0)       # 0 indicates direct evaluation, 4 is auto
                                                   # see figtree.h to learn last 4 integer codes -- DO NOT TOUCH
                                                   # -1 just used to avoid selecting any value

    # v is modified by reference
    nothing
end
# same as fgt but uses slow algorithm (classic matrix-vector multiply)
function dgt!(v::Array{Float64},d::Int,M::Int,N::Int,h::Float64,ε::Float64,
              x::Array{Float64},y::Array{Float64},q::Array{Float64},W::Int)
    # zero-out v output array
    fill!(v,zero(Float64))

    # call c function
    ccall((:figtree,libfigtree),
            Int,
            (Int, Int, Int, Int, Ptr{Float64}, Float64, Ptr{Float64}, Ptr{Float64}, Float64, Ptr{Float64},Int, Int, Int, Int),
            d, N, M, W, x, h, q, y, ε, v, 0, 0, 0, 0)    # 0 indicates direct evaluation, 4 is auto
                                                   # see figtree.h to learn last 4 integer codes -- DO NOT TOUCH
                                                   # -1 just used to avoid selecting any value

    # v is modified by reference
    nothing
end

function mydgt!(v::Array{Float64},d::Int,M::Int,N::Int,h::Float64,ε::Float64,
                x::Array{Float64},y::Array{Float64},q::Array{Float64},W::Int)
    # zero-out v output array
    fill!(v,zero(Float64))

    # generate nodes
    for ti=1:M              # targets
        for tj=1:N          # sources
            dist = sqrt((x[tj] - y[ti])^2)
            v[ti] += q[tj]*GaussKernel(dist,h)
        end
    end

    nothing
end

function GaussKernel(x,h)
    return exp(-x^2/h^2)
end


text:
"""
    struct ConstantImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}

Impurity density model that assumes a constant impurity density everywhere.

## Fields
* `ρ::T`: the constant value of the impurity density.

## Definition in Configuration File

A `ConstantImpurityDensity` is defined in the configuration file through the field `impurity_density`
(of a `semiconductor`) with `name: constant` and a `value` for `ρ`.

An example definition of a constant impurity density looks like this:
```yaml
impurity_density:
  name: constant
  value: 1.0e10 # 1/m³
```
"""
struct ConstantImpurityDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
    ρ::T
end

function get_impurity_density(cdm::ConstantImpurityDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    return cdm.ρ
end

function ImpurityDensity(T::DataType, t::Val{:constant}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, input_units::NamedTuple)
    ρ::T = haskey(dict, "value") ? _parse_value(T, dict["value"], input_units.length^(-3)) : T(0)
    ConstantImpurityDensity{T}( ρ )
end


text:
using FreeType
using Base.Test

const ok = 0

library = Array(FT_Library, 1)
@test ok == FT_Init_FreeType(library)

aface = Array(FT_Face, 1)
@unix_only begin
  @test ok == FT_New_Face(library[1], "/Library/Fonts/Verdana.ttf", 0, aface)
end

pface = aface[1]
@test ok == FT_Select_Charmap(pface, FreeType.FT_ENCODING_UNICODE)

units_per_EM = 0x0800
@test ok == FT_Set_Char_Size(pface, 0, units_per_EM, 72, 72)

glyph_index = FT_Get_Char_Index(pface, 'A')
@test ok == FT_Load_Glyph(pface, glyph_index, FT_LOAD_DEFAULT)

face = unsafe_load(pface)

@test 0x0800 == face.units_per_EM
@test 2059 == face.ascender
@test -430 == face.descender
@test 2489 == face.height
@test FreeType.FT_BBox(-1013,-621,2963,2049) == face.bbox
@test 913 == face.num_glyphs
@test "Verdana" == bytestring(face.family_name)
@test "Regular" == bytestring(face.style_name)
@test 2 == face.num_charmaps

glyph = unsafe_load(face.glyph)
@test FreeType.FT_GLYPH_FORMAT_OUTLINE == glyph.format
@test FreeType.FT_Glyph_Metrics(1408,1536,0,1536,1408,-704,192,1984) == glyph.metrics

outline = glyph.outline
@test 2 == outline.n_contours
@test 11 == outline.n_points

@test [
 FreeType.FT_Vector(1408,0)
 FreeType.FT_Vector(1203,0)
 FreeType.FT_Vector(1067,384)
 FreeType.FT_Vector(341,384)
 FreeType.FT_Vector(205,0)
 FreeType.FT_Vector(0,0)
 FreeType.FT_Vector(566,1536)
 FreeType.FT_Vector(842,1536)
 FreeType.FT_Vector(997,576)
 FreeType.FT_Vector(704,1383)
 FreeType.FT_Vector(409,576)
] == map(n->unsafe_load(outline.points, n), 1:outline.n_points)


text:
using Plots
using Reproduce
using JLD2
using FileIO
using Statistics

include("plot_utils.jl")

function get_lines(data)
    means = []
    step = 10
    for start in 1:step:size(data)[2]
        chunk = data[:, start:min(start + step - 1, length(data)),: ]
        # println(size(chunk))
        chunk_mean = mean(chunk, dims=2)
        push!(means, chunk_mean)
    end

    smoothed = cat(means..., dims=2)[:,:, :]
    sum_lines = sum(smoothed, dims=1)
    num_runs = size(sum_lines)[3]
    mean_line = dropdims(mean(sum_lines, dims=3), dims=(1, 3))
    std_err_line = dropdims(std(sum_lines, dims=3), dims=(1, 3)) / sqrt(num_runs)

    xs = 1:100*step:30000
    return xs, mean_line, std_err_line
end

num_steps = 30000
log_interval = 100

save_file = "plotting/chunlok/generated_plots/reproduce_plot_visitation.svg"
p = plot()

ic = ItemCollection("./Test_Large_Distractor_GPI/")

println(diff(ic))

data = load_results(ic, :visit_order, return_type="array")

min_length = min([length(x) for x in data]...)

println(min_length)
# Trimming extra episodes for plotting sake
data = [x[1:min_length] for x in data]

data = cat(data..., dims = 2)
runs = size(data)[2]

percentages = []
for r in 1:runs
    percentage = get_single_goal_percentage(data[:,r]; step_size=2)
    push!(percentages, percentage)
end

percentages = cat(percentages..., dims = 3)
percentages = dropdims(mean(percentages, dims=3), dims=3)
println(size(percentages))

# sdfsdfsdf
for i in 1:4
    println(size(percentages[i, :]))
    plot!(p, 1:size(percentages)[2], percentages[i, :])
end

# print(error[3,10])
savefig(save_file)

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function LoessInterpolator()
    return LoessInterpolator(())
end

function LoessInterpolator(arg0::jdouble, arg1::jint)
    return LoessInterpolator((jdouble, jint), arg0, arg1)
end

function LoessInterpolator(arg0::jdouble, arg1::jint, arg2::jdouble)
    return LoessInterpolator((jdouble, jint, jdouble), arg0, arg1, arg2)
end

function interpolate(obj::LoessInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", PolynomialSplineFunction, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function smooth(obj::LoessInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "smooth", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function smooth(obj::LoessInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return jcall(obj, "smooth", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end


text:
using Logging

@testset "mr_load integration test" begin
    # TODO: This is far too trivial! should get / generate some better test data.
    cosy = @test_logs (:warn,r"Different values found in MeasYaps") #=
        =# (:warn,r"Sequence.*is not one of the versions I recognize"s) #=
        =# mr_load("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")

    @test size(cosy.lcosy_scans) == (1,1)
    @test size(cosy.ref_scans) == (0,)
    @test step(cosy.t1) == 0.8u"ms"
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
numbers = [convert(Char, string(i)[1]) for i in 0:9]
push!(numbers, '-', 'e', '.')

function lex_string!(lines::Vector{Char})

    if lines[1] == JSON_QUOTE
        if length(lines) > 1
            popnfirst!(lines, 1)
        end
    else
        return nothing
    end # if

    jsonString = ""
    for c in lines
        if c == JSON_QUOTE
            if length(lines) > 1
                popnfirst!(lines, length(jsonString)+1)
            end
            return jsonString
        else
            jsonString *= c
        end # if
    end # for

    error("No end of string found")

end # lex_string

function lex_number!(lines::Vector{Char})
    jsonNum = ""

    for c in lines
        if c in numbers
            jsonNum *= c
        else
            break
        end # if
    end # for

    popnfirst!(lines, length(jsonNum))

    if length(jsonNum) == 0
        return nothing
    end # if

    if '.' in jsonNum
        return parse(Float64, jsonNum)
    end # if

    return parse(Int64, jsonNum)

end # lex_number

function lex_bool!(lines::Vector{Char})

    if length(lines) >= 4
        if lines[:4] == "true"
            popnfirst!(lines,4)
            return true
        elseif lines[:4] == "null"
            popnfirst!(lines,4)
            return nothing
        end
    elseif length(lines) >= 5 && lines[1:5] == "false"
        popnfirst!(lines,5)
        return false
    end # if

    return nothing
end # lex_bool

function lexer(filename::String)
    lines = Vector{Char}(Mmap.mmap(filename, Vector{UInt8}, filesize(filename)))

    tokens = Vector{Any}()

    while length(lines) > 0
        jsonString = lex_string!(lines)
        if jsonString !== nothing
            push!(tokens, jsonString)
            continue
        end # if

        jsonNum = lex_number!(lines)
        if jsonNum !== nothing
            push!(tokens, jsonNum)
            continue
        end # if

        jsonBool = lex_bool!(lines)
        if jsonBool !== nothing
            push!(tokens, jsonBool)
            continue
        end # if

        if lines[1] in JSON_WHITESPACE
            popfirst!(lines)
        elseif lines[1] in JSON_SYNTAX
            push!(tokens, popfirst!(lines))
        end # if
    end # while

    return tokens
end # lexer

max_stars_repo_name: Iaggelis/jjson

max_issues_repo_name: Iaggelis/jjson

max_forks_repo_name: Iaggelis/jjson


text:
using Test, SDCO

@testset "SDCO instance, feasibility, objectives" begin

    @testset "testcase1a sparse - symmetric = $issym" for issym in [true, false]
        pb, z = testcase1a(symmetric = issym);

        @test evaluate(pb.A, z.x) == [1., 1., 1.]

        @test get_primobj(pb, z) == 3.
        @test get_dualobj(pb, z) == 0.

        @test get_primslacks(pb, z) == [0., 0., 0.]
        @test norm(get_dualslacks(pb, z.y, z.s)) == 0

        @test get_primfeaserr(pb, z) == 0.
        @test get_dualfeaserr(pb, z) == 0.
    end
end

@testset "Nesterov Todd step" begin

    @testset "testcase1a sparse - symmetric = $issym" for issym in [true, false]
        pb, z = testcase1a(symmetric = issym)

        @testset "mu = $mu" for mu in [0.2 0.4 0.6]
            dz = NesterovToddstep(pb, z, mu)

            @test norm(dz.x - PointE([ (mu-1)*[0. -1/3 -1/3 ; -1/3 2/3 0.; -1/3 0. 2/3] ], Float64[])) < 1e-15
            @test norm(dz.y - (1-mu) .* [1., 1/3, 1/3]) < 1e-15
            @test norm(dz.s - PointE([ (mu-1)*[1. 1/3 1/3 ; 1/3 1/3 0.; 1/3 0. 1/3] ], Float64[])) < 1e-15
        end
    end

    @testset "testcase1b" begin
        pb, z = testcase1b()

        @testset "mu = $mu" for mu in [1. 0.2 0.4 0.6]
            dz = NesterovToddstep(pb, z, mu)

            @test norm(dz.y - [0.6 * (1-3*mu)]) < 1e-15

            @test norm(dz.s - PointE(SDCO.Dense{Float64}[], Float64[ 0.6*(3*mu-1), 1.2*(3*mu-1) ])) < 1e-15

            @test norm(dz.x - PointE(SDCO.Dense{Float64}[], Float64[ 0.4*mu - 2/15, -0.2*mu + 1/15 ])) < 1e-15
        end
    end

end

@testset "Central path following" begin

    @testset "testcase1a sparse - symmetric = $issym" for issym in [true, false]
        pb, z = testcase1a(symmetric = issym)

        @assert delta(pb, z) < 1/3

        zfinal = NTpredcorr(pb, z; outlev=0)

        epsilon = 1e-10
        @test mu(zfinal) < 1e-10

        @test get_primfeaserr(pb, zfinal) < epsilon
        @test get_dualfeaserr(pb, zfinal) < epsilon
        @test dot(zfinal.x, zfinal.s) < epsilon
    end

    @testset "testcase1b" begin
        pb, z = testcase1b()

        @assert delta(pb, z) < 1/3

        zfinal = NTpredcorr(pb, z; outlev=0)

        epsilon = 1e-10
        @test mu(zfinal) < 1e-10

        @test get_primfeaserr(pb, zfinal) < epsilon
        @test get_dualfeaserr(pb, zfinal) < epsilon
        @test dot(zfinal.x, zfinal.s) < epsilon
    end
end

@testset "Getting central point" begin
    @testset "testcase1c sparse - symmetric = $issym" for issym in [true, false]
        pb, z = testcase1c(symmetric = issym)

        @assert delta(pb, z) > 1/3

        zcentralpath = NTgetcentralpathpoint(pb, z; outlev=0)

        @test delta(pb, zcentralpath) < 1/3
    end
end

@testset "Minimum Matrix Norm, feasible starting point" begin
    @testset "mmn-$p-$q-$r" for (p, q, r) in [ (2, 2, 2), (4, 2, 2) ]
        ε = 1e-10

        pb, z = testcase2(p, q, r)
        pb.options[:opt_ε] = ε
        pb.options[:opt_outlev] = 0

        zfinal = solve(pb, z)

        @test mu(zfinal) < ε

        # @test get_primfeaserr(pb, zfinal) < epsilon
        # @test get_dualfeaserr(pb, zfinal) < epsilon
        # @test dot(zfinal.x, zfinal.s) < epsilon
    end
end


text:
# Packages required for testing
using Test
using Random
using Logging
using LinearAlgebra
using Distributions
using NBInclude

# Package Under Test
using GaussianFilters

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Fix randomness during tests
Random.seed!(0)

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test isapprox(a,b, rtol=rtol, atol=atol)
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        @test isapprox(a, y, rtol=rtol, atol=atol)
    end
end

# Define package tests

@time @testset "GaussianFilter Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")

    @time @testset "Models Tests" begin
        include(joinpath(testdir, "test_models.jl"))
    end

    @time @testset "GaussianFilter GM-PHD Testing" begin
        include(joinpath(testdir, "test_gmphd.jl"))
    end

    @time @testset "GaussianFilter KF Testing" begin
        include(joinpath(testdir, "test_kf.jl"))
    end

    @time @testset "GaussianFilter EKF Testing" begin
        include(joinpath(testdir, "test_ekf.jl"))
    end

    @time @testset "GaussianFilter UKF Testing" begin
        include(joinpath(testdir, "test_ukf.jl"))
    end

    @testset "Notebooks testing" begin
        nbdir = joinpath(dirname(@__DIR__), "notebooks")
        for d in readdir(nbdir)
            if endswith(d, ".ipynb")
                path = joinpath(nbdir, d)
                @testset "$d" begin
                    @nbinclude path
                end
            end
        end
    end

end


text:
"""
    setup(dir)

Initialise a `Publish` project in the given directory `dir`, which is created
if it does not already exist.

When `dir` has a Julia project structure with a `Project.toml` file then
`Publish` is added to it's dependencies list, otherwise, if it is not a Julia
package then a `Project.toml` file is created as well as a `README.md` file.

## Keywords

  - `pkg::Bool=true` can be set to `false` to disable all `Pkg` calls during
    use of `setup`.
"""
function setup(dir; pkg=true)
    isdir(dir) || mkpath(dir)
    cd(dir) do
        _, name = splitdir(pwd())
        ## Project setup.
        proj = "Project.toml"
        if isfile(proj)
            @info "found '$proj'. Using it."
        else
            @info "no '$proj' found. Creating one."
            uuid = string(Base.identify_package("Publish").uuid)
            toml = """
            name = "$name"

            [deps]
            Publish = "$uuid"
            """
            write(proj, toml)
        end
        if pkg
            @info "Activating and instantiating '$proj'."
            Pkg.activate(proj)
            Pkg.instantiate()
            @info "Adding 'Publish' to project dependencies."
            Pkg.add("Publish")
        end
        ## Content setup.
        readme = "README.md"
        if isfile(readme)
            @info "found '$readme'."
        else
            @info "no '$readme' found, creating one."
            write(readme, "# $name\n")
        end
        ## Create github actions workflow.
        action = ".github/workflows/Publish.yml"
        if isfile(action)
            @warn "This project already has a 'Publish.yml' workflow."
        else
            @info "Creating a GitHub Actions workflow for deployment."
            mkpath(dirname(action))
            contents = read(joinpath(@__DIR__, "templates", "Publish.yml"), String)
            open(action, "w") do handle
                write(handle, contents)
            end
        end
        ## Display a message for users to help them start a server.
        path = joinpath(pwd(), proj)
        @info "run `serve($(repr(path)))` to start a server."
    end
end


text:
# imports
using FMI
using FMIFlux
using Flux
using DifferentialEquations: Tsit5
import Plots

simpleFMUPath = joinpath(dirname(@__FILE__), "../model/SpringPendulum1D.fmu")
realFMUPath = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")
println("SimpleFMU path: ", simpleFMUPath)
println("RealFMU path: ", realFMUPath)

tStart = 0.0
tStep = 0.01
tStop = 5.0
tSave = collect(tStart:tStep:tStop)

realFMU = fmiLoad(realFMUPath)
fmiInstantiate!(realFMU; loggingOn=false)
fmiInfo(realFMU)

fmiSetupExperiment(realFMU, tStart, tStop)

fmiEnterInitializationMode(realFMU)
fmiExitInitializationMode(realFMU)

x₀ = fmiGetContinuousStates(realFMU)

vrs = ["mass.s", "mass.v", "mass.a", "mass.f"]
_, realSimData = fmiSimulate(realFMU, tStart, tStop; recordValues=vrs, saveat=tSave, setup=false, reset=false)
fmiPlot(realFMU, vrs, realSimData)

fmiUnload(realFMU)

velReal = collect(data[2] for data in realSimData.saveval)
posReal = collect(data[1] for data in realSimData.saveval)

simpleFMU = fmiLoad(simpleFMUPath)

fmiInstantiate!(simpleFMU; loggingOn=false)
fmiInfo(simpleFMU)

vrs = ["mass.s", "mass.v", "mass.a"]
_, simpleSimData = fmiSimulate(simpleFMU, tStart, tStop; recordValues=vrs, saveat=tSave, reset=false)
fmiPlot(simpleFMU, vrs, simpleSimData)

velSimple = collect(data[2] for data in simpleSimData.saveval)
posSimple = collect(data[1] for data in simpleSimData.saveval)

# loss function for training
function lossSum()
    solution = neuralFMU(x₀, tStart)

    posNet = collect(data[1] for data in solution.u)
    #velNet = collect(data[2] for data in solution.u)

    Flux.Losses.mse(posReal, posNet) #+ Flux.Losses.mse(velReal, velNet)
end

# callback function for training
global counter = 0
function callb()
    global counter += 1

    if counter % 20 == 1
        avgLoss = lossSum()
        @info "Loss [$counter]: $(round(avgLoss, digits=5))   Avg displacement in data: $(round(sqrt(avgLoss), digits=5))"
    end
end

# NeuralFMU setup
numStates = fmiGetNumberOfStates(simpleFMU)

net = Chain(inputs -> fmiEvaluateME(simpleFMU, inputs),
            Dense(numStates, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))

neuralFMU = ME_NeuralFMU(simpleFMU, net, (tStart, tStop), Tsit5(); saveat=tSave);

solutionBefore = neuralFMU(x₀, tStart)
fmiPlot(simpleFMU, solutionBefore)

# train
paramsNet = Flux.params(neuralFMU)

optim = ADAM()
Flux.train!(lossSum, paramsNet, Iterators.repeated((), 300), optim; cb=callb)

# plot results mass.s
solutionAfter = neuralFMU(x₀, tStart)

fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
                 xtickfontsize=12, ytickfontsize=12,
                 xguidefontsize=12, yguidefontsize=12,
                 legendfontsize=8, legend=:topright)

posNeuralFMU = collect(data[1] for data in solutionAfter.u)

Plots.plot!(fig, tSave, posSimple, label="SimpleFMU", linewidth=2)
Plots.plot!(fig, tSave, posReal, label="RealFMU", linewidth=2)
Plots.plot!(fig, tSave, posNeuralFMU, label="NeuralFMU (300 epochs)", linewidth=2)
fig

Flux.train!(lossSum, paramsNet, Iterators.repeated((), 700), optim; cb=callb)
# plot results mass.s
solutionAfter = neuralFMU(x₀, tStart)
posNeuralFMU = collect(data[1] for data in solutionAfter.u)
Plots.plot!(fig, tSave, posNeuralFMU, label="NeuralFMU (1000 epochs)", linewidth=2)
fig

fmiUnload(simpleFMU)