text:
using GeometricIntegrators.CommonFunctions
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solutions: createHDF5
using GeometricIntegrators.TestProblems.Oscillator
using HDF5: HDF5File
using Test

ntime = 10
Δt    = .1
ni    = 5

t0    = 0.
x0    = rand(2)
q0    = rand(1)
z0    = [x0[1], x0[2], x0[1] + x0[2]]
p0    = q0.^2
λ0    = [0.]

t1    = 1.
x1    = rand(2, ni)
q1    = rand(1, ni)
p1    = q1.^2
λ1    = zeros(1, ni)

z1    = rand(3, ni)
z1[3,:] .= z0[1,:] .+ z0[2,:]

tx    = zero(x0)
x2    = rand(2, ni)

ode   = oscillator_ode()
dae   = oscillator_dae()
pode  = oscillator_pode()
pdae  = oscillator_pdae()

h5file = "test.hdf5"

@testset "$(rpad("ODE Solution",80))" begin
    sol = Solution(ode, Δt, ntime)
    @test typeof(sol) <: SolutionODE

    sol0 = Solution(similar(ode, x0), Δt, ntime)
    @test typeof(sol0) <: SolutionODE

    sol1 = Solution(similar(ode, x1), Δt, ntime)
    @test typeof(sol1) <: SolutionODE

    @test sol != sol0
    @test sol != sol1

    set_initial_conditions!(sol, t0, x0)
    get_initial_conditions!(sol, tx, 1)

    # @test sol != sol0
    @test tx == x0

    set_initial_conditions!(sol1, similar(ode, t1, x2))
    get_initial_conditions!(sol1, tx, 1)
    @test tx == x2[:,1]

    # test hdf5 in- and output
    h5 = createHDF5(sol, h5file)
    @test typeof(h5) == HDF5File
    close(h5)
    @test isfile(h5file)

    h5 = createHDF5(sol, h5file; overwrite=false)
    @test typeof(h5) == HDF5File
    close(h5)
    @test isfile(h5file)
    rm(h5file)

    write_to_hdf5(sol, h5file)
    @test isfile(h5file)
    rm(h5file)

    create_hdf5(sol, h5file)
    write_to_hdf5(sol)
    close(sol)

    sol2 = SolutionODE(h5file)
    @test sol != sol2
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(ode, Δt, 20, 2)
    @test sol.nt == 10

    sol = Solution(ode, Δt, 20, 2, 10)
    @test sol.nt == 5
end

@testset "$(rpad("PODE Solution",80))" begin
    psol = Solution(pode, Δt, ntime)
    @test typeof(psol) <: SolutionPODE

    psol0 = Solution(similar(pode, q0, p0), Δt, ntime)
    @test typeof(psol0) <: SolutionPODE

    psol1 = Solution(similar(pode, q1, p1), Δt, ntime)
    @test typeof(psol1) <: SolutionPODE

    @test psol != psol0
    @test psol != psol1

    # test hdf5 in- and output
    create_hdf5(psol, "test.hdf5")
    write_to_hdf5(psol)
    close(psol)

    psol2 = SolutionPODE("test.hdf5")
    @test psol != psol2
    rm("test.hdf5")
end

@testset "$(rpad("DAE Solution",80))" begin
    sol = Solution(dae, Δt, ntime)
    @test typeof(sol) <: SolutionDAE

    sol0 = Solution(similar(dae, z0, λ0), Δt, ntime)
    @test typeof(sol0) <: SolutionDAE

    sol1 = Solution(similar(dae, z1, λ1), Δt, ntime)
    @test typeof(sol1) <: SolutionDAE

    @test sol != sol0
    @test sol != sol1

    # test hdf5 in- and output
    create_hdf5(sol, "test.hdf5")
    write_to_hdf5(sol)
    close(sol)

    sol2 = SolutionDAE("test.hdf5")
    @test sol != sol2
    rm("test.hdf5")
end

@testset "$(rpad("PDAE Solution",80))" begin
    psol = Solution(pdae, Δt, ntime)
    @test typeof(psol) <: SolutionPDAE

    psol0 = Solution(similar(pdae, x0, x0.^2, zero(x0)), Δt, ntime)
    @test typeof(psol0) <: SolutionPDAE

    psol1 = Solution(similar(pdae, x1, x1.^2, zero(x1)), Δt, ntime)
    @test typeof(psol1) <: SolutionPDAE

    @test psol != psol0
    @test psol != psol1

    # test hdf5 in- and output
    create_hdf5(psol, "test.hdf5")
    write_to_hdf5(psol)
    close(psol)

    psol2 = SSolutionPDAE("test.hdf5")
    @test psol != psol2
    rm("test.hdf5")
end


text:
function partition!(data::Array, p::Int64, r::Int64)
	i::Int64 = rand(p:r)
	data[i], data[r] = data[r], data[i]
	x = data[r]
	i = p - 1
	for j in p : r - 1
		if data[j] <= x
			i += 1
			data[i], data[j] = data[j], data[i]
		end
	end
	data[i + 1], data[r] = data[r], data[i + 1]
	return i + 1
end

function quicksort!(data::Array, start::Int64, stop::Int64)
	if stop - start < 1
		return
	end

	p = partition!(data, start, stop)

	quicksort!(data, start, p - 1)
	quicksort!(data, p + 1, stop)
end

function quicksort!(data::Array)
	quicksort(data, 1, length(data))
end

ext: jl


text:
#####################
# Credal circuits
#####################
abstract type CredalΔNode{O} <: DecoratorΔNode{O} end
abstract type CredalLeafNode{O} <: CredalΔNode{O} end
abstract type CredalInnerNode{O} <: CredalΔNode{O} end

mutable struct CredalLiteral{O} <: CredalLeafNode{O}
    origin::O
    data
    bit::Bool
    CredalLiteral(n) = new{node_type(n)}(n, nothing, false)
end

mutable struct Credal⋀{O} <: CredalInnerNode{O}
    origin::O
    children::Vector{<:CredalΔNode{<:O}}
    data
    bit::Bool
    Credal⋀(n, children) = begin
        new{node_type(n)}(n, convert(Vector{CredalΔNode{node_type(n)}},children), nothing, false)
    end
end

mutable struct Credal⋁{O} <: CredalInnerNode{O}
    origin::O
    children::Vector{<:CredalΔNode{<:O}}
    log_thetas::Vector{Float64}
    log_thetas_u::Vector{Float64}
    data
    bit::Bool
    Credal⋁(n, children) = new{node_type(n)}(n, convert(Vector{CredalΔNode{node_type(n)}},children), some_vector(Float64, length(children)), some_vector(Float64, length(children)), nothing, false)
end

const CredalΔ{O} = AbstractVector{<:CredalΔNode{<:O}}

Base.eltype(::Type{CredalΔ{O}}) where {O} = CredalΔNode{<:O}

#####################
# traits
#####################

import LogicCircuits.GateType # make available for extension
import LogicCircuits.node_type

@inline GateType(::Type{<:CredalLiteral}) = LiteralGate()
@inline GateType(::Type{<:Credal⋀}) = ⋀Gate()
@inline GateType(::Type{<:Credal⋁}) = ⋁Gate()

@inline node_type(::CredalΔNode) = CredalΔNode

#####################
# constructors and conversions
#####################

const CredalCache = Dict{ΔNode, CredalΔNode}

function CredalΔ2(circuit::Δ)::CredalΔ
    node2dag(CredalΔ2(circuit[end]))
end

function CredalΔ2(circuit::ΔNode)::CredalΔNode
    f_con(n) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")
    f_lit(n) = CredalLiteral(n)
    f_a(n, cn) = Credal⋀(n, cn)
    f_o(n, cn) = Credal⋁(n, cn)
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, CredalΔNode{node_type(circuit)})
end

