text:
using Test, Onda, Dates

@testset "basic TimeSpan code paths" begin
    t = TimeSpan(Nanosecond(rand(UInt32)))
    @test t == TimeSpan(t)
    @test t == TimeSpan(first(t), last(t))
    @test contains(t, t)
    @test overlaps(t, t)
    @test shortest_timespan_containing([t]) == t
    @test shortest_timespan_containing((t,t,t)) == t
    @test duration(t) == Nanosecond(0)
    @test duration(first(t)) == first(t)
    @test_throws ArgumentError TimeSpan(4, 2)
end

@testset "contains(::TimeSpan...)" begin
    @test contains(TimeSpan(10, 20), TimeSpan(10, 20))
    @test contains(TimeSpan(10, 20), TimeSpan(11, 19))
    @test contains(TimeSpan(11, 20), TimeSpan(11, 19))
    @test contains(TimeSpan(10, 19), TimeSpan(11, 19))
    @test !contains(TimeSpan(10, 20), TimeSpan(11, 21))
    @test !contains(TimeSpan(11, 20), TimeSpan(10, 19))
    @test !contains(TimeSpan(10, 19), TimeSpan(10, 21))
    @test !contains(TimeSpan(11, 19), TimeSpan(10, 20))
end

@testset "overlaps(::TimeSpan...)" begin
    @test overlaps(TimeSpan(10, 20), TimeSpan(10, 20))
    @test overlaps(TimeSpan(10, 20), TimeSpan(11, 19))
    @test overlaps(TimeSpan(11, 20), TimeSpan(11, 19))
    @test overlaps(TimeSpan(10, 19), TimeSpan(11, 19))
    @test overlaps(TimeSpan(10, 20), TimeSpan(11, 21))
    @test overlaps(TimeSpan(11, 20), TimeSpan(10, 19))
    @test overlaps(TimeSpan(10, 19), TimeSpan(10, 21))
    @test overlaps(TimeSpan(11, 19), TimeSpan(10, 20))
    @test overlaps(TimeSpan(10, 20), TimeSpan(20, 30))
    @test overlaps(TimeSpan(20, 30), TimeSpan(10, 20))
    @test !overlaps(TimeSpan(10, 20), TimeSpan(21, 30))
    @test !overlaps(TimeSpan(21, 30), TimeSpan(10, 20))
end

@testset "shortest_timespan_containing(spans)" begin
    @test shortest_timespan_containing([TimeSpan(1, 2),
                                        TimeSpan(5, 10),
                                        TimeSpan(2, 3)]) == TimeSpan(1, 10)
    @test shortest_timespan_containing([TimeSpan(3, 7),
                                        TimeSpan(1, 10),
                                        TimeSpan(2, 5)]) == TimeSpan(1, 10)
end

@testset "time <--> index conversion" begin
    @test_throws ArgumentError time_from_index(200, 0)
    @test time_from_index(100, 1) == Nanosecond(0)
    @test time_from_index(100, 301:600) == TimeSpan(Second(3), Second(6))
    @test time_from_index(100, 101:101) == TimeSpan(Second(1), Second(1))
    @test_throws ArgumentError index_from_time(200, Nanosecond(-1))
    @test index_from_time(100, Nanosecond(0)) == 1
    @test index_from_time(100, TimeSpan(Second(3), Second(6))) == 301:600
    @test index_from_time(100, TimeSpan(Second(1), Second(1))) == 101:101
    # test non-integer sample rates
    rate = 100.66
    ns_per_sample = Onda.nanoseconds_per_sample(rate)
    for i in 1:1000
        t = Nanosecond(ceil(Int, (i - 1) * ns_per_sample))
        @test index_from_time(rate, t) == i
        @test time_from_index(rate, i) == t
    end
end