function CredalΔ(circuit::Δ, cache::CredalCache = CredalCache())

    sizehint!(cache, length(circuit)*4÷3)

    pc_node(::LiteralGate, n::ΔNode) = CredalLiteral(n)
    pc_node(::ConstantGate, n::ΔNode) = error("Cannot construct a probabilistic circuit from constant leafs: first smooth and remove unsatisfiable branches.")

    pc_node(::⋀Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Credal⋀(n, children)
    end

    pc_node(::⋁Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        Credal⋁(n, children)
    end

    map(circuit) do node
        pcn = pc_node(GateType(node), node)
        cache[node] = pcn
        pcn
    end
end

#####################
# methods
#####################

import LogicCircuits: literal, children # make available for extension

@inline literal(n::CredalLiteral)::Lit  = literal(n.origin)
@inline children(n::CredalInnerNode) = n.children

num_parameters(n::Credal⋁) = num_children(n)
num_parameters(c::CredalΔ) = sum(n -> num_parameters(n), ⋁_nodes(c))

"Return the first origin that is a probabilistic circuit node"
prob_origin(n::DecoratorΔNode)::CredalΔNode = origin(n, CredalΔNode)

"Return the first origin that is a probabilistic circuit"
prob_origin(c::DecoratorΔ)::CredalΔ = origin(c, CredalΔNode)

function estimate_credal_parameters2(pc::CredalΔ, data::XData{Bool}, s_idm::Float64; pseudocount::Float64)
    Logical.pass_up_down2(pc, data)
    w = (data isa PlainXData) ? nothing : weights(data)
    estimate_credal_parameters_cached2(pc, w, s_idm; pseudocount=pseudocount)
end

function estimate_credal_parameters_cached2(pc::CredalΔ, w, s_idm::Float64; pseudocount::Float64)
    flow(n) = Float64(sum(sum(n.data)))
    children_flows(n) = sum.(map(c -> c.data[1] .& n.data[1], children(n)))

    if issomething(w)
        flow_w(n) = sum(Float64.(n.data[1]) .* w)
        children_flows_w(n) = sum.(map(c -> Float64.(c.data[1] .& n.data[1]) .* w, children(n)))
        flow = flow_w
        children_flows = children_flows_w
    end

    estimate_credal_parameters_node2(n::CredalΔNode, s_idm::Float64) = ()
    function estimate_credal_parameters_node2(n::Credal⋁, s_idm::Float64)
        if num_children(n) == 1
            n.log_thetas_u .= 0.0
            n.log_thetas .= 0.0
        else
            smoothed_flow = flow(n) + pseudocount
            uniform_pseudocount = pseudocount / num_children(n)
            n.log_thetas .= log.((children_flows(n) .+ uniform_pseudocount)  ./ (smoothed_flow .+ s_idm))
            n.log_thetas_u .= log.(((children_flows(n) .+ uniform_pseudocount) .+ s_idm) ./ (smoothed_flow .+ s_idm))
            # @assert isapprox(sum(exp.(n.log_thetas)), 1.0, atol=1e-6) "Parameters do not sum to one locally"
            # normalize away any leftover error
            # n.log_thetas .- logsumexp(n.log_thetas)
        end
    end

    foreach(estimate_credal_parameters_node2, pc, s_idm)
end

function log_likelihood_per_instance2(pc::CredalΔ, data::XData{Bool})
    Logical.pass_up_down2(pc, data)
    log_likelihood_per_instance_cached(pc, data)
end

function log_likelihood_per_instance_cached(pc::CredalΔ, data::XData{Bool})
    log_likelihoods = zeros(num_examples(data))
    indices = some_vector(Bool, num_examples(data))::BitVector
    for n in pc
         if n isa Credal⋁ && num_children(n) != 1 # other nodes have no effect on likelihood
            foreach(n.children, n.log_thetas) do c, log_theta
                indices = n.data[1] .& c.data[1]
                view(log_likelihoods, indices::BitVector) .+=  log_theta # see MixedProductKernelBenchmark.jl
            end
         end
    end
    log_likelihoods
end

import LogicCircuits: conjoin_like, disjoin_like, literal_like, copy_node, normalize, replace_node # make available for extension

"Conjoin nodes in the same way as the example"
@inline function conjoin_like(example::CredalΔNode, arguments::Vector)
    if isempty(arguments)
        # @assert false "Credalabilistic circuit does not have anonymous true node"
        nothing
    elseif example isa Credal⋀ && children(example) == arguments
        example
    else
        n = conjoin_like(origin(example), origin.(arguments))
        Credal⋀(n, arguments)
    end
end

"Disjoin nodes in the same way as the example"
@inline function disjoin_like(example::CredalΔNode, arguments::Vector)
    if isempty(arguments)
        # @assert false "Probabilistic circuit does not have false node"
        nothing
    elseif example isa Credal⋁ && children(example) == arguments
        example
    else
        n = disjoin_like(origin(example), origin.(arguments))
        # normalize parameters
        thetas = zeros(Float64, length(arguments))
        flag = falses(length(arguments))
        for (i, c) in enumerate(arguments)
            ind = findfirst(x -> x == c, children(example))
            if issomething(ind)
                thetas[i] = exp(example.log_thetas[ind])
                flag[i] = true
            end
        end
        if all(flag)
            thetas = thetas / sum(thetas)
        end
        p = Credal⋁(n, arguments)
        p.log_thetas .= log.(thetas)
        p
    end
end

"Construct a new literal node like the given node's type"
@inline literal_like(::CredalΔNode, lit::Lit) = CredalLiteral(lit)

@inline copy_node(n::Credal⋁, cns) = begin
    orig = copy_node(origin(n), origin.(cns))
    p = Credal⋁(orig, cns)
    p.log_thetas .= copy(n.log_thetas)
    p
end

@inline copy_node(n::Credal⋀, cns) = begin
    orig = copy_node(origin(n), origin.(cns))
    Credal⋀(orig, cns)
end

import LogicCircuits.normalize

@inline normalize(n::Credal⋁, old_n::Credal⋁, kept::Union{Vector{Bool}, BitArray}) = begin
     thetas = exp.(old_n.log_thetas[kept])
     n.log_thetas .= log.(thetas / sum(thetas))
end

function estimate_credal_parameters(pc::CredalΔ, data::XBatches{Bool}, s_idm::Float64; pseudocount::Float64)
    estimate_credal_parameters(AggregateFlowΔ(pc, aggr_weight_type(data)), data, s_idm; pseudocount=pseudocount)
end

function estimate_credal_parameters(afc::AggregateFlowΔ, data::XBatches{Bool}, s_idm::Float64; pseudocount::Float64)
    @assert feature_type(data) == Bool "Can only learn probabilistic circuits on Bool data"
    @assert (afc[end].origin isa CredalΔNode) "AggregateFlowΔ must originate in a CredalΔ"
    collect_aggr_flows(afc, data)
    estimate_credal_parameters_cached(afc, s_idm; pseudocount=pseudocount)
    afc
end

function estimate_credal_parameters(fc::FlowΔ, data::XBatches{Bool}, s_idm::Float64; pseudocount::Float64)
    @assert feature_type(data) == Bool "Can only learn probabilistic circuits on Bool data"
    @assert (prob_origin(afc[end]) isa CredalΔNode) "FlowΔ must originate in a CredalΔ"
    collect_aggr_flows(fc, data)
    estimate_credal_parameters_cached(origin(fc), s_idm; pseudocount=pseudocount)
end

 # turns aggregate statistics into theta parameters
function estimate_credal_parameters_cached(afc::AggregateFlowΔ, s_idm::Float64; pseudocount::Float64)
    foreach(n -> estimate_credal_parameters_node(n, s_idm; pseudocount=pseudocount), afc)
end

estimate_credal_parameters_node(::AggregateFlowΔNode, s_idm::Float64; pseudocount::Float64) = () # do nothing
function estimate_credal_parameters_node(n::AggregateFlow⋁, s_idm::Float64; pseudocount)
    origin = n.origin::Credal⋁
    if num_children(n) == 1
        origin.log_thetas_u .= 0.0
        origin.log_thetas .= 0.0
    else
        smoothed_aggr_flow = (n.aggr_flow + pseudocount)
        uniform_pseudocount = pseudocount / num_children(n)
        origin.log_thetas .= log.( (n.aggr_flow_children .+ uniform_pseudocount) ./ (smoothed_aggr_flow .+ s_idm ))
        origin.log_thetas_u .= log.( ((n.aggr_flow_children .+ uniform_pseudocount) .+ s_idm) ./ (smoothed_aggr_flow .+ s_idm ))
        # @assert isapprox(sum(exp.(origin.log_thetas)), 1.0, atol=1e-6) "Parameters do not sum to one locally: $(exp.(origin.log_thetas)), estimated from $(n.aggr_flow) and $(n.aggr_flow_children). Did you actually compute the aggregate flows?"
        #normalize away any leftover error ()
        # origin.log_thetas .- logsumexp(origin.log_thetas)
    end
end

# compute log likelihood
function compute_log_likelihood(pc::CredalΔ, data::XBatches{Bool})
    compute_log_likelihood(AggregateFlowΔ(pc, aggr_weight_type(data)))
end

# compute log likelihood, reusing AggregateFlowΔ but ignoring its current aggregate values
function compute_log_likelihood(afc::AggregateFlowΔ, data::XBatches{Bool})
    @assert feature_type(data) == Bool "Can only test probabilistic circuits on Bool data"
    collect_aggr_flows(afc, data)
    ll = log_likelihood(afc)
    (afc, ll)
end

# return likelihoods given current aggregate flows.
function log_likelihood(afc::AggregateFlowΔ)
    sum(n -> log_likelihood(n), afc)
end

log_likelihood(::AggregateFlowΔNode) = 0.0
log_likelihood(n::AggregateFlow⋁) = sum(n.origin.log_thetas .* n.aggr_flow_children)

"""
Calculates log likelihood for a batch of fully observed samples.
(Also retures the generated FlowΔ)
"""
function log_likelihood_upper_per_instance(pc::CredalΔ, batch::PlainXData{Bool})
    fc = FlowΔ(pc, num_examples(batch), Bool)
    (fc, log_likelihood_upper_per_instance(fc, batch))
end

function log_likelihood_lower_per_instance(pc::CredalΔ, batch::PlainXData{Bool})
    fc = FlowΔ(pc, num_examples(batch), Bool)
    (fc, log_likelihood_lower_per_instance(fc, batch))
end
"""
Calculate log likelihood per instance for batches of samples.
"""
function log_likelihood_per_instance(pc::CredalΔ, batches::XBatches{Bool})::Vector{Float64}
    mapreduce(b -> log_likelihood_per_instance(pc, b)[2], vcat, batches)
end

function log_prob_upper(pc::CredalΔ, batch::PlainXData{Bool})
    log_likelihood_upper_per_instance(pc, batch)[2]
end

function log_prob_lower(pc::CredalΔ, batch::PlainXData{Bool})
    log_likelihood_lower_per_instance(pc, batch)[2]
end

function log_marginal_upper(pc::CredalΔ, batch::PlainXData{Int8})
    marginal_log_likelihood_upper_per_instance(pc, batch)[2]
end

function log_marginal_lower(pc::CredalΔ, batch::PlainXData{Int8})
    marginal_log_likelihood_lower_per_instance(pc, batch)[2]
end

"""
Calculate log likelihood for a batch of fully observed samples.
(This is for when you already have a FlowΔ)
"""
function log_likelihood_upper_per_instance(fc::FlowΔ, batch::PlainXData{Bool})
    @assert (prob_origin(fc[end]) isa CredalΔNode) "FlowΔ must originate in a CredalΔ"
    pass_up_down(fc, batch)
    log_likelihoods = zeros(num_examples(batch))
    indices = some_vector(Bool, flow_length(fc))::BitVector
    for n in fc
         if n isa DownFlow⋁ && num_children(n) != 1 # other nodes have no effect on likelihood
            origin = prob_origin(n)::Credal⋁
            foreach(n.children, origin.log_thetas) do c, log_theta
                #  be careful here to allow for the Boolean multiplication to be done using & before switching to float arithmetic, or risk losing a lot of runtime!
                # log_likelihoods .+= prod_fast(downflow(n), pr_factors(c)) .* log_theta
                assign_prod(indices, downflow(n), pr_factors(c))
                view(log_likelihoods, indices::BitVector) .+=  log_theta # see MixedProductKernelBenchmark.jl
                # TODO put the lines above in Utils in order to ensure we have specialized types
            end
         end
    end
    log_likelihoods
end

function log_likelihood_lower_per_instance(fc::FlowΔ, batch::PlainXData{Bool})
    @assert (prob_origin(fc[end]) isa CredalΔNode) "FlowΔ must originate in a CredalΔ"
    pass_up_down(fc, batch)
    log_likelihoods = zeros(num_examples(batch))
    indices = some_vector(Bool, flow_length(fc))::BitVector
    for n in fc
         if n isa DownFlow⋁ && num_children(n) != 1 # other nodes have no effect on likelihood
            origin = prob_origin(n)::Credal⋁
            foreach(n.children, origin.log_thetas) do c, log_theta
                #  be careful here to allow for the Boolean multiplication to be done using & before switching to float arithmetic, or risk losing a lot of runtime!
                # log_likelihoods .+= prod_fast(downflow(n), pr_factors(c)) .* log_theta
                assign_prod(indices, downflow(n), pr_factors(c))
                view(log_likelihoods, indices::BitVector) .+=  log_theta # see MixedProductKernelBenchmark.jl
                # TODO put the lines above in Utils in order to ensure we have specialized types
            end
         end
    end
    log_likelihoods
end

"""
Calculate the upper log likelihood for a batch of samples with partial evidence P(e).
(Also returns the generated FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_upper_per_instance(pc::CredalΔ, batch::PlainXData{Int8})
    opts = (flow_opts★..., el_type=Float64, compact⋀=false, compact⋁=false)
    fc = UpFlowΔ(pc, num_examples(batch), Float64, opts)
    (fc, marginal_log_likelihood_upper_per_instance(fc, batch))
end

"""
Calculate the lower log likelihood for a batch of samples with partial evidence P(e).
(Also returns the generated FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_lower_per_instance(pc::CredalΔ, batch::PlainXData{Int8})
    opts = (flow_opts★..., el_type=Float64, compact⋀=false, compact⋁=false)
    fc = UpFlowΔ(pc, num_examples(batch), Float64, opts)
    (fc, marginal_log_likelihood_lower_per_instance(fc, batch))
end

"""
Calculate the upper log likelihood for a batch of samples with partial evidence P(e).
(If you already have a FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_upper_per_instance(fc::UpFlowΔ, batch::PlainXData{Int8})
    @assert (prob_origin(fc[end]) isa CredalΔNode) "FlowΔ must originate in a CredalΔ"
    credal_marginal_upper_pass_up(fc, batch)
    pr(fc[end])
end

"""
Calculate the lower log likelihood for a batch of samples with partial evidence P(e).
(If you already have a FlowΔ)

To indicate a variable is not observed, pass -1 for that variable.
"""
function marginal_log_likelihood_lower_per_instance(fc::UpFlowΔ, batch::PlainXData{Int8})
    @assert (prob_origin(fc[end]) isa CredalΔNode) "FlowΔ must originate in a CredalΔ"
    credal_marginal_lower_pass_up(fc, batch)
    pr(fc[end])
end

"""
Calculate conditional inference upper bound
"""

function conditional_upper(pc::CredalΔ, batch::PlainXData{Int8})
    conditional_upper_per_instance(pc, batch)[2]
end

"""
Calculate conditional inference lower bound
"""

function conditional_lower(pc::CredalΔ, batch::PlainXData{Int8}, mu::Array{Float64,1})
    conditional_lower_per_instance(pc, batch, mu)[2]
end

"""
Calculate conditional inference upper bound per instance
"""
function conditional_upper_per_instance(pc::CredalΔ, batch::PlainXData{Int8})
    opts = (compact⋀=false, compact⋁=false)
    #Initialize the flows as a tuple instead of a single value
    fc = UpFlowΔ(pc, num_examples(batch), Tuple{Float64,Float64,Float64,Float64},opts);
    (fc, conditional_upper_instance(fc, batch))
end

"""
Calculate conditional inference upper bound per instance
"""
function conditional_lower_per_instance(pc::CredalΔ, batch::PlainXData{Int8}, mu::Array{Float64,1})
    opts = (compact⋀=false, compact⋁=false)
    #Initialize the flows as a tuple instead of a single value
    fc = UpFlowΔ(pc, num_examples(batch), Tuple{Float64,Float64,Float64,Float64},opts);
    (fc, conditional_lower_instance(fc, batch, mu))
end

"""
Calculate conditional inference upper bound per instance

"""
function conditional_upper_instance(fc::UpFlowΔ, batch::PlainXData{Int8})
    # @assert (prob_origin(fc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    conditional_upper_pass_up(fc, batch)
    pr(fc[end])
end

"""
Calculate conditional inference lowe bound per instance

"""
function conditional_lower_instance(fc::UpFlowΔ, batch::PlainXData{Int8}, mu::Array{Float64,1})
    # @assert (prob_origin(fc[end]) isa ProbΔNode) "FlowΔ must originate in a ProbΔ"
    conditional_lower_pass_up(fc, batch, mu)
    pr(fc[end])
end
"""
# function check_parameter_integrity(circuit::CredalΔ)
#     for node in filter(n -> GateType(n) isa Credal⋁, circuit)
#         @assert all(θ -> !isnan(θ), node.log_thetas) "There is a NaN in one of the log_thetas"
#     end
#     true
# end

# ##################
# # Sampling from a psdd
# ##################

# """
# Sample from a PSDD without any evidence
# """
# function sample(circuit::CredalΔ)::AbstractVector{Bool}
#     inst = Dict{Var,Int64}()
#     simulate(circuit[end], inst)
#     len = length(keys(inst))
#     ans = Vector{Bool}()
#     for i = 1:len
#         push!(ans, inst[i])
#     end
#     ans
# end

# # Uniformly sample based on the probability of the items
# # and return the selected index
# function sample(probs::AbstractVector{<:Number})::Int32
#     z = sum(probs)
#     q = rand() * z
#     cur = 0.0
#     for i = 1:length(probs)
#         cur += probs[i]
#         if q <= cur
#             return i
#         end
#     end
#     return length(probs)
# end

# function simulate(node::CredalLiteral, inst::Dict{Var,Int64})
#     if positive(node)
#         inst[variable(node.origin)] = 1
#     else
#         inst[variable(node.origin)] = 0
#     end
# end

# function simulate(node::Credal⋁, inst::Dict{Var,Int64})
#     idx = sample(exp.(node.log_thetas))
#     simulate(node.children[idx], inst)
# end
# function simulate(node::Credal⋀, inst::Dict{Var,Int64})
#     for child in node.children
#         simulate(child, inst)
#     end
# end

# """
# Sampling with Evidence from a psdd.
# Internally would call marginal pass up on a newly generated flow circuit.
# """
# function sample(circuit::CredalΔ, evidence::PlainXData{Int8})::AbstractVector{Bool}
#     opts= (compact⋀=false, compact⋁=false)
#     flow_circuit = UpFlowΔ(circuit, 1, Float64, opts)
#     marginal_pass_up(flow_circuit, evidence)
#     sample(flow_circuit)
# end

# """
# Sampling with Evidence from a psdd.
# Assuming already marginal pass up has been done on the flow circuit.
# """
# function sample(circuit::UpFlowΔ)::AbstractVector{Bool}
#     inst = Dict{Var,Int64}()
#     simulate2(circuit[end], inst)
#     len = length(keys(inst))
#     ans = Vector{Bool}()
#     for i = 1:len
#         push!(ans, inst[i])
#     end
#     ans
# end

# function simulate2(node::UpFlowLiteral, inst::Dict{Var,Int64})
#     if positive(node)
#         #TODO I don't think we need these 'grand_origin' parts below
#         inst[variable(grand_origin(node))] = 1
#     else
#         inst[variable(grand_origin(node))] = 0
#     end
# end

# function simulate2(node::UpFlow⋁, inst::Dict{Var,Int64})
#     prs = [ pr(ch)[1] for ch in children(node) ]
#     idx = sample(exp.(node.origin.log_thetas .+ prs))
#     simulate2(children(node)[idx], inst)
# end

# function simulate2(node::UpFlow⋀, inst::Dict{Var,Int64})
#     for child in children(node)
#         simulate2(child, inst)
#     end
# end

# ##################
# # Most Probable Explanation MPE of a psdd
# #   aka MAP
# ##################

# @inline function MAP(circuit::CredalΔ, evidence::PlainXData{Int8})::Matrix{Bool}
#     MPE(circuit, evidence)
# end

# function MPE(circuit::CredalΔ, evidence::PlainXData{Int8})::Matrix{Bool}
#     # Computing Marginal Likelihood for each node
#     fc, lls = marginal_log_likelihood_per_instance(circuit, evidence)

#     ans = Matrix{Bool}(zeros(size(evidence.x)))
#     active_samples = Array{Bool}(ones( num_examples(evidence) ))

#     mpe_simulate(fc[end], active_samples, ans)
#     ans
# end

# """
# active_samples: bool vector indicating which samples are active for this node during mpe
# result: Matrix (num_samples, num_variables) indicating the final result of mpe
# """
# function mpe_simulate(node::UpFlowLiteral, active_samples::Vector{Bool}, result::Matrix{Bool})
#     if positive(node)
#         result[active_samples, variable(node)] .= 1
#     else
#         result[active_samples, variable(node)] .= 0
#     end
# end
# function mpe_simulate(node::UpFlow⋁, active_samples::Vector{Bool}, result::Matrix{Bool})
#     prs = zeros( length(node.children), size(active_samples)[1] )
#     @simd  for i=1:length(node.children)
#         prs[i,:] .= pr(node.children[i]) .+ (node.origin.log_thetas[i])
#     end

#     max_child_ids = [a[1] for a in argmax(prs, dims = 1) ]
#     @simd for i=1:length(node.children)
#         ids = Vector{Bool}( active_samples .* (max_child_ids .== i)[1,:] )  # Only active for this child if it was the max for that sample
#         mpe_simulate(node.children[i], ids, result)
#     end
# end
# function mpe_simulate(node::UpFlow⋀, active_samples::Vector{Bool}, result::Matrix{Bool})
#     for child in node.children
#         mpe_simulate(child, active_samples, result)
#     end
# end


text:
using bamg
using BALUtils

ba = readbal(ARGS[1])
_, jacobian, _ = generate_problem(ba)


text:
using Unitful
using Test
using UnitfulBuckinghamPi
using LinearAlgebra

# Define Unitful paramaters, which can be quantities, units or dimensions
ℓ = u"m"
g = 9.8u"m/s^2"
m = u"g"
T = u"𝐓"
τ = u"s"
θ = u"NoDims"
v = u"m/s"
α = 2
s = "blah"

u = u"m/s"
ρ = u"g/m^3"
μ = u"g/m/s"

@testset "UnitfulBuckinghamPi" begin
    # Set and check parameters
    @setparameters ℓ g m T θ
    @test UnitfulBuckinghamPi.param_symbols == [:ℓ, :g, :m, :T, :θ]
    @test UnitfulBuckinghamPi.param_values == [u"m", 9.8u"m/s^2", u"g", u"𝐓", u"NoDims"]

    # Check adimensional groups as String eltype
    Π = pi_groups(:String)
    @test Π[1] == "g^(1//2)*ℓ^(-1//2)*T^(1//1)"
    @test Π[2] == "θ^(1//1)"

    # Set and check parameters
    @setparameters ℓ g m τ θ
    @test UnitfulBuckinghamPi.param_values ==
        [u"m", 9.8u"m/s^2", u"g", u"s", u"NoDims"]
    # Check adimensional groups as Expr eltype
    Π = pi_groups()
    @test Π[1] == :(g ^ (1 // 2) * ℓ ^ (-1 // 2) * τ ^ (1 // 1))
    @test Π[2] == :(θ ^ (1 // 1))
    # Test evaluating expressions
    @test eval(Π[1]) ≈ 3.1304951684997055
    @test eval(Π[2]) == NoDims

    # Check add parameter and check results
    @addparameters v
    @test UnitfulBuckinghamPi.param_symbols == [:ℓ, :g, :m, :τ, :θ, :v]
    @test UnitfulBuckinghamPi.param_values ==
        [u"m", 9.8u"m/s^2", u"g", u"s", u"NoDims", u"m/s"]
    Π = pi_groups(:Expr)
    @test length(Π) == 3

    # Test parameter of type Number
    @addparameters α
    @test length(UnitfulBuckinghamPi.param_symbols) == 7
    Π = pi_groups()
    @test length(Π) == 4

    # Avoid adding duplicates
    @addparameters α
    @test length(UnitfulBuckinghamPi.param_symbols) == 7

    # Check setting no parameters for an empty list
    @setparameters
    Π = pi_groups()
    @test size(Π) == (0,)

    # Check singularity in the LU decomposition
    @setparameters u ρ μ ℓ
    @test pi_groups() ==
        [:(ρ ^ (1 // 1) * u ^ (1 // 1) * μ ^ (-1 // 1) * ℓ ^ (1 // 1))]
    @setparameters ℓ ρ μ u
    @test pi_groups() ==
        [:(ρ ^ (1 // 1) * μ ^ (-1 // 1) * ℓ ^ (1 // 1) * u ^ (1 // 1))]

    # Test errors
    @setparameters ℓ g m τ θ
    @test_throws ArgumentError pi_groups(:NotImplemented)
    @test_throws MethodError @setparameters 1
    @test_throws ArgumentError @setparameters s
end

@testset "LU with full pivoting" begin
    # Example from LinearAlgebra.lu
    A = [4 3; 6 3]
    F = UnitfulBuckinghamPi.lu_pq(A)
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L ≈ [1.0 0.0; 0.666667 1.0] (atol = 0.00001)
    @test (F.L == L) && (F.U == U) && (F.p == p) && (F.q == q)
    @test L * U == A[p, q]

    # Singular rectangular wide matrix for which regular `LinearAlgebra.lu` fails
    A = reshape(collect(1:12),3,4)
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]
    @test rank(A) == rank(U) == 2

    # Same but with Rational{Int} eltype
    A = convert.(Rational, reshape(collect(1:12),3,4))
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]
    @test (L isa Matrix{Rational{Int64}}) && (U isa Matrix{Rational{Int64}})

    # Tall matrix version
    A = convert.(Rational, reshape(collect(1:12),4,3))
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]
    @test U == [
        12//1   4//1   8//1;
        0//1  -2//1  -1//1;
        0//1   0//1   0//1;
        0//1   0//1   0//1;
    ]

    # Complex{Int} version
    A = complex(reshape(collect(1:12),3,4))
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]
    @test (L isa Matrix{Complex{Float64}}) && (U isa Matrix{Complex{Float64}})

    # Complex{Float} version
    A = complex(reshape(collect(1.0:12.0),3,4))
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]
    @test (L isa Matrix{Complex{Float64}}) && (U isa Matrix{Complex{Float64}})

    # SVD values  very far apart, for which `LinearAlgebra.rank` fails
    A = [1 1; typemax(Int64)//2 1]
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U == A[p, q]

    # Singular case associated with the second Reynolds number example above
    A = convert.(Rational,[-1 0 -2 -1; 0 1 1 1; 1 -3 -1 -1])
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test U == [
        -3//1 -1//1 -1//1 1//1;
        0//1  -2//1 -1//1 -1//1;
        0//1  0//1  1//3  0//1
    ]

    # Complex{Int} version, with round-offs
    A = convert.(Complex,[-1 0 -2 -1; 0 1 1 1; 1 -3 -1 -1])
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test (L isa Matrix{Complex{Float64}}) && (U isa Matrix{Complex{Float64}})
    @test U ≈ [
        -3.0+0.0im -1.0+0.0im  -1.0+0.0im 1.0+0.0im;
        0.0+0.0im  -2.0+0.0im -1.0+0.0im -1.0+0.0im;
        0.0+0.0im   0.0+0.0im  0.333333+0.0im  0.0+0.0im] (atol = 0.00001)

    # Random wide matrix
    A = rand(7,11)
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U ≈ A[p,q]

    # Random tall matrix
    A = rand(11,7)
    L, U, p, q = UnitfulBuckinghamPi.lu_pq(A)
    @test L * U ≈ A[p,q]
end

max_stars_repo_name: rmsrosa/UnitfulBuckinghamPi

max_issues_repo_name: rmsrosa/UnitfulBuckinghamPi

max_forks_repo_name: rmsrosa/UnitfulBuckinghamPi


text:
## L-Shaped Benders Decomposition for Stochastic Programming
## Edward J. Xu <edxu96@outlook.com>
## August 11, 2019

"MILP model for L-Shaped Benders decomposition (LSBD) in matrix form."
mutable struct ModelLSBD
  n_x
  vec_min_y
  vec_max_y
  vec_c
  vec_f
  vec_b
  mat_aCap
  mat_bCap
  solution::Union{Solution, Missing}

  function ModelLSBD(
      vec_min_y::Array{Int64,1}, vec_max_y::Array{Int64,1},
      vec_c::Array{Int64,2}, vec_b::Array{Int64,2},
      vec_f::Array{Int64,2}, mat_aCap::Array{Int64,2},
      mat_bCap::Array{Int64,2}
      )
    checkListLength(vec_min_y, vec_max_y, 'vec_min_y', 'vec_max_y')
    check_col_vec(vec_c, "vec_c")
    check_col_vec(vec_f, "vec_f")
    check_col_vec(vec_b, "vec_b")
    check_mat_match(vec_c, mat_aCap, row, col, "vec_c", "mat_aCap")
    check_mat_match(vec_f, mat_bCap, row, col, "vec_f", "mat_bCap")
    check_mat_match(mat_aCap, mat_bCap, row, row, "mat_aCap", "mat_bCap")
    check_mat_match(mat_aCap, vec_b, row, row, "mat_aCap", "vec_b")
    new(vec_c, mat_aCap, vec_f, mat_bCap, vec_b, missing)
  end
end

"First two variables are updated."
function solve_master!(obj_mas, model_mas, e1_mat, e2, opt_cut::Bool)
  if opt_cut
    @constraint(model_mas, (e1_mat * vec_y)[1] + q >= e2)
  else
    ## Add feasible cut Constraints
    @constraint(model_mas, (e1_mat * vec_y)[1] >= e2)
    @constraint(model_mas, (e1_mat * vec_y)[1] + q >= e2)
  end

  optimize!(model_mas)
  obj_mas = objective_value(model_mas)
end

"First two variables are updated."
function solve_sub!(vec_ubar, obj_sub, vec_ybar, n_constraint, vec_h , mat_t,
    mat_w, vec_c)

  model_sub = Model(with_optimizer(GLPK.Optimizer))
  @variable(model_sub, vec_u[1: n_constraint] >= 0)
  @objective(model_sub, Max, (transpose(vec_h - mat_t * vec_ybar) * vec_u)[1])
  cons_dual = @constraint(model_sub, transpose(mat_w) * vec_u .<= vec_c)
  sol = optimize!(model_sub)

  print("  Sub Problem")
  vec_ubar = value.(vec_u)
  if sol == :OPTIMAL
    bool_sub = true
    obj_sub = objective_value(model_sub)
    vec_result_x = dual(cons_dual)
  elseif sol == :DUAL_INFEASIBLE # ???
    print("Not solved optimally because the feasible set is unbounded.\n")
    bool_sub = false
    obj_sub = objective_value(model_sub)
    vec_result_x = repeat([NaN], length(vec_c))
  elseif sol == :INFEASIBLE
    print("Not solved optimally because of infeasibility. Something is ",
      "wrong.\n")
    bool_sub = false
    obj_sub = NaN
    vec_ubar = hcat(vec_ubar)
    vec_result_x = hcat(repeat([NaN], length(vec_c)))
  end

  return bool_sub, vec_result_x
end

"""
Solve the sub-problem with ray.
  First two variables are updated.
"""
function solve_ray!(vec_ubar, obj_ray, vec_ybar, n_constraint, vec_h, mat_t,
    mat_w)

  model_ray = Model(with_optimizer(GLPK.Optimizer))
  @variable(model_ray, vec_u[1: n_constraint] >= 0)
  @objective(model_ray, Max, 1)
  @constraint(model_ray, (transpose(vec_h - mat_t * vec_ybar) * vec_u)[1] == 1)
  @constraint(model_ray, transpose(mat_w) * vec_u .<= 0)
  optimize!(model_ray)

  print("  Ray Problem\n")
  vec_ubar = value.(vec_u)
  obj_ray = objective_value(model_ray)
end

"""
L-Shaped Benders Decomposition for Stochastic Programming without Integer
  Variables in Second Stage.
"""
function lshaped(; n_x, vec_min_y, vec_max_y, vec_f, vec_pi, mat_c, mat_h,
    mat3_t, mat3_w, epsilon=1e-6, timesIterationMax=100)

  println("Begin L-Shaped Benders Decomposition")
  n_y = length(vec_min_y)
  num_s = length(mat3_t[:, 1, 1])
  n_constraint = length(mat3_w[1, :, 1])
  mod_mas = set_mod_mas(n_y, vec_min_y, vec_max_y)

  let
    ub = Inf  # upper bound
    lb = - Inf  # lower bound

    # initial value of master variables
    mat_uBar = zeros(num_s, n_constraint, 1)
    vec_ybar = zeros(n_y, 1)
    vec_result_x = zeros(num_s, Int8(n_x / num_s), 1)
    obj_sub = 0
    obj_sub_s = zeros(num_s)
    timesIteration = 1

    dict_obj_mas = Dict()
    dict_q = Dict()
    dict_obj_sub = Dict()
    dict_obj_ray = Dict()
    dict_ub = Dict()
    dict_lb = Dict()

    # Must make sure "result_q == obj_sub" in the final iteration
    # while ((ub - lb > epsilon) && (timesIteration <= timesIterationMax))  !!!
    while check_whe_continue(ub, lb, epsilon, result_q, obj_sub,
        timesIteration, timesIterationMax)
      ## 1. Solve sub/ray problem for each scenario
      bool_sub_s = trues(num_s)
      for s = 1: num_s
        bool_sub_s[s], vec_result_x[s, :, :] = solve_sub!(mat_uBar[s, :, :],
          obj_sub_s[s], vec_ybar, n_constraint, mat_h[s, :, :], mat3_t[s, :, :],
          mat3_w[s, :, :], mat_c[s, :, :])
        if !(bool_sub_s[s])
          solve_ray!(mat_uBar[s, :, :], obj_sub_s[s], vec_ybar, n_constraint,
            mat_h[s, :, :], mat3_t[s, :, :], mat3_w[s, :, :])
        end
      end
      obj_sub = (transpose(vec_pi) * obj_sub_s)[1]
      ub = min(ub, obj_sub + (transpose(vec_f) * vec_ybar)[1])

      ## 2. Add optimal cut to master problem
      e1_mat = sum(vec_pi[s] * (transpose(mat_uBar[s, :]) * mat3_t[s, :, :])[1]
        for s = 1: num_s)
      e2 = sum(vec_pi[s] * (transpose(mat_uBar[s, :]) * mat_h[s, :])[1] for
        s = 1: num_s)

      solve_master!(obj_mas, e1_mat, e2, bool_sub_s[1])
      vec_ybar = value.(vec_y)

      ## 3. Compare the bounds and decide whether to stop
      lb = max(lb, obj_mas)
      result_q = value.(q)

      dict_ub[timesIteration] = ub
      dict_lb[timesIteration] = lb
      dict_obj_mas[timesIteration] = obj_mas
      dict_obj_sub[timesIteration] = obj_sub
      dict_q[timesIteration] = result_q

      if bool_sub_s[1]
        println("------------------ Result in $(timesIteration)-th Iteration with Sub ",
            "-------------------\n", "ub: $(round(ub, digits = 5)), ",
            "lb: $(round(lb, digits = 5)), obj_mas: $(round(obj_mas, digits = 5)), ",
            "q: $result_q, obj_sub: $(round(obj_sub, digits = 5)).")
      else
        println("------------------ Result in $(timesIteration)-th Iteration with Ray ",
            "-------------------\n", "ub: $(round(ub, digits = 5)), ",
            "lb: $(round(lb, digits = 5)), obj_mas: $(round(obj_mas, digits = 5)), ",
            "q: $result_q, obj_ray: $(round(obj_ray, digits = 5)).")
      end
      timesIteration += 1
    end

    println("obj_mas: $(objective_value(model_mas))")
    println("-------------------------------------------------------------------------\n",
        "------------------------------ 2/4. Result ------------------------------\n",
        "-------------------------------------------------------------------------")
    println("ub: $(round(ub, digits = 5)), lb: $(round(lb, digits = 5)), ",
        "difference: $(round(ub - lb, digits = 5))")
    println("vec_x: $vec_result_x")
    vec_result_y = value.(vec_y)
    result_q = value.(q)
    println("vec_y: $vec_result_y")
    println("result_q: $result_q")

    println("-------------------------------------------------------------------------\n",
        "------------------------- 3/4. Iteration Result -------------------------\n",
        "-------------------------------------------------------------------------")
    # Initialize
    seq_timesIteration = collect(1: (timesIteration - 1))
    vec_ub = zeros(timesIteration - 1)
    vec_lb = zeros(timesIteration - 1)
    obj_sub_sRay = zeros(timesIteration - 1)
    vec_obj_mas = zeros(timesIteration - 1)
    vec_q = zeros(timesIteration - 1)
    vec_type = repeat(["ray"], (timesIteration - 1))

    #
    for i = 1: (timesIteration - 1)
      vec_obj_mas[i] = round(dict_obj_mas[i], digits = 5)
      vec_ub[i] = round(dict_ub[i], digits = 5)
      vec_lb[i] = round(dict_lb[i], digits = 5)
      vec_q[i] = round(dict_q[i], digits = 5)
      if haskey(dict_obj_sub, i)
        vec_type[i] = "sub"
        obj_sub_sRay[i] = round(dict_obj_sub[i], digits = 5)
      else
        obj_sub_sRay[i] = round(dict_obj_ray[i], digits = 5)
      end
    end
    table_iterationResult = hcat(seq_timesIteration, vec_ub, vec_lb,
      vec_obj_mas, vec_q, vec_type, obj_sub_sRay)
    pretty_table(table_iterationResult,
      ["Seq", "ub", "lb", "obj_mas", "q", "sub/ray", "obj_sub/ray"],
      compact; alignment=:l)
  end

  println("-------------------------------------------------------------------------\n",
      "-------------------------- 4/4. Nominal Ending --------------------------\n",
      "-------------------------------------------------------------------------\n")
end


text:
# while statement
i = 0

while i<5
    println("BOOM!!")
    global i+=1
end

# for statement
for i = 1:5
    println("BOOM!! $i")
end

# for-in statement
arr = ["foo", "bar", "BOOM!!!", "HO"]

for e in arr
    println(e)
end

# nested for (9x9 multiply)
for i = 1:9, j = 1:9
    println("$(i)X$(j)=$(i*j)")
end

# nested for-in
for b ∈ [true, false], e ∈ [1,2,3,4]
    println("$b:$e")
end

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
# Utilities for stocks-related endpoints

import HTTP: request, escapeuri
import JSON: parse

valid_range_values = [
    "max", "5y", "2y", "1y", "ytd", "6m", "3m", "1m", "1mm", "5d", "5dm",
    "date", "dynamic", "next"
]

"""
    get_stocks(conn::Connection, endpoint::String, path_params::Array{Symbol, 1};
        kwargs...)

Get stock prices in various formats
"""
function get_stocks(
    conn::Connection,
    endpoint::String,
    path_params::Array{Symbol, 1};
    kwargs...
)
    kwargs = Dict(kwargs)

    path_dict = Dict()
    query_dict = Dict()
    for (k, v) in kwargs
        if k in path_params
            merge!(path_dict, Dict(k => v))
        else
            merge!(query_dict, Dict(k => v))
        end
    end

    symbol = try
        escapeuri(path_dict[:symbol])
    catch error
        if isa(error, KeyError)
            println("Need a symbol for this endpoint")
        end
    end

    range = get(path_dict, :range, "1m")
    if !(range in valid_range_values)
        throw(ArgumentError("invalid range value"))
    end
    date = get(path_dict, :date, "")

    if isempty(date)
        http_string = "stock/$symbol/$endpoint/$range"
    else
        http_string =  "stock/$symbol/$endpoint/date/$date"
    end

    url = join([conn.url, conn.version, http_string], "/")
    if haskey(query_dict, :symbols)
        query_dict[:symbols] = escapeuri(query_dict[:symbols])
    end
    params = merge(Dict("token" => conn.token), query_dict)
    r = try
        request("GET", url; query = params)
    catch e
        throw(IEXException(e.status))
    end

    return parse(String(r.body))
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
export iauDtdb
"""
An approximation to TDB-TT, the difference between barycentric
dynamical time and terrestrial time, for an observer on the Earth.

The different time scales - proper, coordinate and realized - are
related to each other:

         TAI             <-  physically realized
            :
         offset            <-  observed (nominally +32.184s)
            :
         TT              <-  terrestrial time
            :
   rate adjustment (L_G)   <-  definition of TT
            :
         TCG             <-  time scale for GCRS
            :
   "periodic" terms      <-  iauDtdb  is an implementation
            :
   rate adjustment (L_C)   <-  function of solar-system ephemeris
            :
         TCB             <-  time scale for BCRS
            :
   rate adjustment (-L_B)  <-  definition of TDB
            :
         TDB             <-  TCB scaled to track TT
            :
   "periodic" terms      <-  -iauDtdb is an approximation
            :
         TT              <-  terrestrial time

Adopted values for the various constants can be found in the IERS
Conventions (McCarthy & Petit 2003).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2   double  date, TDB (Notes 1-3)
   ut            double  universal time (UT1, fraction of one day)
   elong         double  longitude (east positive, radians)
   u             double  distance from Earth spin axis (km)
   v             double  distance north of equatorial plane (km)

Returned (function value):
                  double  TDB-TT (seconds)

Notes:

   1. The date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   Although the date is, formally, barycentric dynamical time (TDB),
   the terrestrial dynamical time (TT) can be used with no practical
   effect on the accuracy of the prediction.

   2. TT can be regarded as a coordinate time that is realized as an
   offset of 32.184s from International Atomic Time, TAI.  TT is a
   specific linear transformation of geocentric coordinate time TCG,
   which is the time scale for the Geocentric Celestial Reference
   System, GCRS.

   3. TDB is a coordinate time, and is a specific linear transformation
   of barycentric coordinate time TCB, which is the time scale for
   the Barycentric Celestial Reference System, BCRS.

   4. The difference TCG-TCB depends on the masses and positions of the
   bodies of the solar system and the velocity of the Earth.  It is
   dominated by a rate difference, the residual being of a periodic
   character.  The latter, which is modeled by the present function,
   comprises a main (annual) sinusoidal term of amplitude
   approximately 0.00166 seconds, plus planetary terms up to about
   20 microseconds, and lunar and diurnal terms up to 2 microseconds.
   These effects come from the changing transverse Doppler effect
   and gravitational red-shift as the observer (on the Earth's
   surface) experiences variations in speed (with respect to the
   BCRS) and gravitational potential.

   5. TDB can be regarded as the same as TCB but with a rate adjustment
   to keep it close to TT, which is convenient for many applications.
   The history of successive attempts to define TDB is set out in
   Resolution 3 adopted by the IAU General Assembly in 2006, which
   defines a fixed TDB(TCB) transformation that is consistent with
   contemporary solar-system ephemerides.  Future ephemerides will
   imply slightly changed transformations between TCG and TCB, which
   could introduce a linear drift between TDB and TT;  however, any
   such drift is unlikely to exceed 1 nanosecond per century.

   6. The geocentric TDB-TT model used in the present function is that of
   Fairhead & Bretagnon (1990), in its full form.  It was originally
   supplied by Fairhead (private communications with P.T.Wallace,
   1990) as a Fortran subroutine.  The present C function contains an
   adaptation of the Fairhead code.  The numerical results are
   essentially unaffected by the changes, the differences with
   respect to the Fairhead & Bretagnon original being at the 1e-20 s
   level.

   The topocentric part of the model is from Moyer (1981) and
   Murray (1983), with fundamental arguments adapted from
   Simon et al. 1994.  It is an approximation to the expression
   ( v / c ) . ( r / c ), where v is the barycentric velocity of
   the Earth, r is the geocentric position of the observer and
   c is the speed of light.

   By supplying zeroes for u and v, the topocentric part of the
   model can be nullified, and the function will return the Fairhead
   & Bretagnon result alone.

   7. During the interval 1950-2050, the absolute accuracy is better
   than +/- 3 nanoseconds relative to time ephemerides obtained by
   direct numerical integrations based on the JPL DE405 solar system
   ephemeris.

   8. It must be stressed that the present function is merely a model,
   and that numerical integration of solar-system ephemerides is the
   definitive method for predicting the relationship between TCG and
   TCB and hence between TT and TDB.

References:

   Fairhead, L., & Bretagnon, P., Astron.Astrophys., 229, 240-247
   (1990).

   IAU 2006 Resolution 3.

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Moyer, T.D., Cel.Mech., 23, 33 (1981).

   Murray, C.A., Vectorial Astrometry, Adam Hilger (1983).

   Seidelmann, P.K. et al., Explanatory Supplement to the
   Astronomical Almanac, Chapter 2, University Science Books (1992).

   Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G. & Laskar, J., Astron.Astrophys., 282, 663-683 (1994).

This revision:  2018 January 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauDtdb(double date1, double date2,
#                double ut, double elong, double u, double v)

function iauDtdb(date1::Real, date2::Real,
                 ut::Real, elong::Real,
                 u::Real, v::Real)

   tdb_tt = ccall((:iauDtdb, libsofa_c), Cdouble,
            (Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble),
            convert(Float64, date1), convert(Float64, date2),
            convert(Float64, ut), convert(Float64, elong),
            convert(Float64, u), convert(Float64, v))

   return tdb_tt
end


text:
#
# PivotLDU is a block LDU decomposition of
# [I B;
#  C I]
# where
#
# size(B) = (r1,r2)
# size(C) = (r2,r1)
#

import Base.LinAlg: Factorization, LU

immutable PivotLDU{T,S<:AbstractMatrix} <: Factorization{T}
    B::S
    C::S
    factor::LU{T,S} # LU factorization of I-C*B
    r1::Int
    r2::Int
end

function pivotldufact(A::AbstractMatrix,r1,r2)
    @assert size(A) == (r1+r2,r1+r2)
    B = A[1:r1,r1+1:end]
    C = A[r1+1:end,1:r1]
    PivotLDU(B,C,lufact(I-C*B),r1,r2)
end

function A_ldiv_B1B2!{T<:Number,S}(P::PivotLDU{T,S},b1::AbstractArray{T},b2::AbstractArray{T})
    b2[:] = b2 - P.C*b1
    A_ldiv_B!(P.factor,b2)
    b1[:] = b1 - P.B*b2
    b1,b2
end

function Base.A_ldiv_B!{T<:Number,S}(P::PivotLDU{T,S},b::AbstractVector{T})
    b[1+P.r1:end] -= P.C*b[1:P.r1]
    b[1+P.r1:end] = P.factor\b[1+P.r1:end]
    b[1:P.r1] -= P.B*b[1+P.r1:end]
    b
end

function Base.A_ldiv_B!{T<:Number,S}(P::PivotLDU{T,S},b::AbstractMatrix{T})
    b[1+P.r1:end,:] .-= P.C*b[1:P.r1,:]
    b[1+P.r1:end,:] = P.factor\b[1+P.r1:end,:]
    b[1:P.r1,:] .-= P.B*b[1+P.r1:end,:]
    b
end


text:
using StatsBase

"""
Returns all vowels as an array of a word input as a string
"""
function vowel_checker(word::String)
    vowels = []
    for character in word
        if character in "aeiou"
            push!(vowels,character)
        end
    end
    return vowels
end

"""
Returns an array of double letters, empty if no double letters
"""
function double_letter_checker(word::String)
    double_letter_list = []
    previous_character = word[1]
    for character in word[2:end]
        if character == previous_character
            push!(double_letter_list,string(previous_character,character))
        else
            previous_character = character
        end
    end
    return double_letter_list
end

"""
Returns true if the word contains a naughty substring, false if not
"""
function naughty_string_checker(word::String)
    naughty_strings = ["ab","cd","pq","xy"]
    for naughty_string in naughty_strings
        if occursin(naughty_string,word) == true
            return true
        end
    end
    return false
end

# Open Input
cd(raw"/Users/wapisani/Documents/Programming/Julia/Advent_of_Code/")
open("AoC_2015_5.input") do f
    global words_list = split(read(f,String),'\n')
end

nice_list = []
naughty_list = []
# words_list = ["dvszwmarrgswjxmb"]

for word in words_list
    # Convert substring into string
    word = String(word)
    condition_counter = 0
    vowels = vowel_checker(word)
    if length(vowels) >= 3
        condition_counter += 1
    end

    double_letters = double_letter_checker(word)
    if length(double_letters) >= 1
        condition_counter += 1
    end

    if naughty_string_checker(word) == true
        push!(naughty_list,word)
        continue
    else
        condition_counter += 1
    end

    if condition_counter == 3
        push!(nice_list,word)
    else
        push!(naughty_list,word)
    end

end

num_nice_words = length(nice_list)
println("According to Santa's first set of rules, there are $num_nice_words nice words in Santa's list")

"""
Returns true if the word input as a String contains at least one letter which repeats with exactly one letter between them, false if not
"""
function every_other_letter(word::String)
    for i=2:length(word)-1
        if word[i-1] == word[i+1]
            return true
        end
    end
    return false
end

nice_list_2 = []
# words_list = ["aaa"]
for word in words_list
    word = String(word)
    condition_counter = 0

    double_regex = r"([a-z][a-z]).*\1"
    if occursin(double_regex,word) == true
        condition_counter += 1
    end

    every_other_letter_condition = every_other_letter(word)
    if every_other_letter_condition == true
        condition_counter += 1
    end

    if condition_counter == 2
        push!(nice_list_2,word)
    end
end

num_nice_words_2 = length(nice_list_2)
println("According to Santa's second set of rules, there are $num_nice_words_2 nice words in Santa's list")


text:
using Test

include("BamReader/runtests.jl")


text:
@testset "test frame constructors" begin
    f = M.frame(;name="test")
    @test isa(f, M.PlotlyFrame)
    f = M.PlotlyFrame{Dict{Symbol,Any}}(Dict{Symbol,Any}(:name=>"test"))
    @test isa(f, M.PlotlyFrame)
    @test_logs (:warn, "Frame should have a :name field for expected behavior") frame()
    @test_logs (:warn, "Frame should have a :name field for expected behavior") PlotlyFrame{Dict{Symbol,Any}}(Dict{Symbol,Any}())
end
@testset "_UNDERSCORE_ATTRS" begin
    f = M.frame(;name="test")
    # test setindex!
    f[:paper_bgcolor] = "grey"
    @test haskey(f.fields, :paper_bgcolor)
    @test f.fields[:paper_bgcolor] == "grey"

    # test getindex
    @test f[:paper_bgcolor] == "grey"
    @test f["paper_bgcolor"] == "grey"

    # now do it again with the string form of setindex!
    f = M.Layout()

    # test setindex!
    f["paper_bgcolor"] = "grey"
    @test haskey(f.fields, :paper_bgcolor)
    @test f.fields[:paper_bgcolor] == "grey"

    # test getindex
    @test f[:paper_bgcolor] == "grey"
    @test f["paper_bgcolor"] == "grey"
end

@testset "test setting nested attr" begin
    f = M.frame(;name="test")
    times20 = attr(name="Times", size=20)
    f[:xaxis_titlefont] = times20
    @test isa(f[:xaxis], Dict)
    @test f[:xaxis][:titlefont][:name] == "Times"
    @test f[:xaxis][:titlefont][:size] == 20
end

@testset "test setindex!, getindex methods" begin
    f = M.frame(;name="test", x=1:5, y=1:5, visible=false)

    f[:visible] = true
    @test length(f.fields) == 4
    @test haskey(f.fields, :visible)
    @test f.fields[:visible] == true

    # now try with string. Make sure it updates inplace
    f["visible"] = false
    @test length(f.fields) == 4
    @test haskey(f.fields, :visible)
    @test f.fields[:visible] == false

    # -------- #
    # 2 levels #
    # -------- #
    f[:line, :color] = "red"
    @test length(f.fields) == 5
    @test haskey(f.fields, :line)
    @test isa(f.fields[:line], Dict)
    @test f.fields[:line][:color] == "red"
    @test f["line.color"] == "red"

    # now try string version
    f["line", "color"] = "blue"
    @test length(f.fields) == 5
    @test haskey(f.fields, :line)
    @test isa(f.fields[:line], Dict)
    @test f.fields[:line][:color] == "blue"
    @test f["line_color"] == "blue"

    # now try convenience string dot notation
    f["line.color"] = "green"
    @test length(f.fields) == 5
    @test haskey(f.fields, :line)
    @test isa(f.fields[:line], Dict)
    @test f.fields[:line][:color] == "green"
    @test f[:line_color] == "green"

    # now try symbol with underscore
    f[:(line_color)] = "orange"
    @test length(f.fields) == 5
    @test haskey(f.fields, :line)
    @test isa(f.fields[:line], Dict)
    @test f.fields[:line][:color] == "orange"
    @test f["line.color"] == "orange"

    # now try string with underscore
    f["line_color"] = "magenta"
    @test length(f.fields) == 5
    @test haskey(f.fields, :line)
    @test isa(f.fields[:line], Dict)
    @test f.fields[:line][:color] == "magenta"
    @test f["line.color"] == "magenta"

    # -------- #
    # 3 levels #
    # -------- #
    f[:marker, :line, :color] = "red"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test haskey(f.fields[:marker], :line)
    @test isa(f.fields[:marker][:line], Dict)
    @test haskey(f.fields[:marker][:line], :color)
    @test f.fields[:marker][:line][:color] == "red"
    @test f["marker.line.color"] == "red"

    # now try string version
    f["marker", "line", "color"] = "blue"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test haskey(f.fields[:marker], :line)
    @test isa(f.fields[:marker][:line], Dict)
    @test haskey(f.fields[:marker][:line], :color)
    @test f.fields[:marker][:line][:color] == "blue"
    @test f["marker.line.color"] == "blue"

    # now try convenience string dot notation
    f["marker.line.color"] = "green"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test haskey(f.fields[:marker], :line)
    @test isa(f.fields[:marker][:line], Dict)
    @test haskey(f.fields[:marker][:line], :color)
    @test f.fields[:marker][:line][:color] == "green"
    @test f["marker.line.color"] == "green"

    # now string with underscore notation
    f["marker_line_color"] = "orange"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test haskey(f.fields[:marker], :line)
    @test isa(f.fields[:marker][:line], Dict)
    @test haskey(f.fields[:marker][:line], :color)
    @test f.fields[:marker][:line][:color] == "orange"
    @test f["marker.line.color"] == "orange"

    # now symbol with underscore notation
    f[:(marker_line_color)] = "magenta"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test haskey(f.fields[:marker], :line)
    @test isa(f.fields[:marker][:line], Dict)
    @test haskey(f.fields[:marker][:line], :color)
    @test f.fields[:marker][:line][:color] == "magenta"
    @test f["marker.line.color"] == "magenta"

    # -------- #
    # 4 levels #
    # -------- #
    f[:marker, :colorbar, :tickfont, :family] = "Hasklig-ExtraLight"
    @test length(f.fields) == 6  # notice we didn't add another top level key
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test length(f.fields[:marker]) == 2  # but we did add a key at this level
    @test haskey(f.fields[:marker], :colorbar)
    @test isa(f.fields[:marker][:colorbar], Dict)
    @test haskey(f.fields[:marker][:colorbar], :tickfont)
    @test isa(f.fields[:marker][:colorbar][:tickfont], Dict)
    @test haskey(f.fields[:marker][:colorbar][:tickfont], :family)
    @test f.fields[:marker][:colorbar][:tickfont][:family] == "Hasklig-ExtraLight"
    @test f["marker.colorbar.tickfont.family"] == "Hasklig-ExtraLight"

    # now try string version
    f["marker", "colorbar", "tickfont", "family"] = "Hasklig-Light"
    @test length(f.fields) == 6
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test length(f.fields[:marker]) == 2
    @test haskey(f.fields[:marker], :colorbar)
    @test isa(f.fields[:marker][:colorbar], Dict)
    @test haskey(f.fields[:marker][:colorbar], :tickfont)
    @test isa(f.fields[:marker][:colorbar][:tickfont], Dict)
    @test haskey(f.fields[:marker][:colorbar][:tickfont], :family)
    @test f.fields[:marker][:colorbar][:tickfont][:family] == "Hasklig-Light"
    @test f["marker.colorbar.tickfont.family"] == "Hasklig-Light"

    # now try convenience string dot notation
    f["marker.colorbar.tickfont.family"] = "Hasklig-Medium"
    @test length(f.fields) == 6  # notice we didn't add another top level key
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test length(f.fields[:marker]) == 2  # but we did add a key at this level
    @test haskey(f.fields[:marker], :colorbar)
    @test isa(f.fields[:marker][:colorbar], Dict)
    @test haskey(f.fields[:marker][:colorbar], :tickfont)
    @test isa(f.fields[:marker][:colorbar][:tickfont], Dict)
    @test haskey(f.fields[:marker][:colorbar][:tickfont], :family)
    @test f.fields[:marker][:colorbar][:tickfont][:family] == "Hasklig-Medium"
    @test f["marker.colorbar.tickfont.family"] == "Hasklig-Medium"

    # now string with underscore notation
    f["marker_colorbar_tickfont_family"] = "Webdings"
    @test length(f.fields) == 6  # notice we didn't add another top level key
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test length(f.fields[:marker]) == 2  # but we did add a key at this level
    @test haskey(f.fields[:marker], :colorbar)
    @test isa(f.fields[:marker][:colorbar], Dict)
    @test haskey(f.fields[:marker][:colorbar], :tickfont)
    @test isa(f.fields[:marker][:colorbar][:tickfont], Dict)
    @test haskey(f.fields[:marker][:colorbar][:tickfont], :family)
    @test f.fields[:marker][:colorbar][:tickfont][:family] == "Webdings"
    @test f["marker.colorbar.tickfont.family"] == "Webdings"

    # now symbol with underscore notation
    f[:marker_colorbar_tickfont_family] = "Webdings42"
    @test length(f.fields) == 6  # notice we didn't add another top level key
    @test haskey(f.fields, :marker)
    @test isa(f.fields[:marker], Dict)
    @test length(f.fields[:marker]) == 2  # but we did add a key at this level
    @test haskey(f.fields[:marker], :colorbar)
    @test isa(f.fields[:marker][:colorbar], Dict)
    @test haskey(f.fields[:marker][:colorbar], :tickfont)
    @test isa(f.fields[:marker][:colorbar][:tickfont], Dict)
    @test haskey(f.fields[:marker][:colorbar][:tickfont], :family)
    @test f.fields[:marker][:colorbar][:tickfont][:family] == "Webdings42"
    @test f["marker.colorbar.tickfont.family"] == "Webdings42"

    # error on 5 levels
    @test_throws MethodError f["marker.colorbar.tickfont.family.foo"] = :bar
end

@testset "create plot with frames" begin
    # https://plotly.com/python/animations/
    p = Plot(
        [scatter(x=0:1, y=0:1)],
        Layout(xaxis=attr(range=(0,5), autorange=false), yaxis=attr(range=(0,5), autorange=false), title="start title", updatemenus=[attr(type="buttons", buttons=[attr(label="Play", method="animate", args=[nothing])])]),
        [
            frame(name=:f1, data=[scatter(x=1:2, y=1:2)]),
            frame(name=:f2, data=[scatter(x=1:4, y=1:4)]),
            frame(name=:f3, data=[scatter(x=3:4, y=3:4)], layout=Layout(title_text="End Title")),
        ]
    )
    @test p isa Plot
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using TOML
import TOML: linecol, whitespace, comment, newline, expect, lookup, Parser, parse

using Test
using Dates

macro testval(s, v)
    f = "foo = $s"
    :( @test TOML.get(parse(Parser($f)))["foo"] == $v )
end

macro fail(s...)
    teststr = s[1]
    errstr = ""
    debug = false
    if length(s) == 2
        if isa(s[2], String)
            errstr = s[2]
        elseif isa(s[2], Expr)
            debug = s[2].args[1] == :debug
        end
    end
    if length(s) == 3
        errstr = s[2]
        debug = s[3].args[1] == :debug
    end

    # macro vars
    pvar = :pv
    ppvar = :ppv

    # debugging report
    dbgexp = if debug
        quote
            println("\nTEST FAIL: ", escape_string($teststr))
            if length($pvar.errors) > 0
                println("ERRORS:")
                for e in $pvar.errors
                    println(e)
                end
            end
            if !TOML.isnull($ppvar)
                println("RESULT:")
                println(TOML.get($ppvar))
            end
        end
    else
        :()
    end

    # error for comparison
    errtsts = if !isempty(errstr)
        quote
            @test length($pvar.errors)>0
            @test findlast(m->m==$errstr, map(e->e.msg, $pvar.errors)) != 0
        end
    else
        :()
    end

    quote
        local $pvar = Parser($teststr)
        local $ppvar = parse($pvar)
        $dbgexp
        $errtsts
        @test TOML.isnull($ppvar)
    end
end

macro success(s...)
    teststr = s[1]
    debug = false
    if length(s) == 2 && isa(s[2], Expr)
        debug = s[2].args[1] == :debug
    end

    # macro vars
    pvar = :pv
    ppvar = :ppv

    # debugging report
    dbgexp = if debug
        quote
            println("\nTEST SUCCESS: ", escape_string($teststr))
            if length($pvar.errors) > 0
                println("ERRORS:")
                for e in $pvar.errors
                    println(e)
                end
            end
            if !TOML.isnull($ppvar)
                println("RESULT:")
                println(TOML.get($ppvar))
            end
        end
    else
        :()
    end

    quote
        local $pvar = Parser($teststr)
        local $ppvar = parse($pvar)
        $dbgexp
        @test !TOML.isnull($ppvar)
    end
end

@testset "TOML" begin
@testset "Parser" begin

    @testset "Parser internal functions" begin
        test = """
   #[test]\r
   [test]
   foo = "bar"
"""
        p = Parser(test)
        @test whitespace(p)
        @test position(p) == 4
        @test comment(p)
        @test position(p) == 13
        @test whitespace(p)
        @test position(p) == 16
        @test !expect(p, 'a')
        @test expect(p, '[')
        @test position(p) == 17
    end

    @testset "Lookups" begin
        p = Parser("""hello."world\\t".a.0.'escaped'.value""")
        @testset for (p, s) in zip(TOML.get(lookup(p)), ["hello"; "world\t"; "a"; "0"; "escaped"; "value"])
            @test p == s
        end

        p = Parser("")
        @test TOML.get(lookup(p)) == String[]

        p = Parser("value")
        @test TOML.get(lookup(p)) == String["value"]

        p = Parser("\"\"")
        #TODO: @test get(lookup(p)) == String[""]

    end

    @testset "New line" begin
        @success("""
[project.A]\r
name = \"splay2\"\r
\r
[project]\r
\r
name = \"splay\"\r
version = \"0.1.0\"\r
authors = [\"alex@crichton.co\"]\r
\r
[[lib]]\r
\r
path = \"lib.rs\"\r
name = \"splay\"\r
description = \"\"\"
A Rust implementation of a TAR file reader and writer. This library does not\r
currently handle compression, but it is abstract over all I/O readers and\r
writers. Additionally, great lengths are taken to ensure that the entire\r
contents are never required to be entirely resident in memory all at once.\r
\"\"\"
[[lib.aaa]]\r
bbb = \"aaa\"\r
""")

        @fail("\r")

        @fail("a = [ \r ]")

        @fail("a = \"\"\"\r\"\"\"")

        @fail("a = \"\"\"\\  \r  \"\"\"")

        @success("a = '''\r'''")

        @success("a = '\r'")

        @fail("a = \"\n\"")

        @fail("a = '\n'")

        p = Parser("
            foo = \"\"\"\\\r\n\"\"\"
            bar = \"\"\"\\\r\n   \r\n   \r\n   a\"\"\"
        ")
        res = TOML.get(parse(p))
        @test res["foo"] == ""
        @test res["bar"] == "a"

        @fail("0=0r=false", "expected a newline after a key")

        @fail("""
0=""o=""m=""r=""00="0"q=\"\"\"0\"\"\"e=\"\"\"0\"\"\"
""", "expected a newline after a key")

        @fail("""
[[0000l0]]
0="0"[[0000l0]]
0="0"[[0000l0]]
0="0"l="0"
""", "expected a newline after a key")

        @fail("""
0=[0]00=[0,0,0]t=["0","0","0"]s=[1000-00-00T00:00:00Z,2000-00-00T00:00:00Z]
""", "expected a newline after a key")

        @fail("""
0=0r0=0r=false
""", "expected a newline after a key")

        @fail("""
0=0r0=0r=falsefal=false
""", "expected a newline after a key")

    end

    @testset "Offset->(Line, Col) " begin
        p = Parser("ab\ncde\nf")
        @test linecol(p, 0) == (1,0)
        @test linecol(p, 1) == (1,1)
        @test linecol(p, 4) == (2,1)
        @test linecol(p, 5) == (2,2)
        @test linecol(p, 8) == (3,1)
        @test linecol(p, 9) == (3,0)
    end

    @testset "Strings" begin
        p = Parser("""
bar = "\\U00000000"
key1 = "One\\nTwo"
key2 = \"\"\"One\nTwo\"\"\"
key3 = \"\"\"
One
Two\"\"\"

key4 = "The quick brown fox jumps over the lazy dog."
key5 = \"\"\"
The quick brown \\

  fox jumps over \\
    the lazy dog.\"\"\"
key6 = \"\"\"\\
       The quick brown \\
       fox jumps over \\
       the lazy dog.\\
       \"\"\"

# What you see is what you get.
winpath  = 'C:\\Users\\nodejs\\templates'
winpath2 = '\\\\ServerX\\admin\$\\system32\\'
quoted   = 'Tom "Dubs" Preston-Werner'
regex    = '<\\i\\c*\\s*>'
regex2 = '''I [dw]on't need \\d{2} apples'''
lines  = '''
The first newline is
trimmed in raw strings.
    All other whitespace
    is preserved.
'''
""")
        res = TOML.get(parse(p))

        @test res["bar"]  == "\0"
        @test res["key1"] == "One\nTwo"
        @test res["key2"] == "One\nTwo"
        @test res["key3"] == "One\nTwo"

        msg = "The quick brown fox jumps over the lazy dog."
        @test res["key4"] == msg
        @test res["key5"] == msg
        @test res["key6"] == msg

        @test res["winpath"] == "C:\\Users\\nodejs\\templates"
        @test res["winpath2"] == "\\\\ServerX\\admin\$\\system32\\"
        @test res["quoted"] == """Tom "Dubs" Preston-Werner"""
        @test res["regex"] == """<\\i\\c*\\s*>"""
        @test res["regex2"] == """I [dw]on't need \\d{2} apples"""

        @test res["lines"] == "The first newline is\ntrimmed in raw strings.\n    All other whitespace\n    is preserved.\n"

        @testval("''", "")
        @testval("\"\"", "")
        @testval("\"\"\"\n\n\n\"\"\"", "\n\n")

        #@fail("foo = \"\\uD800\"")

        @fail("foo = \"\\uxx\"")
        @fail("foo = \"\\u\"")
        @fail("foo = \"\\")
        @fail("foo = '")

    end

    @testset "Numbers" begin
        @fail("a = 00")
        @fail("a = -00")
        @fail("a = +00")
        @fail("a = 00.0")
        @fail("a = -00.0")
        @fail("a = +00.0")
        @fail("a = 9223372036854775808")
        @fail("a = -9223372036854775809")

        @fail("a = 0.")
        @fail("a = 0.e")
        @fail("a = 0.E")
        @fail("a = 0.0E")
        @fail("a = 0.0e")
        @fail("a = 0.0e-")
        @fail("a = 0.0e+")
        @fail("a = 0.0e+00")

        @testval("1.0", 1.0)
        @testval("1.0e0", 1.0)
        @testval("1.0e+0", 1.0)
        @testval("1.0e-0", 1.0)
        @testval("1.001e-0", 1.001)
        @testval("2e10", 2e10)
        @testval("2e+10", 2e10)
        @testval("2e-10", 2e-10)
        @testval("2_0.0", 20.0)
        @testval("2_0.0_0e0_0", 20.0)
        @testval("2_0.1_0e1_0", 20.1e10)

        @fail("4")

        @testval("1_0", 10)
        @testval("1_0_0", 100)
        @testval("1_000", 1000)
        @testval("+1_000", 1000)
        @testval("-1_000", -1000)

        @fail("foo = 0_")
        @fail("foo = 0__0")
        @fail("foo = __0")
        @fail("foo = 1_0_")

    end

    @testset "Booleans" begin

        @testval("true", true)
        @testval("false", false)

        @fail("foo = true2")
        @fail("foo = false2")
        @fail("foo = t1")
        @fail("foo = f1")

    end

    @testset "Datetime" begin

        @testval("2016-09-09T09:09:09Z", DateTime(2016,9,9,9,9,9))
        @testval("2016-09-09T09:09:09.0Z", DateTime(2016,9,9,9,9,9))
        @testval("2016-09-09T09:09:09.0+10:00", DateTime(2016,9,9,19,9,9))
        @testval("2016-09-09T09:09:09.012-02:00", DateTime(2016,9,9,7,9,9,12))

        @fail("foo = 2016-09-09T09:09:09.Z", "malformed date literal")
        @fail("foo = 2016-9-09T09:09:09Z", "malformed date literal")
        @fail("foo = 2016-09-09T09:09:09+2:00", "malformed date literal")
        @fail("foo = 2016-09-09T09:09:09-2:00", "malformed date literal")
        @fail("foo = 2016-09-09T09:09:09Z-2:00", "expected a newline after a key")

    end

    @testset "Keys" begin

        p = Parser("
            foo = 3
            foo_3 = 3
            foo_-2--3--r23f--4-f2-4 = 3
            _ = 3
            - = 3
            8 = 8
            \"a\" = 3
            \"!\" = 3
            \"a^b\" = 3
            \"\\\\\\\"\" = 3
            \"character encoding\" = \"value\"
            'ʎǝʞ' = \"value\"
        ")
        res = TOML.get(parse(p))

        @test haskey(res, "foo")
        @test haskey(res, "-")
        @test haskey(res, "_")
        @test haskey(res, "8")
        @test haskey(res, "foo_3")
        @test haskey(res, "foo_-2--3--r23f--4-f2-4")
        @test haskey(res, "a")
        @test haskey(res, "!")
        @test haskey(res, "\\\"")
        @test haskey(res, "character encoding")
        @test haskey(res, "ʎǝʞ")

        @fail("key\n=3")
        @fail("key=\n3")
        @fail("key|=3")
        @fail("\"\"=3")
        @fail("=3")
        @fail("\"\"|=3")
        @fail("\"\n\"|=3")
        @fail("\"\r\"|=3")

    end

    @testset "Tables" begin

        @fail("[]")
        @fail("[.]")
        @fail("[\"\".\"\"]")
        @fail("[a.]")
        @fail("[\"\"]")
        @fail("[!]")
        @fail("[\"\n\"]")
        @fail("[a.b]\n[a.\"b\"]")
        @fail("[']")
        @fail("[''']")
        @fail("['''''']")
        @fail("['\n']")
        @fail("['\r\n']")

        p = Parser("
            [a.\"b\"]
            [\"f f\"]
            [\"f.f\"]
            [\"\\\\\\\"\"]
            ['a.a']
            ['\"\"']
        ")
        res = TOML.get(parse(p))
        @test haskey(res, "a.a")
        @test haskey(res, "f f")
        @test haskey(res, "f.f")
        @test haskey(res, "\\\"")
        @test haskey(res, "\"\"")
        @test haskey(res["a"], "b")

        @test haskey(TOML.get(parse(Parser("[foo]"))), "foo")

        @testset "Inline Tables" begin

            @success("a = {}")
            @success("a = {b=1}")
            @success("a = {   b   =   1    }")
            @success("a = {a=1,b=2}")
            @success("a = {a=1,b=2,c={}}")
            @fail("a = {a=1,}")
            @fail("a = {,}")
            @fail("a = {a=1,a=1}")
            @fail("a = {\n}")
            @fail("a = {")
            @success("a = {a=[\n]}")
            @success("a = {\"a\"=[\n]}")
            @success("a = [\n{},\n{},\n]")

        end

        @testset "Redefinition" begin
            @fail("
                [a]
                foo=\"bar\"
                [a.b]
                foo=\"bar\"
                [a]
            ", "redefinition of table `a`")
            @fail("
                [a]
                foo=\"bar\"
                b = { foo = \"bar\" }
                [a]
            ", "redefinition of table `a`")
            @success("
                [a.b]
                c = {}
                [a]
                d = 1
            ")
            @fail("
                [a.b]
                c = 1
                [a]
                b = {}
            ", "duplicate key `b`")
            @fail("
                [a]
                b = {}
                [a]
            ", "redefinition of table `a`")
        end

        @testset "Nesting" begin
            @fail("
                a = [2]
                [[a]]
                b = 5
            ", "expected type `TOML.Table`, found type `$(Int)`")
            @fail("
                a = 1
                [a.b]
            ", "key `a` was not previously a table")
            @fail("
                a = []
                [a.b]
            ", "array `a` does not contain tables")
            @success("
                a = [{}]
                [a.b]
            ")
            @fail("
                a = []
                [[a.b]]
            ", "array `a` does not contain tables")
            @fail("
                [a]
                b = { c = 2, d = {} }
                [a.b]
                c = 2
            ", "duplicate key `c`")
        end
    end

    @testset "Arrays" begin

        p = Parser("""
[[foo]]
  #…
  [foo.bar]
    #…
[[foo]] # ...
  #…
  [foo.bar]
    #...
""")
        res = TOML.get(parse(p))
        @test haskey(res, "foo")
        arr = res["foo"]
        @test length(arr) == 2
        @test isa(arr[1], TOML.Table)
        @test haskey(arr[1], "bar")
        @test haskey(arr[2], "bar")

        p = Parser("""
[[fruit]]
  name = "apple"
  [fruit.physical]
    color = "red"
    shape = "round"
  [[fruit.variety]]
    name = "red delicious"
  [[fruit.variety]]
    name = "granny smith"
[[fruit]]
  name = "banana"
  [[fruit.variety]]
    name = "plantain"
""")
        res = TOML.get(parse(p))
        @test haskey(res, "fruit")
        fruit = res["fruit"]
        @test length(fruit) == 2
        @test isa(fruit[1], TOML.Table)
        apple = fruit[1]
        banana = fruit[2]

        @test haskey(apple, "name")
        @test apple["name"] == "apple"
        @test haskey(apple, "physical")
        @test apple["physical"]["color"] == "red"
        @test apple["physical"]["shape"] == "round"
        @test apple["variety"][1]["name"] == "red delicious"
        @test apple["variety"][2]["name"] == "granny smith"

        @test haskey(banana, "name")
        @test banana["name"] == "banana"
        @test banana["variety"][1]["name"] == "plantain"

        res = TOML.parse("""
[[products]]
name = "Hammer"
sku = 738594937

[[products]]

[[products]]
name = "Nail"
sku = 284758393
color = "gray"
""")
        @test haskey(res, "products")
        products = res["products"]
        @test isa(products, Array)
        @test length(products) == 3
        @test products[1]["name"] == "Hammer"
        @test products[1]["sku"] == 738594937
        @test length(products[2]) == 0
        @test products[3]["name"] == "Nail"
        @test products[3]["sku"] == 284758393
        @test products[3]["color"] == "gray"

        @fail("""array = [
        #         "This might most likely happen in multiline arrays",
        #         Like here,
        #         "or here,
        #         and here"
        #         ]     End of array comment, forgot the #
        """)
    end

end

@testset "Printer" begin
    res1 = TOML.parse("""
title = "TOML Example"
[owner]
bio = \"\"\"GitHub Cofounder & CEO
Likes tater tots and beer.\"\"\"
dob = 1979-05-27T07:32:00Z # First class dates? Why not?

[database]
server = "192.168.1.1"
ports = [ 8001, 8001, 8002 ]
enabled = true

[servers]

  # You can indent as you please. Tabs or spaces. TOML don't care.
  [servers.alpha]
  ip = "10.0.0.1"
  dc = "eqdc10"

  [servers.beta]
  ip = "10.0.0.2"
  dc = "eqdc10"
  country = "中国" # This should be parsed as UTF-8

[clients]
data = [ ["gamma", "delta"], [1, 2] ] # just an update to make sure parsers support it

  [[products]]
  name = "Hammer"
  sku = 738594937

  [[products]]
  name = "Nail"
  sku = 284758393
  color = "gray"
""")
    io = IOBuffer()
    TOML.print(io, res1)
    seekstart(io)
    res2 = TOML.parse(io)
    @test res1 == res2

end
end


text:
#=
shepplogan.jl
=#

export shepp_logan, ellipse_parameters
export SheppLogan
export SheppLoganEmis
export SheppLoganBrainWeb
export SheppLoganToft
export SouthPark

"""
    EllipsePhantomVersion
Parent type for different versions of ellipse phantoms:
* `SheppLogan` original CT version from Shepp&Logan paper
* `SheppLoganEmis` higher contrast version suitable for emission tomography
* `SheppLoganBrainWeb` integer index version based on brainweb
* `SheppLoganToft` higher contrast version from Toft, 1996
* `SouthPark` for fun
"""
abstract type EllipsePhantomVersion end

"""
    SheppLogan
Original version from :
Larry A Shepp, Benjamin F Logan,
"The Fourier reconstruction of a head section,"
IEEE Transactions on Nuclear Science, 21(3):21-42, June 1974.
[doi](http://doi.org/10.1109/TNS.1974.6499235)

Also in Kak and Slaney 1988 text, p. 255.
[doi](http://doi.org/10.1137/1.9780898719277)
"""
struct SheppLogan <: EllipsePhantomVersion end

"""
    SheppLoganToft
Toft, Peter Aundal & Sørensen, John Aasted
"The Radon transform-theory and implementation,"
Technical University of Denmark (DTU), 1996. Page 201.
https://files.openpdfs.org/1ra51GP6gJO.pdf
"""
struct SheppLoganToft <: EllipsePhantomVersion end

struct SheppLoganEmis <: EllipsePhantomVersion end
struct SheppLoganBrainWeb <: EllipsePhantomVersion end
struct SouthPark <: EllipsePhantomVersion end

"""
    values = shepp_logan_values(::EllipsePhantomVersion)
Return 10 Shepp-Logan ellipse amplitudes for various versions.
"""
shepp_logan_values(::SheppLogan) =
    [2.0, -0.98, -0.02, -0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
shepp_logan_values(::SheppLoganEmis) =
    [1, 1, -2, 2, 3, 4, 5, 6, 1, 1]
shepp_logan_values(::SheppLoganBrainWeb) =
    [1, 0, 2, 3, 4, 5, 6, 7, 8, 9] # brainweb uses index 1-10
shepp_logan_values(::SheppLoganToft) =
    [1.0, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

"""
    phantom = Ellipse(n × 6 AbstractMatrix)
Return vector of Ellipse objects, one for each row of input matrix.
"""
function Ellipse(params::AbstractMatrix{<:RealU})
    size(params,2) == 6 || throw("ellipses need 6 parameters")
    return [Ellipse(params[n,:]) for n in 1:size(params,1)]
end

function shepp_logan(case::EllipsePhantomVersion; kwargs...)
    return Ellipse(ellipse_parameters(case; kwargs...))
end

"""
    image = shepp_logan(M, [N,], case, options...)
Convenience method for generating `M×N` samples of Shepp-Logan phantoms.

In
* `M::Int` : horizontal size
* `N::Int` : vertical size, defaults to `M`
* `case::EllipsePhantomVersion = SheppLogan()`

# Options
* `oversample::Int = 3` (usually)
* `yflip::Bool = true` (reverse `y` samples for convenience.)
* `kwargs...` remaining options passed to `ellipse_parameters` for parameters.

# Out
* `image` : `M × N` matrix

The default here is 3× over-sampling along both axes (9 samples per pixel),
except for the `SheppLoganBrainWeb` phantom that consists of integer indices.
"""
function shepp_logan(
    M::Int, N::Int,
    case::EllipsePhantomVersion = SheppLogan() ;
    oversample::Int = case == SheppLoganBrainWeb() ? 1 : 3,
    yflip::Bool = true,
    kwargs...
)
    ob = shepp_logan(case ; kwargs...)
    x = LinRange(-0.5, 0.5, M)
    y = LinRange(-0.5, 0.5, N)
    if yflip
        y = reverse(y)
    end
    return oversample > 1 ?
        phantom(x, y, ob, oversample) :
        phantom(x, y, ob) # type unstable because of over-sampling
end

shepp_logan(M::Int, case::EllipsePhantomVersion = SheppLogan(); kwargs...) =
    shepp_logan(M, M, case ; kwargs...)

"""
    params = ellipse_parameters(case::Symbol; fovs::NTuple{2}, u::Tuple, disjoint::Bool)

By default the first four columns are unitless "fractions of field of view",
so columns 1,3 are scaled by `xfov` and columns 2,4 are scaled by `yfov`,
where `(xfov, yfov) = fovs`.
The optional 3-tuple `u` specifies scaling and/or units:
* columns 1-4 (center, radii) are scaled by `u[1]` (e.g., mm),
* column 5 (angle) is scaled by `u[2]` (e.g., `1` or `°`),
* column 6 (value) is scaled by `u[3]` (e.g., `1/cm`) for an attenuation map.
If `disjoint==true` then the middle ellipse positions are adjusted to avoid overlap.
"""
function ellipse_parameters(
    case::EllipsePhantomVersion = SheppLogan() ;
    fovs::NTuple{2,RealU} = (1,1),
    u::NTuple{3,Any} = (1,1,1), # unit scaling
    disjoint::Bool = false,
)
    (xfov, yfov) = fovs

    # The "Number" here is to enable units
    params = Number[ # original CT version
    0       0       0.92    0.69    90    2     # skull
    0       -0.0184 0.874   0.6624  90    -0.98 # brain
    0.22    0       0.31    0.11    72    -0.02 # right big
    -0.22   0       0.41    0.16   108    -0.02 # left big
    0       0.35    0.25    0.21    90    0.01 # top
    0       0.1     0.046   0.046    0    0.01 # middle high
    0       -0.1    0.046   0.046    0    0.01 # middle low
    -0.08   -0.605  0.046   0.023    0    0.01 # bottom left
    0       -0.605  0.023   0.023    0    0.01 # bottom center
    0.06    -0.605  0.046   0.023   90    0.01 # bottom right
    ]

    if disjoint # tweak so that middle ellipses are disjoint
        params[4,1] += -0.04
        params[5,2] += 0.07
    end

    params[:,[1,3]] .*= xfov/2
    params[:,[2,4]] .*= yfov/2
    params[:,5] .*= π/180
    params[:,6] = shepp_logan_values(case)

    return Number[params[:,1:4] * u[1] params[:,5] * u[2] params[:,6] * u[3]]
end

"""
    params = ellipse_parameters(SouthPark() ; fovs::NTuple{2,Number} = (100,100))
Ellipse parameters for "South Park" phantom.
"""
function ellipse_parameters(::SouthPark ; fovs::NTuple{2,RealU} = (100,100))
    params = Float64[
        0 0 85 115 0 100
        0 -60 30 20 0 -80 # mouth
        30 20 25 35 30 20 # eyes
        -30 20 25 35 -30 20
        35 25 7 7 0 -100 # pupils
        -15 25 7 7 0 -100
        0 75 60 15 0 -50  # hat
    ]
    params[:,[1,3]] .*= fovs[1] / 256
    params[:,[2,4]] .*= fovs[2] / 256
    params[:,5] .*= π/180
    return params
end


text:
# We just define some convenience function to make the experience feel a bit more polished.
# These are probably not too useful in practice though.

const UNARY_FUNCS = (:-, :~, :conj, :abs,
                  :sin, :cos, :tan, :sinh, :cosh, :tanh,
                  :asin, :acos, :atan, :asinh, :acosh, :atanh,
                  :sec, :csc, :cot, :asec, :acsc, :acot,
                  :sech, :csch, :coth, :asech, :acsch, :acoth,
                  :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
                  :sind, :tand, :acosd, :acotd, :acscd, :asecd,
                  :asind, :atand, :rad2deg, :deg2rad,
                  :log, :log2, :log10, :log1p, :exponent, :exp,
                  :exp2, :expm1, :cbrt, :sqrt, :erf,
                  :erfc, :erfcx, :erfi, :dawson, :ceil, :floor,
                  :trunc, :round, :significand)

const BINARY_FUNCS  = (:.+, :.-,:.*, :./, :.\, :.^,
                    :+, :-,
                   :min, :max,
                   :div, :mod)

const BOOLEAN_BINARY_FUNCS = (:.==, :.!=, :.<, :.<=, :.>, :.>=)

const REDUCTION_FUNCS = ((:maximum, :max),
                         (:minimum,  :min),
                         (:sum, :.+),
                         (:prod, :.*))

####################
# Binary functions #
####################

macro blockwise_binary(func, Tres)
    esc(quote
        @assert A.block_sizes == B.block_sizes
        block_array_new = similar(A, $Tres)
        @inbounds for I in eachindex(A.blocks)
            block_array_new.blocks[I] = $func(A.blocks[I], B.blocks[I])
        end
        return block_array_new
    end)
end

macro blockwise_binary_pseudo(func, Tres)
    esc(quote
        @assert A.block_sizes == B.block_sizes
        return PseudoBlockArray($func(A.blocks, B.blocks), copy(A.block_sizes))
    end)
end

macro blockwise_binary_scalar(func, Tres, mode)
    esc(quote
        block_array_new = similar(block_array, $Tres)
        @inbounds for I in eachindex(block_array.blocks)
            if $mode == 1
                block_array_new.blocks[I] = $func(A, block_array.blocks[I])
            else
                block_array_new.blocks[I] = $func(block_array.blocks[I], A)
            end
        end
        return block_array_new
    end)
end

macro blockwise_binary_scalar_pseudo(func, Tres, mode)
    esc(quote
        if $mode == 1
            return PseudoBlockArray($func(A, block_array.blocks), copy(block_array.block_sizes))
        else
            return PseudoBlockArray($func(block_array.blocks, A), copy(block_array.block_sizes))
        end
    end)
end

for f in BINARY_FUNCS
    @eval begin
        Base.$f{T1, T2, N}(A::BlockArray{T1, N}, B::BlockArray{T2, N}) = @blockwise_binary($f, Base.promote_op($f, T1, T2))
        Base.$f{T1, T2, N}(A::PseudoBlockArray{T1, N}, B::PseudoBlockArray{T2, N}) = @blockwise_binary_pseudo($f, Base.promote_op($f, T1, T2))

        # These generate ambiquous errors when called
        #Base.$f{T1 <: Number, T2, N}(A::T1, block_array::BlockArray{T2, N}) = @blockwise_binary_scalar($f, Base.promote_op($f, T1, T2), 1)
        #Base.$f{T1, T2 <: Number, N}(block_array::BlockArray{T1, N}, A::T2) = @blockwise_binary_scalar($f, Base.promote_op($f, T1, T2), 2)
        #Base.$f{T1 <: Number, T2, N}(A::T1, block_array::PseudoBlockArray{T2, N}) = @blockwise_binary_scalar_pseudo($f, Base.promote_op($f, T1, T2), 1)
        #Base.$f{T1, T2 <: Number, N}(block_array::PseudoBlockArray{T1, N}, A::T2) = @blockwise_binary_scalar_pseudo($f, Base.promote_op($f, T1, T2), 2)
    end
end

for f in BOOLEAN_BINARY_FUNCS
    @eval begin
        Base.$f{T1, T2, N}(A::BlockArray{T1, N}, B::BlockArray{T2, N}) = @blockwise_binary($f, Bool)
        Base.$f{T1, T2, N}(A::PseudoBlockArray{T1, N}, B::PseudoBlockArray{T2, N}) = @blockwise_binary_pseudo($f, Bool)
    end
end

###################
# Unary functions #
###################

macro blockwise_unary(func, Tres)
    esc(quote
        block_array_new = similar(block_array, $Tres)
        @inbounds for I in eachindex(block_array.blocks)
            block_array_new.blocks[I] = $func(block_array.blocks[I])
        end
        return block_array_new
    end)
end

macro blockwise_unary_pseudo(func, Tres)
    esc(quote
        return PseudoBlockArray($func(block_array.blocks), copy(block_array.block_sizes))
    end)
end

for f in UNARY_FUNCS
    @eval begin
        Base.$f{T}(block_array::BlockArray{T}) = @blockwise_unary($f, typeof($f(one(T))))
        Base.$f{T}(block_array::PseudoBlockArray{T}) = @blockwise_unary_pseudo($f, typeof($f(one(T))))
    end
end

##############
# Reductions #
##############

macro blockwise_reduction(func, merging_func)
    esc(quote
        blocks = block_array.blocks
        @assert length(blocks) > 0
        s = $func(blocks[1])
        @inbounds for i in 2:length(blocks)
            s = $merging_func(s, $func(blocks[i]))
        end
        return s
    end)
end

macro blockwise_reduction_pseudo(func, merging_func)
    esc(quote
        return $func(block_array.blocks)
    end)
end

for (f, f_merge) in REDUCTION_FUNCS
    @eval begin
        Base.$f(block_array::BlockArray) = @blockwise_reduction($f, $f_merge)
        Base.$f(block_array::PseudoBlockArray) = @blockwise_reduction_pseudo($f, $f_merge)
    end
end


text:
include("snoop_bot_config.jl")

snoop_bot(
    botconfig,
    joinpath(@__DIR__, "precompile_script.jl"),
)


text:
using DashBase
using JSON2

@testset "components creation" begin
    test_comp = Component("html_div", "Div", "dash_html_components",
     [:children, :id, :n_clicks],
     [Symbol("data-"), Symbol("aria-")];
     id = 10,
     n_clicks = 1
     )

    @test get_name(test_comp) == "html_div"
    @test get_type(test_comp) == "Div"
    @test get_namespace(test_comp) == "dash_html_components"
    @test test_comp.id == 10
    @test test_comp.n_clicks == 1
    test_comp.var"data-id" = 20
    @test test_comp.var"data-id" == 20

    @test sort(propertynames(test_comp)) == sort([:children, :id, :n_clicks])
    @test isnothing(test_comp.children)

    #anavailable property
    @test_throws ErrorException test_comp.key
    @test_throws ErrorException test_comp.key = 20
    @test_throws ErrorException Component("html_div", "Div", "dash_html_components",
     [:children, :id, :n_clicks],
     [Symbol("data-"), Symbol("aria-")];
     id = 10,
     key = 1
     )
    json = JSON2.write(test_comp)
    res = JSON2.read(json)
    @test keys(res) == (:type, :namespace, :props)
    @test sort(collect(keys(res.props))) == sort([:id, :n_clicks, Symbol("data-id")])
    @test res.props.id == 10
    @test res.props.n_clicks == 1
    @test res.props.var"data-id" == 20
end

@testset "empty wilds" begin
    test_comp = Component("html_div", "Div", "dash_html_components",
     [:children, :id, :n_clicks],
     Symbol[];
     id = 10,
     n_clicks = 1
     )

    @test get_name(test_comp) == "html_div"
    @test get_type(test_comp) == "Div"
    @test get_namespace(test_comp) == "dash_html_components"
    @test test_comp.id == 10
    @test test_comp.n_clicks == 1

    @test_throws ErrorException test_comp.var"data-id" = 20
    @test_throws ErrorException test_comp.aaaaa
end

max_stars_repo_name: waralex/DashBase

max_issues_repo_name: waralex/DashBase

max_forks_repo_name: waralex/DashBase


text:
using Test, Documenter, Hubbard

DocMeta.setdocmeta!(Hubbard, :DocTestSetup, :(using Hubbard); recursive=true)
doctest(Hubbard)

max_stars_repo_name: jburgy/hubbard

max_issues_repo_name: jburgy/hubbard

max_forks_repo_name: jburgy/hubbard


text:
using StatsBase, Statistics, OnlineStats, Dates

function latency_profiler(chin::Channel; chout::Channel, nworkers::Integer, qlower::Real=0.1, qupper::Real=0.9, buffersize::Integer=1000, windowsize::Dates.AbstractTime=Second(10))
    @info "latency_profiler task started"

    # maintain a window of latency samples for each worker
    ws = [MovingTimeWindow(windowsize, valtype=@NamedTuple{comp::Float64, comm::Float64}, timetype=Time) for _ in 1:nworkers]

    # utility function for updating the correct window with a new latency measurement
    function fitwindow(v)
        0 < v.worker <= nworkers || @error "Expected v.worker to be in [1, $nworkers], but it is $(v.worker)"
        isnan(v.comm) || isnan(v.comp) || fit!(ws[v.worker], (v.timestamp, @NamedTuple{comp::Float64,comm::Float64}((v.comp, v.comm))))
        return
    end

    # helper functions for computing the mean and variance over the samples in a window
    buffer = zeros(buffersize)
    function processwindow(w::MovingTimeWindow, key::Symbol)
        key == :comp || key == :comm || throw(ArgumentError("key must be either :comp or :comm, but is $key"))

        # populate the buffer
        i = 0
        n = 0
        for (_, t) in value(w)
            v = t[key]
            if !isnan(v)
                buffer[i+1] = v
                i = mod(i + 1, buffersize)
                n += 1
            end
        end
        n = min(n, buffersize)

        # return NaNs if there are no values
        if n == 0
            return NaN, NaN
        end

        # compute quantile indices
        sort!(view(buffer, 1:n))
        il = max(1, round(Int, n*qlower))
        iu = min(buffersize, round(Int, n*qupper))

        # compute mean and variance over the values between qlower and qupper
        vs = view(buffer, il:iu)
        m = mean(vs)
        v = var(vs, mean=m, corrected=true)
        m, v
    end

    # process incoming latency samples
    while isopen(chin)

        # consume all values currently in the channel
        try
            vin = take!(chin)
            fitwindow(vin)
        catch e
            if e isa InvalidStateException
                break
            else
                rethrow()
            end
        end
        while isready(chin)
            try
                vin = take!(chin)
                fitwindow(vin)
            catch e
                if e isa InvalidStateException
                    break
                else
                    rethrow()
                end
            end
        end

        for i in 1:nworkers

            # for each worker, compute the mean and variance of the values in the window
            comp_mean, comp_var = processwindow(ws[i], :comp)
            comm_mean, comm_var = processwindow(ws[i], :comm)
            if isnan(comp_mean) || isnan(comp_var) || isnan(comm_mean) || isnan(comm_var)
                continue
            end

            # push the computed statistics into the output channel
            vout = @NamedTuple{worker::Int,comp_mean::Float64,comp_var::Float64,comm_mean::Float64,comm_var::Float64}((i,comp_mean,comp_var,comm_mean,comm_var))
            try
                push!(chout, vout)
            catch e
                if e isa InvalidStateException
                    break
                else
                    rethrow()
                end
            end
        end
    end
    @info "latency_profiler task finished"
end

function process_stats(chin::Channel)
    @info "process_stats task started"
    for v in chin
        @info v
    end
    @info "process_stats task finished"
end

function get_test_values(df)
    df = filter(:jobid => (x)->x==1, df)
    df.time, df.compute_latency_worker_4
end

function main(df, jobid=17, worker=4, nworkers=10)

    df = filter(:jobid => (x)->x==jobid, df)
    sort!(df, :iteration)
    ts = df.time
    comps = df[:, "compute_latency_worker_$worker"]
    comms = df[:, "latency_worker_$worker"] .- comps

    nwait = 1
    ps = fill(240, nworkers)
    θs = fill(1/nworkers, nworkers)
    min_processed_fraction = sum(θs ./ ps) * nwait / nworkers

    # channel to push recorded latency into
    latency_channel = Channel{@NamedTuple{worker::Int,timestamp::Time,comp::Float64,comm::Float64}}(Inf)
    stats_channel = Channel{@NamedTuple{worker::Int,comp_mean::Float64,comp_var::Float64,comm_mean::Float64,comm_var::Float64}}(Inf)
    partitions_channel = Channel{@NamedTuple{worker::Int,p::Int}}(Inf)

    t1 = Threads.@spawn latency_profiler(latency_channel; chout=stats_channel, nworkers)
    t2 = Threads.@spawn load_balancer(stats_channel; chout=partitions_channel, ps0=ps, θs, min_processed_fraction, nwait)
    # t2 = Threads.@spawn process_stats(stats_channel)

    # while !istaskstarted(t1) || !istaskstarted(t2)
    #    sleep(1e-3)
    # end

    timestamps = Nanosecond.(round.(Int, ts .* 1e9)) .+ Time(0)
    i = 0
    for i in 1:size(df, 1)
        for j in 1:nworkers
            comp = df[i, "compute_latency_worker_$j"]
            comm = df[i, "latency_worker_$j"] .- comp
            if !isnan(comp) && !isnan(comm)
                v = @NamedTuple{worker::Int,timestamp::Time,comp::Float64,comm::Float64}((j, timestamps[i], comp, comm))
                push!(latency_channel, v)
            end
        end

        while isready(partitions_channel)
            try
                vin = take!(partitions_channel)
                @info "ps[$(vin.worker)] = $(vin.p)"
                ps[vin.worker] = vin.p
            catch e
                if e isa InvalidStateException
                    break
                else
                    rethrow()
                end
            end
        end

        sleep(1e-6)

        # if i > 100
        #     break
        # end
    end

    println("coord finished")

    # wait for all latency values to be processed
    while !istaskfailed(t1) && !istaskfailed(t2) && isready(latency_channel)
        sleep(0.1)
    end

    # close the channels and wait fo the tasks to finish
    close(latency_channel)
    close(stats_channel)
    wait(t1)
    wait(t2)

    return
end