text:
#=
P23 (**) Extract a given number of randomly selected elements from a list.
The selected items shall be returned in a list.
Example:
* (rnd-select '(a b c d e f g h) 3)
(E D A)

Hint: Use the built-in random number generator and the result of problem P20.
=#

function rnd_select(coll, n)
    [coll[rand(1:length(coll))] for _ in 1:n]
end

coll = [:a, :b, :c, :d, :e, :f, :g, :h]
rnd_selected = rnd_select(coll, 3)
@assert length(rnd_selected) == 3
@assert all(x -> x in coll, rnd_selected)

println("Tests passed: JL-23.jl")

ext: jl


text:
using PackageCompiler

pkgs = ["NeuralPDE", "Flux", "ModelingToolkit", "GalacticOptim", "Optim", "DiffEqFlux",
        "Quadrature", "Cuba", "CUDA", "QuasiMonteCarlo"]
create_sysimage(pkgs; sysimage_path="PINNSysimage.so",
                      precompile_execution_file="run-slim.jl")

max_stars_repo_name: jpsamaroo/DaggerFluxPINNs

max_issues_repo_name: jpsamaroo/DaggerFluxPINNs

max_forks_repo_name: jpsamaroo/DaggerFluxPINNs


text:
Threads.nthreads()

using LinearAlgebra, Random, Distributions
using MixedModelsBLB, BenchmarkTools, Profile
using KNITRO

N = 10^2
reps = 20000
p = 20 # number of fixed effects
q = 2 # number of random effects
βtrue  = ones(p)
σ²true = 1.5
σtrue  = sqrt(σ²true)
Σtrue  = [1.0 0; 0 3.0] # Matrix(Diagonal([2.0; 1.2; 1.0]))
Ltrue  = Matrix(cholesky(Symmetric(Σtrue)).L)
fenames = vcat("Intercept", "x" .* string.([1:1:(p-1);]))
renames = ["Intercept", "z1"]

# generate data
obsvec = Vector{blblmmObs{Float64}}(undef, N)
@inbounds for i in 1:N
    # initialize arrays
    X = Matrix{Float64}(undef, reps, p)
    X[:, 1] = ones(reps)
    Z = Matrix{Float64}(undef, reps, q)
    Z[:, 1] = ones(reps)
    storage_q = Vector{Float64}(undef, q)
    re_storage = Vector{Float64}(undef, q)
    y = Vector{Float64}(undef, reps)
    randn!(y) # y = standard normal error
    @views randn!(X[:, 2:p])
    BLAS.gemv!('N', 1., X, βtrue, σtrue, y) # y = Xβ + σtrue * standard normal error
    randn!(storage_q)
    BLAS.gemv!('N', 1., Ltrue, storage_q, 0., re_storage)
    @views randn!(Z[:, 2:q]) #Distributions.rand!(Normal(), Z[:, 2:q])
    BLAS.gemv!('N', 1., Z, re_storage, 1., y) # y = Xβ + Zα + error
    # y = X * βtrue .+ Z * (Ltrue * randn(q)) .+ σtrue * randn(ns[i])
    obsvec[i] = blblmmObs(y, X, Z)
end

# construct two models
m1 = blblmmModel(obsvec, fenames, renames, N, true);
m2 = blblmmModel(obsvec, fenames, renames, N, false);

# solver = KNITRO.KnitroSolver(outlev=0)
solver = Ipopt.IpoptSolver(print_level=0, max_iter=100, mehrotra_algorithm = "yes", warm_start_init_point = "yes", warm_start_bound_push = 1e-9)

@benchmark MixedModelsBLB.fit!($m1, solver) setup = init_ls!(m1)
# median time:      6.535 s

@benchmark MixedModelsBLB.fit!($m2, solver) setup = init_ls!($m2)
# median time:      685.603 ms

function update_logl_multithreaded!(
    m::blblmmModel{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: BlasReal
    Threads.@threads for obs in m.data
        loglikelihood!(obs, m.β, m.σ², m.ΣL, needgrad, needhess)
    end
end
function update_logl!(
    m::blblmmModel{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: LinearAlgebra.BlasReal
    for obs in m.data
        loglikelihood!(obs, m.β, m.σ², m.ΣL, needgrad, needhess)
    end
end

BLAS.set_num_threads(2)
@benchmark MixedModelsBLB.update_logl_multithreaded!($m1, false, false) setup = init_ls!($m1)
# median time:      1.271 ms
BLAS.set_num_threads(1)
@benchmark MixedModelsBLB.update_logl_multithreaded!($m1, false, false) setup = init_ls!($m1)
# median time:      1.913 ms

Profile.clear()
@profile MixedModelsBLB.update_logl_multithreaded!(m1, false, false)
Profile.print(format = :flat)

BLAS.set_num_threads(2)
@benchmark update_logl!($m1, false, false) setup = init_ls!($m1)
# median time:      99.543 μs
BLAS.set_num_threads(1)
@benchmark update_logl!($m1, false, false) setup = init_ls!($m1)
# median time:      98.884 μs

Profile.clear()
@profile update_logl!(m2, false, false)
Profile.print(format = :flat)

# # This is the bottle neck
# @profile for (i, id) in enumerate(subset_id)
#     obsvec[i] = datatable_cols |>
#         TableOperations.filter(x -> Tables.getcolumn(x, :id) == id) |>
#         Tables.columns |>
#         myobs(feformula, reformula)
# end
# Profile.print(format = :flat)


text:
# Gated-Feedback LSTM

type GFLSTMLayer # a single layer of a multi-layer RNN
        # cell's input gate params
        wix::NNMatrix
        wih::NNMatrix
        bi::NNMatrix
        # cell's forget gate parameters
        wfx::NNMatrix
        wfh::NNMatrix
        bf::NNMatrix
        # cell's out gate parameters
        wox::NNMatrix
        woh::NNMatrix
        bo::NNMatrix
        # cell's write parameters
        wcx::NNMatrix
        wgh::Array{NNMatrix, 1}
        wch::Array{NNMatrix, 1}
        wgx::Array{NNMatrix, 1}
        bg::Array{NNMatrix, 1}
        bc::NNMatrix
    function GFLSTMLayer(prevsize::Int, hiddensizes::Array{Int, 1}, layer::Int, std::Float64)
        ###  gate parameters ###
        # cell's input gate params
        hiddensize = hiddensizes[layer]
        totalsize = sum(hiddensizes)
        wix = randNNMat(hiddensize, prevsize, std)
        wih = randNNMat(hiddensize, hiddensize, std)
        bi = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        # cell's forget gate parameters
        wfx = randNNMat(hiddensize, prevsize, std)
        wfh = randNNMat(hiddensize, hiddensize, std)
        bf = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        # cell's out gate parameters
        wox = randNNMat(hiddensize, prevsize, std)
        woh = randNNMat(hiddensize, hiddensize, std)
        bo = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        ### cell's write parameters ###
        wcx = randNNMat(hiddensize, prevsize, std)
        layers = length(hiddensizes)
        wgx = Array(NNMatrix, layers)
        wgh = Array(NNMatrix, layers)
        wch = Array(NNMatrix, layers)
        bg = Array(NNMatrix, layers)
        for d in 1:layers
            wgh[d] = randNNMat(hiddensize, totalsize, std)
            wch[d] = randNNMat(hiddensize, hiddensize, std)
            wgx[d] = randNNMat(hiddensize, prevsize, std)
            bg[d] = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))
        end
        bc = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))

        new(wix,wih,bi,　wfx,wfh,bf,　wox,woh,bo,　wcx,wgh,wch,wgx,bg,bc)
    end
end

type GFLSTM <: Model
    hdlayers::Array{GFLSTMLayer,1} # holds variable number of hidden layers
    whd::NNMatrix # hidden to decoder weights & gradients
    bd::NNMatrix  # bias of hidden to decoder layer
    matrices::Array{NNMatrix,1} # used by solver - holds references to each of matrices in model
    hiddensizes::Array{Int, 1}
    function GFLSTM(inputsize::Int, hiddensizes::Array{Int, 1}, outputsize::Int, std::Float64=0.08)
        hdlayers = Array(GFLSTMLayer, length(hiddensizes))
        matrices = Array(NNMatrix, 0)
        for d in 1:length(hiddensizes)
            prevsize = d == 1? inputsize : hiddensizes[d-1]
            layer = GFLSTMLayer(prevsize, hiddensizes, d, std)
            hdlayers[d] = layer
            # input gate
            push!(matrices, layer.wix)
            push!(matrices, layer.wih)
            push!(matrices, layer.bi)
            # forget gate
            push!(matrices, layer.wfx)
            push!(matrices, layer.wfh)
            push!(matrices, layer.bf)
            # output gate
            push!(matrices, layer.wox)
            push!(matrices, layer.woh)
            push!(matrices, layer.bo)
            # cell params
            push!(matrices, layer.wcx)
            for d2 in 1:length(hiddensizes)
                push!(matrices, layer.wgh[d2])
                push!(matrices, layer.wch[d2])
                push!(matrices, layer.wgx[d2])
                push!(matrices, layer.bg[d2])
            end
            push!(matrices, layer.bc)
        end
        whd = randNNMat(outputsize, hiddensizes[end], std)
        bd = NNMatrix(outputsize, 1, zeros(outputsize,1), zeros(outputsize,1))
        push!(matrices, whd) # bias of hidden to hidden
        push!(matrices, bd) # bias of hidden to hidden
#         println("bh.w size($(size(bd.w)))")
        new(hdlayers, whd, bd, matrices, hiddensizes)
    end
end

function forwardprop(g::Graph, model::GFLSTM, x, prev)

    # forward prop for a single tick of Gated Feedback LSTM
    # g is graph to append ops to
    # model contains GFLSTM parameters
    # x is 1D column vector with observation
    # prev is a struct containing hidden and cell from previous iteration

    prevhd, prevcell, _ = prev # extract previous hidden and cell states from the tuple
    hiddenprevs = Array(NNMatrix,0)
    cellprevs = Array(NNMatrix,0)
    if length(prevhd) == 0
        for hdsize in model.hiddensizes
            push!(hiddenprevs, NNMatrix(hdsize,1))
            push!(cellprevs, NNMatrix(hdsize,1))
        end
    else
        hiddenprevs = prevhd
        cellprevs = prevcell
    end

    hidden = Array(NNMatrix,0)
    cell = Array(NNMatrix,0)
    hstar = concat(g, hiddenprevs...)
    layers = length(model.hiddensizes)
    for d in 1:layers
        input = d == 1 ? x : hidden[d-1]
        hdprev = hiddenprevs[d]
        cellprev = cellprevs[d]
        # cell's input gate params
        wix = model.hdlayers[d].wix
        wih = model.hdlayers[d].wih
        bi = model.hdlayers[d].bi
        # cell's forget gate parameters
        wfx = model.hdlayers[d].wfx
        wfh = model.hdlayers[d].wfh
        bf = model.hdlayers[d].bf
        # cell's out gate parameters
        wox = model.hdlayers[d].wox
        woh = model.hdlayers[d].woh
        bo = model.hdlayers[d].bo
        # cell's write parameters
        wcx = model.hdlayers[d].wcx
        wgh = model.hdlayers[d].wgh
        wgx = model.hdlayers[d].wgx
        wch = model.hdlayers[d].wch
        bg = model.hdlayers[d].bg
        bc = model.hdlayers[d].bc

        # input gate
        h0 = mul(g, wix, input)
        h1 = mul(g, wih, hdprev)
        inputgate = sigmoid(g, add(g, h0, h1, bi))

        # forget gate
        h2 = mul(g, wfx, input)
        h3 = mul(g, wfh, hdprev)
        forgetgate = sigmoid(g, add(g, h2, h3, bf))

        # output gate
        h4 =mul(g, wox, input)
        h5 = mul(g, woh, hdprev)
        outputgate = sigmoid(g, add(g, h4, h5, bo))

        # global reset gates
        gr = Array(NNMatrix, layers)
        @inbounds for gd in 1:layers
            hg = mul(g, wgx[gd], input)
            hu = mul(g, wgh[gd], hstar)
            gr[gd] = sigmoid(g, add(g, hg, hu, bg[gd]))
        end

        # write operation on cells
        h6 = mul(g, wcx, input)
        h = Array(NNMatrix, layers)
        @inbounds for hd in 1:layers
            hi = mul(g, wch[hd], hiddenprevs[hd])
            h[hd] = eltmul(g, gr[hd], hi)
        end
        cellwrite = tanh(g, add(g, h6, bc, h...))

        # compute new cell activation
        retaincell = eltmul(g, forgetgate, cellprev) # what do we keep from cell
        writecell = eltmul(g, inputgate, cellwrite)  # what do we write to cell
        cell_d = add(g, retaincell, writecell)        # new cell contents

        # compute hidden state as gated, saturated cell activations
        hidden_d = eltmul(g, outputgate, tanh(g, cell_d))

        push!(hidden,hidden_d)
        push!(cell, cell_d)
    end

    # one decoder to outputs at end
    output = add(g, mul(g, model.whd, hidden[end]),model.bd)

    # return cell memory, hidden representation and output
    return hidden, cell, output
end

max_issues_repo_name: Andy-P/RecurrentNN

max_forks_repo_name: Andy-P/RecurrentNN


text:
include("runtests-prerequisite.jl")

mutable struct MyNewStruct <: PostgresORM.IEntity

  actor_id::Union{Missing,Int32}
  first_name::Union{Missing,String}

  MyNewStruct(args::NamedTuple) = MyNewStruct(;args...)
  MyNewStruct(;
        actor_id = missing,
        first_name = missing,
 ) = (
            x = new(missing, missing);
            x.actor_id = actor_id;
            x.first_name = first_name;
            return x
            )

end

# ############################################### #
# TEST 1: the ORM has its function 'get_id_props()' #
# ############################################### #
module MyNewStructORM
    data_type = Main.MyNewStruct
    Main.PostgresORM.get_orm(x::Main.MyNewStruct) = return(MyNewStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :actor_id => "actor_id",
         :first_name => "first_name"
      )
    get_id_props() = return [:actor_id,:first_name]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_get_ids_props_names(o::IEntity)` with a ORM with a function `get_id_props()``" begin
    myNewStruct = MyNewStruct(actor_id = 3,first_name="toto")
    PostgresORM.get_orm(myNewStruct)
    @test PostgresORM.Controller.util_get_ids_props_names(myNewStruct) == [:actor_id,:first_name]
end

# ####################################################### #
# TEST 2: The ORM doesn't have a function 'get_id_props()'
#            but has a id_property defined as a vector
# ####################################################### #
module MyNewStructORM
    data_type = Main.MyNewStruct
    Main.PostgresORM.get_orm(x::Main.MyNewStruct) = return(MyNewStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :actor_id => "actor_id",
         :first_name => "first_name"
      )
    id_property = [:actor_id,:first_name]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_get_ids_props_names(o::IEntity)` with a ORM without a function `get_id_props()``" begin
    myNewStruct = MyNewStruct(actor_id = 3,first_name="toto")
    PostgresORM.get_orm(myNewStruct)
    @test PostgresORM.Controller.util_get_ids_props_names(myNewStruct) == [:actor_id,:first_name]
end

# ####################################################### #
# TEST 3: The ORM doesn't have a function 'get_id_props()'
#            but has a id_property defined as a scalar
# ####################################################### #
module MyNewStructORM
    data_type = Main.MyNewStruct
    Main.PostgresORM.get_orm(x::Main.MyNewStruct) = return(MyNewStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :actor_id => "actor_id",
         :first_name => "first_name"
      )
    id_property = :actor_id

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_get_ids_props_names(o::IEntity)` with a ORM without a function `get_id_props()``" begin
    myNewStruct = MyNewStruct(actor_id = 3,first_name="toto")
    orm_module = PostgresORM.get_orm(myNewStruct)
    @test PostgresORM.Controller.util_get_ids_props_names(myNewStruct) == [:actor_id]
end


text:
struct Node{NBranches,T,Values<:Union{Nothing,MutableNTuple{NBranches,T}}}
    count::Threads.Atomic{Int}
    sense::Threads.Atomic{Bool}
    winner::NTuple{2,Threads.Atomic{Int}}
    values::NTuple{2,Values}
    waiters::NTuple{2,NTuple{NBranches,OneWayCondition}}
end
# TODO: Clarify why two sets of `values` are required: different rounds can have
# different sets of tasks involved. If not, remove it.

function Node{NBranches,T}() where {NBranches,T}
    if Base.issingletontype(T)
        Values = Nothing
        values = (nothing, nothing)
    else
        Values = MutableNTuple{NBranches,T}
        values = (Values(), Values())
    end
    return Node{NBranches,T,Values}(
        Threads.Atomic{Int}(0),
        Threads.Atomic{Bool}(true),
        (Threads.Atomic{Int}(0), Threads.Atomic{Int}(0)),
        values,
        ntuple(_ -> ntuple(_ -> OneWayCondition(), NBranches), 2),
    )
end

struct TreeBarrier{NBranches,T,Op,Values} <: SyncBarriers.TreeBarrier{NBranches,T}
    n::Int
    local_sense::Vector{Bool}  # TODO: pad
    nodes::Vector{Node{NBranches,T,Values}}
    values::MutableNTuple{2,T}
    op::Op
end

TreeBarrier{NBranches}(n) where {NBranches} = TreeBarrier{NBranches,Nothing}(right, n)

function TreeBarrier{NBranches,T}(op::Op, n::Integer) where {NBranches,T,Op}
    @argcheck NBranches isa Integer && NBranches > 1 && T isa Type && n > 0
    local_sense = [true for _ in 1:n]
    nnodes = foldl_leaf_to_root((_, x) -> x.stop, nothing, n, Val(NBranches), 1)
    nodes = [Node{NBranches,T}() for _ in 1:nnodes]
    values = MutableNTuple{2,T}()
    return TreeBarrier(n, local_sense, nodes, values, op)::TreeBarrier{NBranches,T,Op}
end

function foldl_leaf_to_root(
    rf,
    acc,
    n::Integer,
    ::Val{NBranches},
    i::Integer,
) where {NBranches}
    offset = 0
    width = n  # number of nodes for each depth
    while true
        width, nlast = divrem(width, NBranches)
        width += nlast > 0
        q, r = divrem(i - 1, NBranches)
        i = q + 1
        if nlast > 0 && i == width
            branches = 1:nlast
        else
            branches = 1:NBranches
        end
        acc = rf(
            acc,
            (
                inode = offset + i,       # node = nodes[inode]
                iself = r + 1,            # node.values[_][iself] is my slot
                branches = branches,      # node.values[_][brancehs] is me + siblings
                offset = offset,
                stop = offset + width,
                width = width,
            ),
        )
        acc isa Break && break
        width == 1 && break
        offset += width
    end
    return acc
end

function recurse_leaf_to_root(
    rf,
    acc,
    n::Integer,
    ::Val{NBranches},
    i::Integer,
) where {NBranches}
    function rec(acc, offset, width)
        width, nlast = divrem(width, NBranches)
        width += nlast > 0
        q, r = divrem(i - 1, NBranches)
        i = q + 1
        if nlast > 0 && i == width
            branches = 1:nlast
        else
            branches = 1:NBranches
        end
        isroot = width == 1
        x = (
            inode = offset + i,       # node = nodes[inode]
            iself = r + 1,            # node.values[_][iself] is my slot
            branches = branches,      # node.values[_][brancehs] is me + siblings
            offset = offset,
            stop = offset + width,
            isroot = isroot,
        )
        rf(acc, x) do acc
            acc isa Break && return acc
            isroot && return acc
            rec(acc, offset + width, width)
        end
    end
    return rec(acc, 0, n)
end

function _reduce_arrive!(
    handle::BarrierHandle{<:TreeBarrier{NBranches,T}},
    value,
    ::Val{ShouldWait},
    spin = nothing,
) where {NBranches,T,ShouldWait}
    value = convert(T, value)
    barrier = handle.barrier
    i = handle.i
    s = barrier.local_sense[handle.i] = !barrier.local_sense[handle.i]
    return recurse_leaf_to_root(
        Some(value),
        barrier.n,
        Val(NBranches),
        i,
    ) do recurse, acc, x
        node = barrier.nodes[x.inode]
        vals = node.values[s+1]
        if vals !== nothing
            vals[x.iself] = something(acc)
        end
        if Threads.atomic_add!(node.count, 1) == length(x.branches) - 1
            if !ShouldWait
                node.winner[s+1][] = x.iself
            end
            if vals !== nothing
                a = vals[1]
                for j in x.branches[2:end]
                    a = barrier.op(a, vals[j])
                end
                if x.isroot
                    barrier.values[s+1] = a
                end
                acc = Some(a)
            end
            acc0 = acc
            acc = recurse(acc)
            if acc isa Break
                return acc
            end
            if !ShouldWait
                if x.isroot
                    for node in barrier.nodes
                        node.winner[!s+1][] = 0
                    end
                end
            end
            node.count[] = 0
            node.sense[] = s
            for (j, waiter) in pairs(node.waiters[s+1])
                if j != x.iself
                    notify(waiter)
                end
            end
            return acc
        else
            if ShouldWait
                sense = node.sense
                waitif(() -> sense[] != s, node.waiters[s+1][x.iself], spin)
                return nothing
            else
                return Break(nothing)
            end
        end
    end
end

function SyncBarriers.reduce!(
    handle::BarrierHandle{<:TreeBarrier{NBranches,T}},
    value,
    spin::Union{Nothing,Integer} = nothing,
) where {NBranches,T}
    acc = _reduce_arrive!(handle, value, Val(true), spin)::Union{Nothing,Some}
    barrier = handle.barrier
    s = barrier.local_sense[handle.i]
    if acc isa Some
        @_assert barrier.values[s+1] === something(acc)
        return something(acc)
    else
        return barrier.values[s+1]
    end
end

function SyncBarriers.reduce_arrive!(
    handle::BarrierHandle{<:TreeBarrier{NBranches,T}},
    value,
) where {NBranches,T}
    acc = _reduce_arrive!(handle, value, Val(false))
    if acc isa Break
        return nothing
    else
        return acc::Some
    end
end

function SyncBarriers.depart!(
    handle::BarrierHandle{<:TreeBarrier{NBranches}},
    spin::Union{Integer,Nothing} = nothing,
) where {NBranches}
    barrier = handle.barrier
    s = barrier.local_sense[handle.i]
    recurse_leaf_to_root(nothing, barrier.n, Val(NBranches), handle.i) do recurse, _, x
        node = barrier.nodes[x.inode]
        if node.winner[s+1][] == x.iself
            if x.isroot  # already notified during arrive
                return Break(nothing)
            end
            acc = recurse(nothing)
            acc isa Break && return acc

            node.count[] = 0
            node.sense[] = s
            for (j, waiter) in pairs(node.waiters[s+1])
                if j != x.iself
                    notify(waiter)
                end
            end
        else
            sense = node.sense
            waitif(() -> sense[] != s, node.waiters[s+1][x.iself], spin)
        end
        return nothing
    end
    return barrier.values[s+1]
end


text:
parse_datetime(str) = DateTime(split(str, ".")[1], "yyyy-mm-ddTHH:MM:SS.")
tostring(dt::DateTime) = Dates.format(dt, "yyyy-mm-ddTHH:MM:SS.")
clean_json(response) = JSON.parse(replace(readstring(response), r"<\/?[^>]+>|( )"," "))
function writeall(file, str)
    open(file, "w") do io
        write(io, str)
    end
end


text:
"""
    run(env, cicfg::CIService; kwargs...)
"""
function run(env = ENV,
             cicfg::CIService=auto_detect_ci_service(;env=env);
             merge_new_packages::Bool,
             merge_new_versions::Bool,
             new_package_waiting_period,
             new_jll_package_waiting_period,
             new_version_waiting_period,
             new_jll_version_waiting_period,
             registry::String,
             #
             tagbot_enabled::Bool=false,
             #
             authorized_authors::Vector{String},
             authorized_authors_special_jll_exceptions::Vector{String},
             #
             additional_statuses::AbstractVector{<:AbstractString} = String[],
             additional_check_runs::AbstractVector{<:AbstractString} = String[],
             #
             error_exit_if_automerge_not_applicable::Bool = false,
             #
             master_branch::String = "master",
             master_branch_is_default_branch::Bool = true,
             suggest_onepointzero::Bool = true,
             #
             registry_deps::Vector{<:AbstractString} = String[],
             api_url::String="https://api.github.com",
             check_license::Bool=false,
             # A list of public Julia registries (repository URLs)
             # which will be checked for UUID collisions in order to
             # mitigate the dependency confusion vulnerability. See
             # the `dependency_confusion.jl` file for details.
             public_registries::Vector{<:AbstractString} = String[],
             read_only::Bool=false)::Nothing
    all_statuses = deepcopy(additional_statuses)
    all_check_runs = deepcopy(additional_check_runs)
    push!(all_statuses, "automerge/decision")
    unique!(all_statuses)
    unique!(all_check_runs)
    api = GitHub.GitHubWebAPI(HTTP.URI(api_url))

    registry_head = directory_of_cloned_registry(cicfg; env=env)

    # Figure out what type of build this is
    run_pr_build = conditions_met_for_pr_build(cicfg; env=env, master_branch=master_branch)
    run_merge_build = conditions_met_for_merge_build(cicfg; env=env, master_branch=master_branch)

    if !(run_pr_build || run_merge_build)
        throw_not_automerge_applicable(
            AutoMergeWrongBuildType,
            "Build not determined to be either a PR build or a merge build. Exiting.";
            error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable,
        )
        return nothing
    end

    # Authentication
    key = run_pr_build || !tagbot_enabled ? "AUTOMERGE_GITHUB_TOKEN" : "AUTOMERGE_TAGBOT_TOKEN"
    auth = my_retry(() -> GitHub.authenticate(api, env[key]))
    whoami = my_retry(() -> username(api, cicfg; auth=auth))
    @info("Authenticated to GitHub as \"$(whoami)\"")
    registry_repo = my_retry(() -> GitHub.repo(api, registry; auth=auth))

    if run_pr_build
        pr_number = pull_request_number(cicfg; env=env)
        pr_head_commit_sha = current_pr_head_commit_sha(cicfg; env=env)
        pull_request_build(api,
                           pr_number,
                           pr_head_commit_sha,
                           registry_repo,
                           registry_head;
                           auth = auth,
                           authorized_authors = authorized_authors,
                           authorized_authors_special_jll_exceptions = authorized_authors_special_jll_exceptions,
                           error_exit_if_automerge_not_applicable = error_exit_if_automerge_not_applicable,
                           master_branch = master_branch,
                           master_branch_is_default_branch = master_branch_is_default_branch,
                           suggest_onepointzero = suggest_onepointzero,
                           whoami = whoami,
                           registry_deps = registry_deps,
                           check_license = check_license,
                           public_registries = public_registries,
                           read_only = read_only)
    else
        always_assert(run_merge_build)
        cron_or_api_build(api,
                          registry_repo;
                          auth = auth,
                          authorized_authors = authorized_authors,
                          authorized_authors_special_jll_exceptions = authorized_authors_special_jll_exceptions,
                          merge_new_packages = merge_new_packages,
                          merge_new_versions = merge_new_versions,
                          new_package_waiting_period = new_package_waiting_period,
                          new_jll_package_waiting_period = new_jll_package_waiting_period,
                          new_version_waiting_period = new_version_waiting_period,
                          new_jll_version_waiting_period = new_jll_version_waiting_period,
                          whoami = whoami,
                          all_statuses = all_statuses,
                          all_check_runs = all_check_runs,
                          read_only = read_only)
    end
    return nothing
end


text:
import Base: IteratorSize, HasLength, IsInfinite, length

struct FiniteTimeTrigger <: AbstractFiniteTrigger
    time_fire_at::Dates.Time
    n::Int
end

struct InfiniteTimeTrigger <: AbstractInfiniteTrigger
    time_fire_at::Dates.Time
end

"""
    TimeTrigger(t::Dates.Time[, n=number_of_times])

A trigger which should trigger a job daily at a given time.

# Optional parameter
- `n=1`: trigger once
- `n=-1` (default): trigger every day indefinitely
- `n=value`: trigger just a number of times
"""
function TimeTrigger(time_fire_at; n=-1)
    if n < 0
        InfiniteTimeTrigger(time_fire_at)
    else
        FiniteTimeTrigger(time_fire_at, n)
    end
end

"""
    Trigger(t::Dates.Time[, n=number_of_times])

Return an `TimeTrigger` which should trigger a job daily at a given time (once, a finite number of times or indefinitely).
"""
Trigger(t::Dates.Time; kwargs...) = TimeTrigger(t; kwargs...)

function get_next_dt_fire(trigger::Union{FiniteTimeTrigger,InfiniteTimeTrigger}, dt_previous_fire, dt_now)
    d = Date(dt_now)
    t = Dates.Time(dt_now)
    if trigger.time_fire_at > t
        d + trigger.time_fire_at
    else
        d + trigger.time_fire_at + Dates.Day(1)
    end
end


text:
export Scene

"""

    Scene

The `Scene` type contains the root node of the object hierarchy, a list of the light
sources in the scene, and the global bounding box.

"""
struct Scene{T<:Aggregate, S<:Spectrum}
    aggregate::T
    lights::Array{LightSource,1}
    bounds::BoundingBox
    background::S
end

Scene(agg::Aggregate, lights::Array{LightSource,1}) = Scene(agg, lights, world_bounds(agg), nolight)
Scene(agg::Aggregate, lights::Array{LightSource,1}, background::Spectrum) = Scene(agg, lights, world_bounds(agg), background)

function Scene(primitives::Array{GeometricPrimitive}, lights::Array{LightSource,1})
    return Scene(BVHAccelerator(primitives), lights)
end

function Scene(primitives::Array{GeometricPrimitive}, lights::Array{LightSource,1}, background::Spectrum)
    return Scene(BVHAccelerator(primitives), lights, background)
end

"""
    world_bounds(S::Scene)

Get the bounding box of the entire scene (this just asks for the bounding box of the
scene's acceleration aggregate.)

"""
world_bounds(S::Scene) = scene.bounds

"""
    intersect(r::Ray, S::Scene)

Find intersection of a ray in the scene. Just delegates the call to the scene's
acceleration aggregate's `intersect` method. Returns an `Intersection`.

"""
intersect(r::Ray, S::Scene) = intersect(r, S.aggregate)

"""
    intersectP(r::Ray, S::Scene)

A quicker true/false ray intersect function for a scene. This just delegates the call to
the scene's acceleration aggregate's `intersectP`.

"""
intersectP(r::Ray, S::Scene) = intersectP(r, S.aggregate)

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using UnicodePlots

scatterplot(randn(50), randn(50),
    title = "My Scatterplot")

myPlot = densityplot([1:100], randn(100), border=:dotted)

lineplot(rand(3), rand(3), title = "My Lineplot")

max_stars_repo_path: Chapter 5/ch05-unicode.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 5/ch05-unicode.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 5/ch05-unicode.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
#this file holds functions related to data processing

#reads in a CSV type
function readSingleTypeCSV(path::String; colType=Union{Missing, String},
    streamer::Function = I_CSV_STREAMER)::DataFrame

  #count the number of cols
  iStream = streamer(path)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)
  nCols::Int = size(df,2)

  #set the type
  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  df = CSV.read(iStream, strings=:raw, categorical=false, types=types)
  close(iStream)

  return df
end

#same as the read version only returns nothing and just checks for errors
function validateSingleTypeCSV(path::String; colType=Union{Missing, String},
  streamer::Function = I_CSV_STREAMER)::Nothing

  #count the number of cols

  iStream = streamer(path)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)
  nCols::Int = size(df,2)

  #set the type

  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  CSV.validate(iStream, strings=:raw, categorical=false, types=types)
  close(iStream)

end

#this fucntion is almsot a script which analyzes data
function preProcessData(; refreshData::Bool=true, validateDownloads::Bool = true,
  nccsPath::String = NCCS_PATH,  refreshMetadata::Bool = true, reserializeData::Bool = true)::Nothing

  #downloadNCCSCSVs(NCCS_PATH, refreshData=refreshData, #=delay=0.001,=# validate=validateDownloads)
  metadata::Vector{MetaDictionary} = getMetadata(refreshMetadata = refreshMetadata)
  #=loadNCCSFiles(metadata, reserializeData=reserializeData, serializeOnly = true, mode=:parallel)
  summarizeNCCSFile(testFile)=#
  #writeSummaryNCCSTables(metadata[25:27])

  #indexDF::DataFrame = readMetaIndex()
  #metaDict::MetaDictionary = MetaDictionary(indexDF[1,:name])
  #println("fileName: $(metaDict.fileName)")
  future = @spawnat workers()[1] willFail(metadata[1])

  fetch(future)

  #=println("fileName: $(metadata[1].fileName)")
  willFail(metadata[1])=#

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function AR_exit_T(MeName::String, MeOName::String, MeOFrac::Real, NMe::Real, NMeO::Real,
    InertName::String, Ts_in::Real, Tg_in = 300::Real)
    if MeOFrac > 1
        println("The metal fraction ($MeOFrac) must be smaller than 1! ")
    else
        # mole of metal and inert
        M_Me = 1000*MeOFrac/Mm(MeOName)*NMe/NMeO
        M_inert = 1000*(1 - MeOFrac)/Mm(InertName)
        # calculate needed air
        M_O2 = M_Me/NMe/0.95 # consider 95% O2 conversion
        M_N2 = M_O2*0.79/0.21
        #
        Hinlet = M_O2*GasEnthalpy("O2", Tg_in) + M_N2*GasEnthalpy("N2", Tg_in) +
        M_Me*SolidsEnthalpy(MeName, Ts_in) + M_inert*SolidsEnthalpy(InertName, Ts_in)
        T_exit = Ts_in # initial exit temperature
        ΔT = 1.0
        while abs(ΔT) > 0.1
            dH = Hinlet - (M_Me/NMe*NMeO*SolidsEnthalpy(MeOName, T_exit) +
            M_inert*SolidsEnthalpy(InertName, T_exit) +
            M_O2*0.05*GasEnthalpy("O2", T_exit) +
            M_N2*GasEnthalpy("N2", T_exit))
            if isnan(dH)
                println("error during calculation")
                return NaN
            end
            ΔT = dH/(M_Me/NMe*NMeO*SolidsCp(MeOName, T_exit) +
            M_inert*SolidsCp(InertName, T_exit) +
            M_O2*0.05*GasCp("O2", T_exit) +
            M_N2*GasCp("N2", T_exit) )

            T_exit = T_exit + ΔT
        end
        return T_exit
    end
end


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

JLLWrappers.@generate_main_file_header("OpenSSH")
JLLWrappers.@generate_main_file("OpenSSH", UUID("9bd350c2-7e96-507f-8002-3f2e150b4e1b"))
end  # module OpenSSH_jll


text:
export inst_codes, inst_code

@doc """
    inst_codes(S::GphysData)

Get the instrument code of each channel in `S`.

    inst_code(S::GphysData, i::Integer)

Get the instrument code of channel `i`.

inst_code(C::GphysChannel)

Get the instrument code from `C.id`.

Assumes each ID ends with an alphanumeric three-digit channel code and that
ID fields are separated by periods: for example, two channels with IDs
"XX.YYY.00.EHZ" and "_.YHY" each have an instrument code of 'H'.

Channel codes less than two characters long (e.g. "Z" in "AA.BBB.CC.Z") are ignored.

SEED channel codes of seismic and seismoacoustic data (for which operations
like `detrend!` and `taper!` are sane) include D, G, H, J, L, M, N, P, Z.

SEED channel codes of seismometers (for which `translate_resp!` and
`remove_resp!` are sane) are H, J, L, M, N, P, Z.
""" inst_codes
function inst_codes(S::GphysData)
  N = S.n
  codes = Array{Char, 1}(undef, N)
  fill!(codes, '\0')
  @inbounds for i = 1:N
    id = S.id[i]
    L = length(id)
    for j = L:-1:1
      if id[j] == '.'
        j > L-2 && break
        setindex!(codes, id[j+2], i)
        break
      end
    end
  end
  return codes
end

@doc (@doc inst_codes)
function inst_code(S::GphysData, i::Integer)
  N = S.n
  id = S.id[i]
  L = length(id)
  c = '\0'
  for j = L:-1:1
    if id[j] == '.'
      j > L-2 && break
      c = id[j+2]
      break
    end
  end
  return c
end

@doc (@doc inst_codes)
function inst_code(C::GphysChannel)
  L = length(C.id)
  c = '\0'
  for j = L:-1:1
    if C.id[j] == '.'
      j > L-2 && break
      c = C.id[j+2]
      break
    end
  end
  return c
end


text:
#///////////////////////////////////////
#// File Name: predictor.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/28
#// Description: Ado prediction models for Risk Sensitive Stochastic SAC
#///////////////////////////////////////

using LinearAlgebra
using PyCall
using Distributions

abstract type Predictor end;

# Trajectron Predictor
struct TrajectronPredictorParameter <: Parameter
    prediction_steps::Int64    # Prediction steps for ado robots (Trajectron)
    num_samples::Int64         # Number of samples for ado simulation (Trajectron)
    use_robot_future::Bool     # Whether Trajectron is conditioned on the robot future.
    deterministic::Bool        # Use mode-mode sampling for deterministic prediction.
    rng_seed_py::Int64         # rng seed (Trajectron)
end

mutable struct TrajectronPredictor <: Predictor
    param::TrajectronPredictorParameter
    trajectron::PyObject       # Trajectron model
    edge_addition_filter::PyObject
    edge_removal_filter::PyObject
end

function TrajectronPredictor(param::TrajectronPredictorParameter,
                             model_dir::String,
                             conf_file::String,
                             device::PyObject;
                             verbose=true)
    if !occursin("rob", model_dir)
        incl_robot_node = false;
    else
        incl_robot_node = true;
    end
    if param.deterministic
        @assert param.num_samples == 1 "num_samples has to be 1 for deterministic prediction!"
    end

    # Trajectron arguments
    py"""
    args = dict({
        # Model Parameters
        "offline_scene_graph": "no",                    # whether to precompute the scene graphs offline, options are "no" and "yes"
        "dynamic_edges": "yes",                         # whether to use dynamic edges or not, options are "no" and "yes"
        "edge_state_combine_method": "sum",             # the method to use for combining edges of the same type
        "edge_influence_combine_method": "attention",   # the method to use for combining edge influences
        "edge_addition_filter": [0.25, 0.5, 0.75, 1.0], # what scaling to use for edges as they're created
        "edge_removal_filter": [1.0, 0.0],              # what scaling to use for edges as they're removed
        "incl_robot_node": $(incl_robot_node),          # whether to include a robot node in the graph or simply model all agents
        "map_encoding": False,                          # whether to use map encoding or not
        "no_edge_encoding": False,                      # whether to use neighbors edge encoding
        "eval_device": $(device),                       # what device to use during evaluation
        "seed": $(param.rng_seed_py),                   # manual seed to use
    })
    """
    trajectron, hyperparams, log =
        pycall(py"initialize_model", NTuple{3, PyObject},
               pycall(py"EasyDict", PyObject, py"args"), model_dir, conf_file);

    if verbose
        for string in log
            println(string)
        end
    end

    return TrajectronPredictor(param, trajectron,
                               pycall(hyperparams.get, PyObject, "edge_addition_filter"),
                               pycall(hyperparams.get, PyObject, "edge_removal_filter"))
end

function initialize_scene_graph!(predictor::TrajectronPredictor,
                                 scene_loader::TrajectronSceneLoader)

    if predictor.param.use_robot_future
        @assert scene_loader.param.incl_robot_node "Using robot future without spawning robot in the scene!"
    end
    if scene_loader.param.incl_robot_node
        @assert predictor.param.use_robot_future "Spawning robot in the scene without using robot future!"
    end
    scene_loader.online_scene.calculate_scene_graph(attention_radius=scene_loader.eval_env.attention_radius,
                                                    edge_addition_filter=predictor.edge_addition_filter,
                                                    edge_removal_filter=predictor.edge_removal_filter);
    online_env = py"Environment"(node_type_list=scene_loader.eval_env.node_type_list,
                                 standardization=scene_loader.eval_env.standardization,
                                 scenes=[scene_loader.online_scene],
                                 attention_radius=scene_loader.eval_env.attention_radius,
                                 robot_type=scene_loader.eval_env.robot_type)
    predictor.trajectron.set_environment(online_env, scene_loader.curr_time_idx - 1)
end

function sample_future_ado_positions!(predictor::TrajectronPredictor,
                                      ado_state_dict::Dict,
                                      robot_present_and_future::Union{Nothing, Array{Float64, 3}}=nothing)
    z_mode, gmm_mode = false, false
    if predictor.param.deterministic
        z_mode, gmm_mode = true, true # mode-mode sampling
    end
    if !isnothing(robot_present_and_future)
        @assert predictor.param.use_robot_future "Do not provide robot_present_and_future when not using it!"
    end
    if predictor.param.use_robot_future
        @assert !isnothing(robot_present_and_future) "robot_present_and_future must not be nothing!"
    end
    if !isnothing(robot_present_and_future)
        state_dict = Dict{PyObject, Vector{Float64}}();
        robot_state = robot_present_and_future[1, 1, :];
        robot_node = predictor.trajectron.env.scenes[1].robot;
        state_dict[robot_node] = robot_state;
        for key in keys(ado_state_dict)
            state_dict[key] = ado_state_dict[key]
        end
    else
        state_dict = ado_state_dict;
    end
    maps = nothing; # TODO: Implement maps if necessary.
    preds_py = predictor.trajectron.incremental_forward(state_dict, maps=maps,
                                                        prediction_horizon=predictor.param.prediction_steps,
                                                        num_samples=predictor.param.num_samples,
                                                        z_mode=z_mode,
                                                        gmm_mode=gmm_mode,
                                                        robot_present_and_future=robot_present_and_future)

    preds_py = preds_py[2]; # using samples only
    outputs_dict = Dict{String, Array{Float64, 3}}();
    for key in keys(preds_py)
        if isnothing(robot_present_and_future) && predictor.param.use_robot_future
            outputs_dict[pybuiltin("str")(key)] = preds_py[key][1, :, :, :]
        else
            outputs_dict[pybuiltin("str")(key)] =
                reshape(permutedims(preds_py[key], [2, 1, 3, 4]), :,
                        size(preds_py[key], 3), size(preds_py[key], 4))
        end
    end
    return outputs_dict
end

# Oracle Predictor (for Data Scenes Only)
struct OraclePredictorParameter <: Parameter
    prediction_steps::Int64
    dummy_pos::Float64 # dummy variable to fill in missing prediction values
end

function OraclePredictorParameter(prediction_steps::Int64)
    dummy_pos = 500.;
    OraclePredictorParameter(prediction_steps, dummy_pos)
end

mutable struct OraclePredictor <: Predictor
    param::OraclePredictorParameter
    eval_scene::PyObject
    state_def::PyObject
    curr_time_idx::Int64
    max_time_idx::Int64
end

function OraclePredictor(param::OraclePredictorParameter,
                         scene_loader::TrajectronSceneLoader)
    @assert scene_loader.curr_time_idx + param.prediction_steps <= scene_loader.max_time_idx
    return OraclePredictor(param, scene_loader.eval_scene, scene_loader.state_def,
                           scene_loader.curr_time_idx, scene_loader.max_time_idx);
end

function sample_future_ado_positions!(predictor::OraclePredictor,
                                      ado_pos_dict::Dict)
    # existing ado agents
    full_ado_list =
        predictor.eval_scene.get_nodes_clipped_at_time(collect(predictor.curr_time_idx:
                                                               predictor.curr_time_idx + predictor.param.prediction_steps),
                                                       predictor.state_def);
    # initialize outputs_dict
    outputs_dict = Dict{String, Array{Float64, 3}}();
    for key in full_ado_list
        if !in(key, keys(ado_pos_dict))
            continue;
        end
        outputs_dict[pybuiltin("str")(key)] =
            Array{Float64, 3}(undef, 1, predictor.param.prediction_steps, 2);
    end

    # sweep through pos_dict at different timesteps and fill in outputs_dict
    for ii = 1:predictor.param.prediction_steps
        input_dict = predictor.eval_scene.get_clipped_input_dict(predictor.curr_time_idx + ii,
                                                                 predictor.state_def)
        input_dict = convert_nodes_to_str(input_dict);
        ado_intersect = intersect(keys(input_dict), keys(outputs_dict))
        for key in ado_intersect
            outputs_dict[key][1, ii, 1] = input_dict[key][1];
            outputs_dict[key][1, ii, 2] = input_dict[key][2];
        end
        ado_missing = setdiff(keys(outputs_dict), ado_intersect)
        for ado in ado_missing
            outputs_dict[ado][1, ii, 1] = predictor.param.dummy_pos
            outputs_dict[ado][1, ii, 2] = predictor.param.dummy_pos
        end
    end
    predictor.curr_time_idx += 1;
    return outputs_dict
end

# Gaussian Predictor
struct GaussianPredictorParameter <: Parameter
    #dto::Float64
    prediction_steps::Int64
    num_samples::Int64
    deterministic::Bool   # Use mean estimate for deterministic prediction
    rng::AbstractRNG
end

mutable struct GaussianPredictor <: Predictor
    param::GaussianPredictorParameter
    dto::Float64
    ado_vel_dict::Dict{String, <:MvNormal}

    function GaussianPredictor(param::GaussianPredictorParameter,
                               dto::Float64,
                               ado_vel_dict::Dict{String, <:MvNormal})
        if param.deterministic
            @assert param.num_samples == 1 "num_samples has to be 1 for deterministic prediction!"
        end
        return new(param, dto, ado_vel_dict)
    end
end

function sample_future_ado_positions!(predictor::GaussianPredictor,
                                      ado_pos_dict::Dict{String, Vector{Float64}})
    @assert keys(predictor.ado_vel_dict) == keys(ado_pos_dict) ||
            collect(keys(predictor.ado_vel_dict)) == ["Any"]

    outputs_dict = Dict{String, Array{Float64, 3}}();
    for key in keys(ado_pos_dict)
        outputs_dict[key] =
            Array{Float64, 3}(undef, predictor.param.num_samples,
                              predictor.param.prediction_steps, 2);
        if collect(keys(predictor.ado_vel_dict)) == ["Any"]
            d = predictor.ado_vel_dict["Any"];
        else
            d = predictor.ado_vel_dict[key];
        end
        if predictor.param.deterministic
            for ii = 1:predictor.param.num_samples
                outputs_dict[key][ii, :, :] =
                    Array(transpose(hcat(
                        [d.μ for ii = 1:predictor.param.prediction_steps]...
                    )));
            end
        else
            for ii = 1:predictor.param.num_samples
                outputs_dict[key][ii, :, :] =
                    Array(transpose(hcat(
                        [rand(predictor.param.rng, d) for ii = 1:predictor.param.prediction_steps]...
                    )));
            end
        end
        outputs_dict[key] = cumsum(outputs_dict[key], dims=2).*
                            predictor.dto;
        outputs_dict[key][:, :, 1] .+= ado_pos_dict[key][1];
        outputs_dict[key][:, :, 2] .+= ado_pos_dict[key][2];
    end
    return outputs_dict
end


text:
module DomainDecomposition
using SparseArrays
using LinearAlgebra
using jInv.Mesh
using jInv.Utils
import jInv.LinearSolvers.AbstractSolver
import jInv.LinearSolvers.solveLinearSystem
import jInv.LinearSolvers.solveLinearSystem!
import jInv.LinearSolvers.copySolver;
using Multigrid
using Multigrid.ParallelJuliaSolver
using Distributed
using KrylovMethods

const DDIndType = UInt32;

export DomainDecompositionOperatorConstructor,DomainDecompositionParam,DomainDecompositionPreconditionerParam,getDomainDecompositionParam

mutable struct DomainDecompositionPreconditionerParam
	sub_problem_param
	i				:: Array{Int64,1}
	A_i				:: SparseMatrixCSC
	DirichletMass	:: Vector
	Ainv			:: AbstractSolver
end

mutable struct DomainDecompositionParam <: AbstractSolver
	PrecParams 			::Union{Array{DomainDecompositionPreconditionerParam,1},Array{RemoteChannel,1}}
	GlobalIndices		::Array{Array{DDIndType}}
	Mesh				::RegularMesh
	numDomains			::Array{Int64,1}
	overlap   			::Array{Int64,1}
	getIndicesOfCell	::Function
	Ainv				::AbstractSolver
	getSubDomainMass	::Function
	workers				::Array{Int64}
	constructor			::Any
	out					::Int64
	doClear				::Int
	nFac				::Int
	facTime				::Real
	nSolve				::Int
	solveTime			::Real
end
function getDomainDecompositionParam(Mesh,numDomains,overlap,getIndicesOfCell,Ainv::AbstractSolver,getMass::Function = identity)
	return DomainDecompositionParam((DomainDecompositionPreconditionerParam)[],[],Mesh,numDomains,overlap,getIndicesOfCell,Ainv,getMass,(Int64)[],getEmptyCtor(),0,0,0,0.0,0,0.0);
end
mutable struct DomainDecompositionOperatorConstructor
	problem_param	::Any
	getSubParams  	::Function
	getOperator		::Function
	getDirichletMass::Function
end

function getEmptyCtor()
	return DomainDecompositionOperatorConstructor(0,identity,identity,identity);
end

include("DDIndices.jl");
include("DDService.jl")
include("DDSerial.jl");
include("DDParallel.jl");

import Base.isempty
function isempty(p::DomainDecompositionParam)
	return isempty(p.PrecParams);
end

import jInv.LinearSolvers.copySolver;
function copySolver(s::DomainDecompositionParam)
	# copies absolutely what's necessary.
	error("TODO");
end

# import jInvUtils.clear!
function clear!(s::DomainDecompositionParam)
	for k=1:length(s.DDPreconditioners)
		clear!(s.DDPreconditioners[k]);
	end
	s.doClear = 0;
end

import jInv.LinearSolvers.setupSolver
function setupSolver(A::SparseMatrixCSC,DDparam::DomainDecompositionParam)
	return setupDDSerial(sparse(A'),DDparam);
end

import jInv.LinearSolvers.solveLinearSystem!;
function solveLinearSystem!(At,B,X,param::DomainDecompositionParam,doTranspose=0)

	if param.doClear==1
		clear!(param);
	end
	# build preconditioner
	if isempty(param)
		setupDDSerial(At,param);
	end
	flag = 0; rnorm = 0.0; iter = 0; resvec = [];
	if !isempty(B)
		if issparse(B)
			B = full(B);
		end
		if size(B,2) == 1
			B = vec(B);
		end
		if norm(B) == 0.0
			X[:] .= 0.0;
			return X, param,flag,rnorm,iter,resvec;
		end
		Prec = r->solveDDSerial(At,r,zeros(eltype(X),size(X)),param,1,doTranspose)[1];
		# Prec = r->solveGSDDSerial(At,r,zeros(eltype(X),size(X)),param,1,doTranspose)[1];
		X,flag,rnorm,iter,resvec = KrylovMethods.fgmres(getAfun(At,zeros(eltype(X),size(X)),4),B,5,tol = 1e-6,maxIter = 100,M = Prec, x = X,out=1,flexible=true);
		# solveDD(At,B,X,param,doTranspose)
	end
	return X,param,flag,rnorm,iter,resvec
end

end


text:
## --- Digitize image data

    """
    ```julia
    (x,dx,y,dy) = digitize_plotmarkers(img, marker_color, xlims, ylims; atol=0.16)
    ```
    Calculate approximate x and y positions and uncertainties for colored
    markers in an image
    """
    function digitize_plotmarkers(img, marker_color, xlims, ylims; atol=0.16)

        # Test for approximate equality in color to marker
        t = isapprox.(marker_color,img,atol=atol)

        # Figure out our image dimensions
        xrows = size(t,2)
        yrows = size(t,1)
        xlims = collect(xlims)
        ylims = collect(ylims)

        # Allocate index arrays
        imin = Array{Float64}(undef,round(Int,xrows/2))
        imax = Array{Float64}(undef,round(Int,xrows/2))
        jmin = Array{Float64}(undef,round(Int,xrows/2))
        jmax = Array{Float64}(undef,round(Int,xrows/2))

        # Fill index arrays
        found = false
        markernumber = 0
        for j=1:xrows
            if any(t[:,j])
                list = findall(t[:,j])
                if ~found
                    markernumber += 1
                    imin[markernumber] = minimum(list)
                    imax[markernumber] = maximum(list)
                    jmin[markernumber] = j
                    jmax[markernumber] = j
                else
                    imin[markernumber] = min(imin[markernumber],minimum(list))
                    imax[markernumber] = max(imax[markernumber],maximum(list))
                    jmax[markernumber] = j
                end
                found = true
            else
                found = false
            end
        end

        # Return only the filled indices
        imin = imin[1:markernumber]
        imax = imax[1:markernumber]
        jmin = jmin[1:markernumber]
        jmax = jmax[1:markernumber]

        # Calculate x and y positions from indices
        y = ylims[2] - (imin+imax)/2 * nanrange(ylims) / yrows
        dy = (imax-imin)/2 * nanrange(ylims) / yrows
        x = (jmin+jmax)/2 * nanrange(xlims) / xrows + xlims[1]
        dx = (jmax-jmin)/2 * nanrange(xlims) / xrows

        return x, dx, y, dy
    end
    export digitize_plotmarkers

    """
    ```julia
    (x,y) = digitize_plotline(img, line_color, xlims, ylims; atol=0.16)
    ```
    Calculate approximate x and y positions for a colored line in an image
    """
    function digitize_plotline(img, line_color, xlims, ylims; atol=0.16)
        # Test for approximate equality in color to marker
        t = isapprox.(line_color,img,atol=atol)

        # Figure out our image dimensions
        xrows = size(t,2)
        yrows = size(t,1)
        xlims = collect(xlims)
        ylims = collect(ylims)

        # Calculate x for each column
        x = cntr(range(xlims[1], xlims[2], length=xrows+1))

        # y as a function of i-position in image
        # (note: images are typically flipped)
        yᵢ(i) = ylims[2] - i * nanrange(ylims) / yrows

        # Calculate y, defaulting to NaN if no matches
        y = fill(NaN, xrows)
        for j = 1:xrows
            if any(t[:,j])
                list = findall(t[:,j])
                y[j] = yᵢ(nanmean(list))
            else
                y[j] = NaN
            end
        end

        return x, y
    end
    export digitize_plotline

## --- End of File


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This code depends upon sudoku_common.jl and sudoku_permute.jl

using Random

"""
    symbol_swap!(puzzle,value_1,value_2)

Within a SolvablePuzzle interchange the roles of two integers.
"""
function symbol_swap!(puzzle::SolvablePuzzle,value_1::Integer,value_2::Integer)
    rank = get_rank(puzzle)
    rank_squared = rank*rank
    p = BitVector(undef,rank_squared) # reused in loop

    for row = 1:rank_squared
        for col = 1:rank_squared
            p .= false
            # Swap possible value roles
            p.chunks[1] = puzzle.grid[row,col].possibilities
            possible_value_1 = p[value_1]
            possible_value_2 = p[value_2]
            # Change the contents if the bits are different
            if possible_value_1 != possible_value_2
                p[value_2] = possible_value_1
                p[value_1] = possible_value_2
                puzzle.grid[row,col] = PuzzleEntry(p)
            end # end switch on required change
        end # end column loop
    end # end row loop
end

"""
    random_permutation!(puzzle,rng)

Apply a random permutation to a SolvablePuzzle and return the permutation.

Permutation is reported as a variable-length sequence of integers.
"""
function random_permutation!(puzzle::SolvablePuzzle,rng::AbstractRNG)
    rank = get_rank(puzzle)
    # Symbols to select permutation type
    #undefined = 0 # less than a valid choice
    #column = 1
    #block_column = 2
    #row = 3
    #block_row = 4
    #mirror_horizontal = 5
    #mirror_vertical = 6
    #symbol_swap = 7
    #error = 8 # greater than a valid choice
    permutation_type = rand(rng,1:7)
    if permutation_type == 1 || permutation_type == 3
        block = rand(rng,1:rank)
        sub_row_col_a = rand(rng,1:rank)
        sub_row_col_b = 1 + ((sub_row_col_a - 1 + rand(rng,1:(rank-1))) % rank)
        if permutation_type == 1
            col_permute!(puzzle.grid,block,sub_row_col_a,sub_row_col_b)
        else
            row_permute!(puzzle.grid,block,sub_row_col_a,sub_row_col_b)
        end
        return [permutation_type,block,sub_row_col_a,sub_row_col_b]
    elseif permutation_type == 2 || permutation_type == 4
        source_block = rand(rng,1:rank)
        destination_block = 1 + ((source_block - 1 + rand(rng,1:(rank-1))) % rank)
        if permutation_type == 2
            col_block_permute!(puzzle.grid,source_block,destination_block)
        else
            row_block_permute!(puzzle.grid,source_block,destination_block)
        end
        return [permutation_type,source_block,destination_block]
    elseif permutation_type == 5
        mirror_horizontal!(puzzle.grid)
        return [permutation_type]
    elseif permutation_type == 6
        mirror_vertical!(puzzle.grid)
        return [permutation_type]
    elseif permutation_type == 7
        rank_squared = rank*rank
        symbol_1 = rand(rng,1:rank_squared)
        symbol_2 = 1 + ((symbol_1 - 1 + rand(rng,1:(rank_squared-1))) % rank_squared)
        symbol_swap!(puzzle,symbol_1,symbol_2)
        return [permutation_type,symbol_1,symbol_2]
    end
end

"""
    random_solution(rank,rng,permutation_count)

Given a rank, an RNG, and a number of permutations, make new solution.
"""
function random_solution(rank::Integer,rng::AbstractRNG,permutation_count::Integer)
    result = SolvablePuzzle(rank)
    for p = 1:permutation_count
       random_permutation!(result,rng)
    end
    return result
end

"""
    random_puzzle(rank,rng,permutation_count,removals)

Construct a randomized solution, puzzle pair.
"""
function random_puzzle(rank::Integer,rng::AbstractRNG,permutation_count::Integer,removals::Integer)
    rank_squared = rank*rank
    if removals > rank_squared*rank_squared
        throw(DomainError("Cannot remove more entries than rank dictates"))
    elseif rank_squared*rank_squared - removals <= 3*rank_squared
        print("WARNING: high number of removals requested for rank")
    end
    solution = random_solution(rank,rng,permutation_count)
    puzzle = deepcopy(solution)
    for i = 1:removals
        # Randomly grab rows and columns to remove
        row = rand(rng,1:rank_squared)
        col = rand(rng,1:rank_squared)
        # If we hit a duplicate removal, jump around until we find an intact entry
        while get_value(puzzle.grid[row,col]) == 0
            row = rand(rng,1:rank_squared)
            col = rand(rng,1:rank_squared)
        end
        set_unknown(puzzle,row,col)
    end
    return [solution,puzzle]
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
## Various helper functions for working with latitude and longitude
const Location2D = SVector{2, Float64}
const LatLonCoords = NamedTuple{(:lat, :lon), Tuple{Float64,Float64}}

"""
Return the vector form of a latitude-longitude point [lat, lon]
"""
function convert_to_vector(c::LatLonCoords)
    return Location2D(c.lat, c.lon)
end

struct EuclideanLatLongMetric <: Metric
end

const AgentTask = NamedTuple{(:origin, :dest), Tuple{Int64, Int64}}
const AvgSpeed = 10.0  # 6 m/s for drones and cars
const MaxHop = 3

## Types for Aerial MAPF problem

"""
MAPF State for aerial routes. Thin wrapper around GTG vertices, with just drone time.
"""
@with_kw struct AerialMAPFState <: MAPFState
    idx::Int64
    ground_transit_idx::Int64   # Which TRANSIT GRAPH vertex ID
    timeval::Float64 = 0.0  # Only relevant for transit points
    distval::Float64 = 0.0
    considered_cars::Set{Int64} = Set{Int64}()
end

@enum ActionType Fly=1 Stay=2
struct GroundTransitAction <: MAPFAction
    action::ActionType
end

@with_kw struct GroundTransitConflict <: MAPFConflict
    overlap_gids::Vector{Int64}          # Which ground IDs are overburdened
    aerial_agent_ids::Vector{Vector{Int64}}
end

@with_kw struct GroundTransitConstraint <: MAPFConstraints
    avoid_gid_vertex_set::Set{Int64}    = Set{Int64}()
end
Base.isempty(gtc::GroundTransitConstraint) = isempty(gtc.avoid_gid_vertex_set)
MultiAgentPathFinding.get_empty_constraint(::Type{GroundTransitConstraint}) = GroundTransitConstraint()

@with_kw struct GroundTransitVertexConstraint <: MAPFConstraints
    avoid_gtg_vertex_set::Set{Int64} = Set{Int64}()
end
Base.isempty(gtvc::GroundTransitVertexConstraint) = isempty(gtvc.avoid_gtg_vertex_set)
MultiAgentPathFinding.get_empty_constraint(::Type{GroundTransitVertexConstraint}) = GroundTransitVertexConstraint()

"""
For the graph of ground time-stamped waypoints. Must be mutable as cars can wait for drones
"""
@with_kw mutable struct GroundTransitState
    idx::Int64                      # Only for reverse lookup
    ground_agent_id::Int64 = 0      # When 0, indicates aerial start-goal
    vertex_along_path::Int64 = 0    # Goes from 1 to path length
    road_vtx_id::Int64              # Vertex ID in road graph
    timeval::Float64 = 0.0          # Global time of ground vehicle
end

Graphs.vertex_index(g::SimpleVListGraph{GroundTransitState}, v::GroundTransitState) = v.idx

## Types for ground MAPF
@with_kw struct GroundMAPFState <: MAPFState
    road_vtx_id::Int64
    is_aerial::Bool
end
Base.isequal(s1::GroundMAPFState, s2::GroundMAPFState) = (s1.road_vtx_id == s2.road_vtx_id && s1.is_aerial == s2.is_aerial)

@with_kw mutable struct GMAPFStateInfo
    idx::Int64
    gval::Float64
    aerial_ids::Set{Int64}
end

@with_kw struct GroundMAPFAction <: MAPFAction
    edge_id::LightGraphs.Edge
    edge_aerial_ids::Set{Int64}
end

@with_kw struct GroundMAPFConflict <: MAPFConflict
    ground_id_pair::Tuple{Int64,Int64}
    conflicting_edge_aerial_ids::Vector{Tuple{LightGraphs.Edge,Set{Int64}}}
end

@with_kw struct GroundMAPFConstraint <: MAPFConstraints
    avoid_edge_copy_set::Dict{LightGraphs.Edge,Set{Int64}} = Dict{LightGraphs.Edge,Set{Int64}}()
end
Base.isempty(gpc::GroundMAPFConstraint) = isempty(gpc.avoid_edge_copy_set)
MultiAgentPathFinding.get_empty_constraint(::Type{GroundMAPFConstraint}) = GroundMAPFConstraint()

"""
General path state type for both agent types
"""
@with_kw mutable struct AgentPathState
    road_vtx_id::Int64
    timeval::Float64
    coord_agent_id::Set{Int64} = Set{Int64}()     # Of the other type
end

"""
General type to store the paths for agents, agnostic to ground/aerial
"""
@with_kw mutable struct AgentPathInfo
    path_states::Vector{AgentPathState} = AgentPathState[]
    total_dist::Float64                 = 0.0
    total_time::Float64                 = 0.0
end

"""
Weight-discounted aerial-annotated edge copies
"""
@with_kw mutable struct EdgeCopyInfo
    aerial_id_to_weight::Dict{Int64, Float64} = Dict{Int64, Float64}()
end

"""
The overall struct for the full MAPF problem.
"""
@with_kw mutable struct CoordinatedMAPFEnv <: MAPFEnvironment
    ground_task_list::Vector{AgentTask}
    aerial_task_list::Vector{AgentTask}
    road_graph::LightGraphs.AbstractGraph
    road_graph_wts::AbstractMatrix{Float64}
    location_list::Vector{Location2D}
    alpha_weight_distance::Float64
    car_capacity::Int64
    ground_paths::Vector{AgentPathInfo} = AgentPathInfo[]
    aerial_paths::Vector{AgentPathInfo} = AgentPathInfo[]
    ground_transit_graph::SimpleVListGraph{GroundTransitState} = SimpleVListGraph{GroundTransitState}()
    aerial_mapf_graph::SimpleVListGraph{AerialMAPFState} = SimpleVListGraph{AerialMAPFState}()
    ground_transit_weights::Dict{Tuple{Int64,Int64}, Float64} = Dict{Tuple{Int64,Int64}, Float64}()
    ground_transit_expanded_paths::Dict{Tuple{Int64,Int64}, Vector{Int64}} = Dict{Tuple{Int64,Int64}, Vector{Int64}}()
    gtg_ground_path_start_ends::Vector{Tuple{Int64,Int64}} = Tuple{Int64,Int64}[]
    gtg_drone_start_goal_idxs::Vector{Tuple{Int64,Int64}} = Tuple{Int64,Int64}[]
    next_gtg_goal_idx::Int64 = 0
    next_amapfg_goal_idx::Int64=0
    num_global_conflicts::Int64 = 0
    threshold_global_conflicts::Int64   = typemax(Int64)
    gtg_idx_gval::Dict{Int64, Float64} = Dict{Int64, Float64}()
    edge_to_copy_info::Dict{LightGraphs.Edge, EdgeCopyInfo} = Dict{LightGraphs.Edge, EdgeCopyInfo}()
    ground_mapf_graph::SimpleVListGraph{GroundMAPFState} = SimpleVListGraph{GroundMAPFState}()
    unique_gmapf_states::Dict{GroundMAPFState,GMAPFStateInfo} = Dict{GroundMAPFState,Tuple{Int64,Float64}}()
    next_gmapfg_goal_idx::Int64 = 0
    gtg_idx_to_aerial_ids::Dict{Int64,Vector{Int64}} = Dict{Int64,Vector{Int64}}()
end


text:
# https://github.com/JuliaLang/julia/pull/29679
if VERSION < v"1.1.0-DEV.472"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end


text:
macro define(name, body)
  quote
    macro $(esc(name))()
      esc($(Expr(:quote, body)))
    end
  end
end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
"""
$(TYPEDEF)

Abstract type for impedance system.
"""
abstract type AbstractImpedanceSystem{Tv <: Number} end

mutable struct ImpedanceSystem{Tv} <: AbstractImpedanceSystem{Tv}
    """
    Nonzero pattern of time domain system matrix
    """
    sysnzval::AbstractVector{Complex{Tv}}

    """
    Discretization grid
    """
    grid

    """
    Derivative of storage term
    """
    storderiv::AbstractMatrix{Tv}

    """
    Complex matrix of impedance system
    """
    matrix::AbstractMatrix{Complex{Tv}}

    """
    Right hand side of impedance system
    """
    F::AbstractMatrix{Complex{Tv}}

    """
    Stationary state
    """
    U0::AbstractMatrix{Tv}

    ImpedanceSystem{Tv}() where Tv =new()
end

"""
$(SIGNATURES)

Construct impedance system from time domain system `sys` and steady state solution `U0`
under the assumption of a periodic perturbation of species `excited_spec` at  boundary `excited_bc`.

"""
function ImpedanceSystem(sys::AbstractSystem{Tv,Ti}, U0::AbstractMatrix, excited_spec,excited_bc) where {Tv,Ti}
    residual=copy(U0)

    # Ensure that sys.matrix contains the jacobian at U0
    # Moreover, here we use the fact that for large time step sizes,
    # the contribution from the time derivative (which should not belong to the
    # main part of the impedance matrix) is small. We also pass the steady state
    # value as the "old  timestep" value.
    # An advantage of this approach is the fact that this way, we get the
    # nonzero pattern for the iω term right.
    eval_and_assemble(sys,U0,U0,residual,1.0e30)

    this=ImpedanceSystem{Tv}()
    this.grid=sys.grid
    # catch the nonzero values of the system matrix
    this.sysnzval=complex(nonzeros(sys.matrix))
    m,n=size(sys.matrix)

    # initialize storage derivative with zero
    this.storderiv=spzeros(Tv,m,n)

    # create sparse matrix with the same nonzero pattern of original matrix
    this.matrix=SparseMatrixCSC(m,n,
                                colptrs(sys.matrix),
                                rowvals(sys.matrix),
                                copy(this.sysnzval)
                                )
    this.U0=U0

    # Initialize right hand side of impedance system
    this.F=unknowns(Complex{Float64},sys)

    matrix=this.matrix
    storderiv=this.storderiv
    F=this.F
    grid=sys.grid

    physics=sys.physics

    # Prepare calculation of derivative of the storage part
    data=physics.data
    node=Node{Tv,Ti}(sys)
    bnode=BNode{Tv,Ti}(sys)
    nspecies=num_species(sys)
    nodeparams=(node,)
    bnodeparams=(bnode,)
    if isdata(data)
        nodeparams=(node,data,)
        bnodeparams=(bnode,data,)
    end

    storagewrap=function(y::AbstractVector, u::AbstractVector)
        y.=0
        physics.storage(y,u,nodeparams...)
    end

    bstoragewrap=function(y::AbstractVector, u::AbstractVector)
        y.=0
        physics.bstorage(y,u,bnodeparams...)
    end

    UK=Array{Tv,1}(undef,nspecies)
    Y=Array{Tv,1}(undef,nspecies)

    F.=0.0

    # structs holding diffresults for storage
    result_s=DiffResults.DiffResult(Vector{Tv}(undef,nspecies),Matrix{Tv}(undef,nspecies,nspecies))

    geom=grid[CellGeometries][1]
    csys=grid[CoordinateSystem]
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]

    bgeom=grid[BFaceGeometries][1]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]

    # Interior cell loop for building up storage derivative
    for icell=1:num_cells(grid)
        for inode=1:num_nodes(geom)
            _fill!(node,cellnodes,cellregions,inode,icell)
            @views UK[1:nspecies]=U0[:,node.index]

            # Evaluate & differentiate storage term at U0
            ForwardDiff.jacobian!(result_s,storagewrap,Y,UK)
            jac_stor=DiffResults.jacobian(result_s)

            # Sort it into storderiv matrix.
            K=node.index
            for idof=_firstnodedof(F,K):_lastnodedof(F,K)
                ispec=_spec(F,idof,K)
                for jdof=_firstnodedof(F,K):_lastnodedof(F,K)
                    jspec=_spec(F,jdof,K)
                    _addnz(storderiv,idof,jdof,jac_stor[ispec,jspec],sys.cellnodefactors[inode,icell])
                end
            end

        end
    end

    # Boundary face loop for building up storage derivative
    # and right hand side contribution from boundary condition
    for ibface=1:num_bfaces(grid)
        ibreg=bfaceregions[ibface]
        for ibnode=1:num_nodes(bgeom)
            _fill!(bnode,bfacenodes,bfaceregions,ibnode,ibface)
            @views UK[1:nspecies]=U0[:,bnode.index]

            # Set right hand side for excited boundary conditions
            # We don't need to put the penalty term to storderiv
            # as it is already in the main part of the matrix
            if ibreg==excited_bc
                for ispec=1:nspecies
                    if ispec!=excited_spec
                        continue
                    end
                    idof=dof(F,ispec,bnode.index)
                    if idof>0
                        fac=sys.boundary_factors[ispec,ibreg]
                        if fac==Dirichlet
                            F[ispec,bnode.index]+=fac
                        else
                            F[ispec,bnode.index]+=fac*sys.bfacenodefactors[ibnode,ibface]
                        end
                    end
                end
            end

            if isdefined(physics, :bstorage) # should involve only bspecies
                # Evaluate & differentiate storage term
                ForwardDiff.jacobian!(result_s,bstoragewrap,Y,UK)
                res_bstor=DiffResults.value(result_s)
                jac_bstor=DiffResults.jacobian(result_s)

                K=bnode.index
                for idof=_firstnodedof(F,K):_lastnodedof(F,K)
                    ispec=_spec(F,idof,K)
                    for jdof=_firstnodedof(F,K):_lastnodedof(F,K)
                        jspec=_spec(F,jdof,K)
                        _addnz(storderiv,idof,jdof,jac_bstor[ispec,jspec],sys.bfacenodefactors[ibnode,ibface])
                    end
                end
            end
        end
    end
    return this
end

"""
$(TYPEDSIGNATURES)

Create a vector of unknowns of the impedance system
"""
unknowns(this::ImpedanceSystem{Tv}) where Tv=copy(this.F)

"""
$(TYPEDSIGNATURES)

Solve the impedance system for given frequency `ω`.
"""
function solve!(UZ::AbstractMatrix{Complex{Tv}},this::ImpedanceSystem{Tv}, ω) where Tv
    iω=ω*1im
    matrix=this.matrix
    storderiv=this.storderiv
    grid=this.grid
    # Reset the matrix nonzero values
    nzval=nonzeros(matrix)
    nzval.=this.sysnzval
    # Add the ω dependent term to main diagonal.
    # This makes the implicit assumption that storderiv does not
    # introduce new matrix entries.
    U0=this.U0
    for inode=1:num_nodes(grid)
        for idof=_firstnodedof(U0,inode):_lastnodedof(U0,inode)
            for jdof=_firstnodedof(U0,inode):_lastnodedof(U0,inode)
                _addnz(matrix,idof,jdof,storderiv[idof,jdof],iω)
            end
        end
    end
    # Factorize + solve
    lufact=LinearAlgebra.lu(matrix)
    ldiv!(values(UZ),lufact,values(this.F))
end

"""
    $(TYPEDSIGNATURES)

    Calculate the derivative of the scalar measurement functional at steady state U0

    Usually, this functional is  a test function integral.  Initially,
    we assume that its value depends on all unknowns of the system.

"""
function measurement_derivative(sys::AbstractSystem,measurement_functional,U0)

    # Create a sparse 1×ndof matrix assuming that the functional
    # depends on all unknowns in the system
    ndof=num_dof(sys)
    colptr=[i for i in 1:ndof+1]
    rowval=[1 for i in 1:ndof]
    nzval=[1.0 for in in 1:ndof]
    jac=SparseMatrixCSC(1,ndof,colptr,rowval,nzval)

    # See https://github.com/JuliaDiff/SparseDiffTools.jl

    # Color the matrix for automtic differentiation
    colors = matrix_colors(jac)

    # Use Julia automatic differentiation for the calculation of the Jacobian
    forwarddiff_color_jacobian!(jac, measurement_functional, values(U0), colorvec = colors)

    # Drop any zero entries
    dropzeros!(jac)

    return jac
end

function freqdomain_impedance(isys, # frequency domain system
                              ω,    # frequency
                              U0 ,  # steady state slution
                              excited_spec,excited_bc,excited_bcval,
                              dmeas_stdy, # Derivative of steady state part of measurement functional
                              dmeas_tran  # Derivative of transient part of the measurement functional
                              )

    iω=1im*ω
    # solve impedance system
    UZ=unknowns(isys)
    solve!(UZ,isys,ω)

    # obtain measurement in frequency  domain
    m_stdy=dmeas_stdy*values(UZ)
    m_tran=dmeas_tran*values(UZ)

    # Calculate complex measurement
    z=m_stdy[1]+iω*m_tran[1]

    # return impedance (in fact, impedance is the reciprocal value here, this needs
    # to be changed)
    return z
end


text:
# SPDX-License-Identifier: X11
# 2020-11-10
# Frog 2 (100pt)

const INF = 10^9 + 5

function main()
  N,K = parse.(Int, (split ∘ readline)())
  H = parse.(Int, (split ∘ readline)())

  DP = fill(INF, N)
  DP[1] = 0

  for i ∈ 2:N
    for j ∈ 1:K
      i <= j && break
      DP[i] = min(DP[i], abs(H[i - j] - H[i]) + DP[i - j])
    end
  end

  println(DP[N])
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Pkg

Pkg.add("JuMP")
Pkg.add("Test")
Pkg.add("Compat")
Pkg.add("MathOptInterface")
Pkg.add("GLPK")
Pkg.add("GLPKMathProgInterface")
Pkg.add("LinearFractional")

Pkg.add("Cbc")
Pkg.add("Clp")
Pkg.add("ECOS")
Pkg.add("Ipopt")
Pkg.add("NLopt")
Pkg.add("SCS")
# Pkg.add("CPLEX") # Needs working installation
# Pkg.add("Gurobi") # Needs working installation
Pkg.add("AmplNLWriter")
# Pkg.add("CoinOptServices")

Pkg.add("JSON")

Pkg.update()

max_stars_repo_name: rodoufu/knapsack-jl

max_issues_repo_name: rodoufu/knapsack-jl

max_forks_repo_name: rodoufu/knapsack-jl


text:
using Drawing

using Base.Test
using SHA: sha1
using Images
using Colors

include("cairo.jl")

MAX_DIFF = 8

function compare(name)
    print("$(name): ")
    im1 = load(name)
    im2 = load("target/$(name)")
    @test size(im1) == size(im2)
    for c in zip(im1, im2)
        c1, c2 = c
        d = colordiff(c1, c2)
        if d >= MAX_DIFF
            println("$c1 $c2 $(colordiff(c1, c2))")
            @test colordiff(c1, c2) < MAX_DIFF
        end
    end
    println("ok")
end

function axes()
    move(0,0)
    line(1,0)
    move(0,0)
    line(0,1)
end

include("basics.jl")
include("shapes.jl")
include("text.jl")
include("errors.jl")
include("colours.jl")
include("towers.jl")


text:
module ClimateGNSS

# Main file for the ClimateERA module that downloads and processes ECMWF
# reanalysis data.

## Modules Used
using Dates, DelimitedFiles, Printf
using NCDatasets
using GeoRegions

## Exporting the following functions:
export
        gnssroot, gnssresort, gnssanalysis, gnssfol,
        eosload, eosloadstation

## Including other files in the module

include("startup.jl")
include("resort.jl")
#include("analysis.jl")
include("eossg.jl")
#include("nglunavco.jl")

end # module


text:
function siteind(ψ::InfiniteMPS, n::Integer)
  return uniqueind(ψ[n], ψ[n - 1], ψ[n + 1])
end

siteind(ψ::InfiniteCanonicalMPS, n::Integer) = siteind(ψ.AL, n)

function default_middle_range(N::Integer, ncenter::Integer)
  return round(Int, N / 2 - ncenter / 2 + 1, RoundDown):round(
    Int, N / 2 + ncenter / 2, RoundDown
  )
end

# Get the closest infinite MPS approximation to the finite MPS
function infinitemps_approx(
  ψ::MPS;
  nsites=1,
  nrange=default_middle_range(length(ψ), nsites),
  nsweeps=10,
  outputlevel=0,
)
  N = length(ψ)

  # nsites are the number of site in the unit cell
  # of the infinite MPS
  @assert length(nrange) == nsites

  s = siteinds(ψ)

  # The site we are interested in now
  ψ = orthogonalize(ψ, first(nrange))
  χ = maximum(n -> linkdim(ψ, n), nrange)

  site_tags = ITensors.commontags(siteinds(ψ))

  # The sites of the infinite MPS
  A = InfiniteMPS(
    [settags(s[nrange[n]], addtags(site_tags, "n=$n")) for n in 1:length(nrange)]; space=χ
  )
  randn!.(A)
  ψ∞ = orthogonalize(A, :)

  # Site indices of the infinite MPS
  # corresponding to those of the finite MPS
  inf_range = (-first(nrange) + 2):(-first(nrange) + N + 1)

  s∞ = [siteind(ψ∞, n) for n in inf_range]

  ψ = replace_siteinds(ψ, s∞)

  for sweep in 1:10
    # XXX: these are just used for testing
    TL = ITensor(1)
    for n in 1:N
      TL = TL * ψ[n] * ψ∞.AL[inf_range[n]]
    end
    TR = ITensor(1)
    for n in reverse(1:N)
      TR = TR * ψ[n] * ψ∞.AR[inf_range[n]]
    end

    L0 = randomITensor(linkind(ψ∞.AL, first(inf_range) - 1 => first(inf_range)))
    vals, vecs, _ = eigsolve(transpose(ITensorMap(ψ∞.AL[inf_range], ψ)), L0)
    L0 = vecs[1]
    λL0 = vals[1]
    LTL = L0 * TL
    LTL = replaceinds(LTL, inds(LTL) => inds(L0))

    RN = randomITensor(linkind(ψ∞.AR, last(inf_range) => last(inf_range) + 1))
    vals, vecs, _ = eigsolve(ITensorMap(ψ∞.AR[inf_range], ψ), RN)
    RN = vecs[1]
    λRN = vals[1]

    TRR = TR * RN
    TRR = replaceinds(TRR, inds(TRR) => inds(RN))

    L = OffsetVector(Vector{ITensor}(undef, N + 1), -1)
    R = OffsetVector(Vector{ITensor}(undef, N + 1), -1)
    L[0] = L0
    R[N] = RN

    for n in 1:N
      L[n] = L[n - 1] * ψ∞.AL[inf_range[n]] * dag(ψ[n])
    end
    for n in reverse(0:(N - 1))
      R[n] = R[n + 1] * ψ∞.AR[inf_range[n + 1]] * dag(ψ[n + 1])
    end

    AL = Vector{ITensor}(undef, nsites)
    AR = Vector{ITensor}(undef, nsites)
    C = Vector{ITensor}(undef, nsites)
    for n′ in eachindex(nrange)
      n = nrange[n′]
      C1 = L[n - 1] * R[n - 1]
      C2 = L[n] * R[n]
      C12 = L[n - 1] * dag(ψ[n]) * R[n]
      U, P = polar(C12 * dag(C2), uniqueinds(C12, C2))
      u = commoninds(U, P)
      U = noprime(U, u)
      AL[n′] = U
      # XXX: compute AR from a right polar decomposition of C2
      AR[n′] = replacetags(U, "Left" => "Right")
      C[n′] = C2 / norm(C2)
    end
    ψL = InfiniteMPS(AL)
    ψR = InfiniteMPS(AR)
    ψC = InfiniteMPS(C)
    ψ∞ = InfiniteCanonicalMPS(ψL, ψC, ψR)
    # This is a measure of the overlap
    if outputlevel > 0
      @show abs(λL0)^(1 / N), abs(λRN)^(1 / N)
    end
  end
  return ψ∞
end


text:
"""
	eulermethod(f, t0, x0; kwargs...)

Solves the differential equation x'(t) = f(t, x(t)) using the euler method.

# Arguments
- `f`: Derivative function
- `t0`: Initial time
- `x0`: Initial position

# Keywords
- `h`: step of the method
- `iter`: number of iterations

# Returns
- `Array{typeof(x0 + h*f(t0, x0))}`: the approximated solution
"""
function eulermethod(f, t0, x0; h = 1e-3, iter = 10000)
	t = t0
	x = Array{typeof(x0 + h*f(t0, x0))}(undef, iter)
	x[1] = x0
	for j = 2:iter
		x[j] = x[j-1] + h*f(t, x[j-1])
		t += h
	end
	return x
end

"""
	targetting_eulermethod(f, dist, t0, x0; kwargs...)

Solves the differential equation x'(t) = f(t, x(t))
using a modification of the euler method to stop close to a given target.
The algorithm is similar to the original euler method,
but the function dist(x) specifies the signed distance from the new point to the target,
if it is less than 0 then it is discarded and the method's time step divided by two,
otherwise it becomes part of the solution.
The algorithm stops when the distance is positive and less than the admitted error
or when the maximum iterations limit is reached.

# Arguments
- `f`: Derivative function
- `dist`: Function for the signed distance to the target
- `t0`: Initial time
- `x0`: Initial position

# Keywords
- `h`: step of the method
- `iter`: number of iterations

# Returns
- `Array{typeof(t0)}, Array{typeof(x0 + h*f(t0, x0))}`: times and values obtained
"""
function targetting_eulermethod(f, dist, t0, x0; h = 1e-3, error = 1e-4, maxit = 100000)
	t = typeof(t0+h)[t0]
	x = Array{typeof(x0 + h*f(t0, x0)), 1}()
	push!(x, x0)
	while dist(last(x)) >= error && maxit > 0
		nx = last(x) + h*f(last(t), last(x))
		d = dist(nx)
		if (d < zero(d))
			h /= 2
		else
			push!(t, last(t)+h)
			push!(x, nx)
		end
		maxit -= 1
	end
	return t, x
end

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
using DataStructures

struct ZerosNearContourException <: Exception
    biggerBox::Tuple
    count:: Int
    potentialLocations::Queue{Tuple{Tuple{Any,Any},Any}}
end

Base.showerror(io::IO, e::ZerosNearContourException) = print(io, "The algorithm has terminated at the box", e.biggerBox, ", containing ", e.count,  " zeros. The remaining list of boxes that have yet to be examined, and the number of zeros they contain is: ", e.potentialLocations, "This is because there are zeros too near the contour during subdivision for the algorithm to work properly at higher precision levels. You can try using different values of arguments (eg. upper_left, lower_right, err, increment) for better results.")

"""
    countZeros(upper_left, lower_right, f::Function, increment=0.01, dif=0.5)::Int

Count zeros of a function in a given rectangular domain

# Example
```julia-repl
julia> countZeros(-5 + 5im, 5 - 5im, x->(x - exp(1) * im)^3 * (x - exp(1))^2)
5
```
"""
function countZeros(a, b, f::Function, increment=0.01, dif=0.5)::Int
    realB = real.(b)
    imagB = imag.(b)

    # Throw exception for invalid input
    if ((imagB >= imag.(a)) || (realB <= real.(a)))
        throw(ErrorException("invalid parameters, please ensure the first parameter is the upper left corner, and the second parameter is the lower right corner"))
    end

    count = 0

    # bottom right to top right
    while imag.(b) < imag.(a)
        angle1 = angle(f(b))
        b += increment * im
        angle2 = angle(f(b))

        if (angle1 - angle2 > 2 * pi - dif)
            count += 1
        elseif (angle1 - angle2 < -2 * pi + dif)
            count -= 1
        end

    end

    # top right to top left
    b = realB + (imag.(a)) * im

    while real.(b) > real.(a)
        angle1 = angle(f(b))
        b -= increment
        angle2 = angle(f(b))

        if (angle1 - angle2 > 2 * pi - dif)
            count += 1
        elseif (angle1 - angle2 < -2 * pi + dif)
            count -= 1
        end

    end

    # top left to bottom left
    while imag.(a) > imagB
        angle1 = angle(f(a))
        a -= increment * im
        angle2 = angle(f(a))

        if (angle1 - angle2 > 2 * pi - dif)
            # println("border, angle1, angle2: ", border, " ", angle1, " ", angle2)
            # println("a, b", a, " ", b)
            count += 1
        elseif (angle1 - angle2 < -2 * pi + dif)
            count -= 1
        end

    end

    # bottom left to bottom right
    a = real.(a) + imagB * im

    while real.(a) < realB
        angle1 = angle(f(a))
        a += increment
        angle2 = angle(f(a))

        if (angle1 - angle2 > 2 * pi - dif)
            count += 1
        elseif (angle1 - angle2 < -2 * pi + dif)
            count -= 1
        end

    end

    return count
end

function getError(a, b)
    (abs(imag.(a) - imag.(b)) + abs(real.(b) - real.(a))) / 2
end

function findZerosWithSubdivisionAux(a, b, f::Function, dif, q::Queue{Tuple{Tuple{Any,Any},Any}}, locations::Array{Tuple{Any,Any}}, err, increment)
    # Check the number of zeros in the initial box. Check whether the sum of sub boxes matches it.

    totalZeros = countZeros(a, b, f, increment, dif)

    if (totalZeros == 0)
        return locations
    end

    realA = real.(a)
    imagA = imag.(a)
    realB = real.(b)
    imagB = imag.(b)
    centerRe = (real.(a) + real.(b)) / 2
    centerIm = (imag.(a) + imag.(b)) / 2
    center = centerRe + centerIm * im

    lowerRightBox = (center, b)
    upperRightBox = (centerRe + imagA * im, realB + centerIm * im)
    upperLeftBox = (a, center)
    lowerLeftBox = (realA + centerIm * im, centerRe + imagB * im)

    count1 = countZeros(lowerRightBox[1], lowerRightBox[2], f, increment, dif)
    count2 = countZeros(upperRightBox[1], upperRightBox[2], f, increment, dif)
    count3 = countZeros(upperLeftBox[1], upperLeftBox[2], f, increment, dif)
    count4 = countZeros(lowerLeftBox[1], lowerLeftBox[2], f, increment, dif)

    if totalZeros != count1 + count2 + count3 + count4
        # this means that a zero is near the contour, stop and just return the existing locations
        throw(ZerosNearContourException((a,b), totalZeros, q))
    end

    # Only enqueue boxes with zeros inside them

    if (count1 != 0)
        enqueue!(q, (lowerRightBox, count1))
    end

    if (count2 != 0)
        enqueue!(q, (upperRightBox, count2))
    end

    if (count3 != 0)
        enqueue!(q, (upperLeftBox, count3))
    end

    if (count4 != 0)
        enqueue!(q, (lowerLeftBox, count4))
    end

    while !isempty(q)
        (x, _) = dequeue!(q)

        if (getError(x[1], x[2]) < err)
            # Error is sufficiently small, this box joins the finalised list of locations
            push!(locations, x)
        elseif (getError(x[1], x[2]) >= err)
            # Further subdivision
            findZerosWithSubdivisionAux(x[1], x[2], f, dif, q, locations, err, increment)
        end

    end

    return locations

end

"""
    findZerosWithSubdivision(upper_left, lower_right, f::Function, err=0.01, increment = err/100, dif=0.5)

Find locations of zeros of a function within a given rectangular domain

Each line of result contains the upper left, lower right corner of the rectangular box that contains the zero

# Arguments
- `upper_left`: a complex number that is the upper left corner of rectangular domain
- `lower_right`: a complex number that is the lower right corner of rectangular domain
- `f::Function`: the function
- `err=0.001`: error of the location of each zero
- `increment=err/100`: step size of the `countZeros` function
- `dif=0.5`: when the difference between ``arg(f(z))`` at 2 consecutive points of
evaluation by the `countZeros` function is greater than `2pi-dif`, a jump point is registered

# Example
```julia-repl
julia> function exp_sum(x)
    return (1+im)*exp(x) + (3.5+2*im) + (3-im)*exp(-x)
end
exp_sum (generic function with 1 method)

julia> findZerosWithSubdivision(-10 + 10im, 10 - 10im, exp_sum, 0.01)
7-element Vector{Tuple{Any, Any}}:
 (0.859375 - 3.0859375im, 0.869140625 - 3.095703125im)
 (0.859375 - 9.375im, 0.869140625 - 9.384765]625im)
 (0.859375 + 9.47265625im, 0.869140625 + 9.462890625im)
 (0.859375 + 3.193359375im, 0.869140625 + 3.18359375im)
 (-0.068359375 + 1.9921875im, -0.05859375 + 1.982421875im)
 (-0.068359375 + 8.271484375im, -0.05859375 + 8.26171875im)
 (-0.068359375 - 4.287109375im, -0.05859375 - 4.296875im)
```
"""
function findZerosWithSubdivision(upper_left, lower_right, f::Function, err=0.001, increment = err/100, dif=0.5)
    q = Queue{Tuple{Tuple{Any,Any},Any}}()
    init = Array{Tuple{Any,Any}}(undef, 0)
    locations = findZerosWithSubdivisionAux(upper_left, lower_right, f::Function, dif, q, init, err, increment)
    return locations
end