text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ic = "exp"
ks_data = CSV.read((@__DIR__)*"/../../data/errors/ks_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_"*ic*".csv",DataFrame)
ks_data.DG_limit2 = tmp.DG_limit

ks_data_lwr = CSV.read((@__DIR__)*"/../../data/errors/ks_lwr_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_lwr_"*ic*".csv",DataFrame)
ks_data_lwr.DG_limit2 = tmp.DG_limit

ks_data_upr = CSV.read((@__DIR__)*"/../../data/errors/ks_upr_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/ks_upr_"*ic*".csv",DataFrame)
ks_data_upr.DG_limit2 = tmp.DG_limit

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot, :dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
colours = [1;3;2;4;5]
_names = ["DG", "DG lim", "Unif","QBD-RAP", "DG lin lim"]

for (c,col) in enumerate(names(ks_data))
    if c<5
        plot!(log10.(1:2:21),log10.(ks_data[:,col]),
            ribbon=(log10.(ks_data[:,col])-log10.(ks_data_lwr[:,col]),log10.(ks_data_upr[:,col])-log10.(ks_data[:,col])),
            label=_names[c],linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    else
        plot!(log10.((1:2:21).+1),log10.(ks_data[:,col]),
            ribbon=(log10.(ks_data[:,col])-log10.(ks_data_lwr[:,col]),log10.(ks_data_upr[:,col])-log10.(ks_data[:,col])),
            label=_names[c],linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    end
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.215,0.635))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
ic = "exp"
l1_cdf_data = CSV.read((@__DIR__)*"/../../data/errors/l1_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_"*ic*".csv",DataFrame)
l1_cdf_data.DG_limit2 = tmp.DG_limit

l1_cdf_data_lwr = CSV.read((@__DIR__)*"/../../data/errors/l1_lwr_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_lwr_"*ic*".csv",DataFrame)
l1_cdf_data_lwr.DG_limit2 = tmp.DG_limit

l1_cdf_data_upr = CSV.read((@__DIR__)*"/../../data/errors/l1_upr_"*ic*".csv",DataFrame)
tmp = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/errors/l1_upr_"*ic*".csv",DataFrame)
l1_cdf_data_upr.DG_limit2 = tmp.DG_limit
_names = ["DG", "DG lim", "Unif","QBDRAP", "DG lin lim"]

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot, :dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross, :diamond]
colours = [1;3;2;4;5]

for (c,col) in enumerate(names(l1_cdf_data))
    if c < 5
        plot!(log10.(1:2:21),log10.(l1_cdf_data[:,col]),
            ribbon=(log10.(l1_cdf_data[:,col])-log10.(l1_cdf_data_lwr[:,col]),log10.(l1_cdf_data_upr[:,col])-log10.(l1_cdf_data[:,col])),
            label=_names[c],linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    else
        plot!(log10.((1:2:21).+1),log10.(l1_cdf_data[:,col]),
            ribbon=(log10.(l1_cdf_data[:,col])-log10.(l1_cdf_data_lwr[:,col]),log10.(l1_cdf_data_upr[:,col])-log10.(l1_cdf_data[:,col])),
            label=_names[c],linestyle=linestyles_vec[c],
            # marker=markerstyles_vec[c],
            linewidth=2,color=colours[c])
    end
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.215,0.635))
# plot!(ylims=ylims(plot!()).+(-0.18,0.0),xlims=xlims(plot!()).+(-0.13,0.0))
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")

end
p()

function p()
colours = [1;3;2;4]
@fullwidth_plot_defaults()
function ht_cdf(sims)
    n_sims=length(sims.t)
    function F(x,i)
        i_idx = sims.φ.==i
        Fxi = sum(sims.t[i_idx].<=x)/n_sims
        return Fxi
    end
    return F
end
sim_data = CSV.read((@__DIR__)*"/../../data/exp/sims.csv",DataFrame)
data = CSV.read((@__DIR__)*"/../../data/exp/order_21_model_dg1.csv",DataFrame)
T_idx = only(findall(data.t.==10.0))
sim_cdf = CSV.read((@__DIR__)*"/../../data/exp/sims_cdf_evaluated.csv",DataFrame).phase_2[1:T_idx]
linestyles_vec = [:solid,:dash,:dashdot, :dot]
plot(data.t[1:T_idx],data.phase_2[1:T_idx],label="DG - 21",color=colours[1],linestyle=linestyles_vec[1],linewidth=2)
data = CSV.read((@__DIR__)*"/../../data/exp/order_21_model_qbdrap4.csv",DataFrame)
plot!(data.t[1:T_idx],data.phase_2[1:T_idx],label="QBDRAP - 21",color=colours[4],linestyle=linestyles_vec[4],linewidth=2)
data = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/exp/order_21_model_dg1.csv",DataFrame)
plot!(data.t[1:T_idx],data.phase_2[1:T_idx],label="DG lin lim - 21",color=colours[4],linestyle=linestyles_vec[4],linewidth=2)
plot!(data.t[1:T_idx],sim_cdf,label="Sim",color=5,linestyle=:dash,linewidth=2)
plot!(xlabel="t")
plot!(ylabel="Probability")
plot!(title="CDF of first hitting time in phase 2")
plot!(legend=(0.7,0.35))
savefig((@__DIR__)*"/cdf_order21DG_and_sims.pdf")
end
p()

function p()
colours = [1;3;2;4;5]
@fullwidth_plot_defaults()
data_dg = CSV.read((@__DIR__)*"/../../data/exp/order_21_model_dg1.csv",DataFrame)
data_qbdrap = CSV.read((@__DIR__)*"/../../data/exp/order_21_model_qbdrap4.csv",DataFrame)
data_dg2 = CSV.read("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_linearDGwLimiter/hitting_times/data/exp/order_21_model_dg1.csv",DataFrame)
linestyles_vec = [:solid,:dash,:dashdot, :dot]
h = data_dg.t[2]-data_dg.t[1]
plot((data_dg.t[1:end-1]+data_dg.t[2:end])/2.0,diff(data_dg.phase_1)./h,label="DG - Phase 1",color=colours[1],linestyle=linestyles_vec[1],linewidth=2)
plot!((data_dg.t[1:end-1]+data_dg.t[2:end])/2.0,diff(data_dg.phase_2)./h,label="DG - Phase 2",color=colours[1],linestyle=linestyles_vec[2],linewidth=2)
plot!((data_dg2.t[1:end-1]+data_dg2.t[2:end])/2.0,diff(data_dg2.phase_1)./h,label="DG lin lim - Phase 1",color=colours[5],linestyle=linestyles_vec[1],linewidth=2)
plot!((data_dg2.t[1:end-1]+data_dg2.t[2:end])/2.0,diff(data_dg2.phase_2)./h,label="DG lin lim - Phase 2",color=colours[5],linestyle=linestyles_vec[2],linewidth=2)
plot!((data_qbdrap.t[1:end-1]+data_qbdrap.t[2:end])/2.0,diff(data_qbdrap.phase_1)./h,label="QBD-RAP - Phase 1",color=colours[4],linestyle=linestyles_vec[1],linewidth=2)
plot!((data_qbdrap.t[1:end-1]+data_qbdrap.t[2:end])/2.0,diff(data_qbdrap.phase_2)./h,label="QBD-RAP - Phase 2",color=colours[4],linestyle=linestyles_vec[2],linewidth=2)
plot!(xlabel="t")
plot!(ylabel="Probability")
plot!(title="First hitting time PDFs")
plot!(xlims=(0.0,4/3))
plot!(legend=(0.77,0.95))
savefig((@__DIR__)*"/pdf_order21.pdf")
end
p()


text:
"""
    initialize(::Type{TransmissionRates},
                    states::DiseaseStates,
                    pop::Population,
                    rf::RiskFunctions{T},
                    rp::RiskParameters{T}) where T <: DiseaseStateSequence

Initialize `TransmissionRates` for a `Population` with specified `DiseaseStates`.
`TransmissionRates` are an array storage of rates of disease transmission between individuals.
`TransmissionRates` are specific to the defined `RiskFunctions`, and their parametrization in `RiskParameters`.
"""
function initialize(::Type{TransmissionRates},
                    states::DiseaseStates,
                    pop::Population,
                    rf::RiskFunctions{T},
                    rp::RiskParameters{T}) where T <: DiseaseStateSequence
  n_ids = length(states)
  tr = TransmissionRates(n_ids)
  for i in findall(states .== Ref(State_S))
    # External exposure
    #tr.external[i] = rf.susceptibility(rp.susceptibility, pop, i) * rf.sparks(rp.sparks, pop, i)
    tr.external[i] = rf.sparks(rp.sparks, pop, i)
    # Internal exposure
    susceptibility = rf.susceptibility(rp.susceptibility, pop, i)
    for k in findall(states .== Ref(State_I))
      tr.internal[k, i] = susceptibility *
                          rf.infectivity(rp.infectivity, pop, i, k) *
                          rf.transmissibility(rp.transmissibility, pop, k)
    end
  end
  @debug "Initialization of $T TransmissionRates complete" ∑external = sum(tr.external) ∑internal = sum(tr.internal)
  return tr
end

"""
    function initialize(::Type{EventRates},
                        tr::TransmissionRates,
                        states::DiseaseStates,
                        pop::Population,
                        rf::RiskFunctions{T},
                        rp::RiskParameters{T}) where T <: DiseaseStateSequence

Initialize `EventRates`, which represent the rates for all competing events in an epidemic (e.g. removal of individual x; exposure of individual y)
"""
function initialize(::Type{EventRates},
                    tr::TransmissionRates,
                    states::DiseaseStates,
                    pop::Population,
                    rf::RiskFunctions{T},
                    rp::RiskParameters{T}) where T <: DiseaseStateSequence
  n_ids = length(states)
  rates = EventRates{T}(n_ids)
  for i = 1:n_ids
    if states[i] == State_S
      if State_E ∈ T
        rates.exposure[i] = sum(tr, i)
      elseif State_E ∉ T
        rates.infection[i] = sum(tr, i)
      end
    elseif states[i] == State_E
      rates.infection[i] = rf.latency(rp.latency, pop, i)
    elseif states[i] == State_I
      if State_R ∈ T
        rates.removal[i] = rf.removal(rp.removal, pop, i)
      end
    end
  end
  @debug "Initialization of $T EventRates complete" rates = [sum(rates[k]) for k in convert(DiseaseStates, T)[2:end]]
  return rates
end

ext: jl


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
root = 0
rank = MPI.Comm_rank(comm)
sz   = MPI.Comm_size(comm)
isroot = MPI.Comm_rank(comm) == root

for T in Base.uniontypes(MPI.MPIDatatype)

    # Allocating
    A = ArrayType{T}(fill(T(rank+1), 4))
    C = MPI.Gather(A, root, comm)
    if isroot
        @test C isa ArrayType{T}
        @test C == ArrayType{T}(repeat(1:sz, inner=4))
    end

    # Allocating, object
    C = MPI.Gather(T(rank+1), root, comm)
    if isroot
        @test C isa Vector{T}
        @test C == Vector{T}(1:sz)
    end

    # Allocating, explicit length
    A = ArrayType{T}([MPI.Comm_rank(comm) + 1, 0])
    C = MPI.Gather(A, 1, root, comm)
    if isroot
        @test C isa ArrayType{T}
        @test C == ArrayType{T}(1:MPI.Comm_size(comm))
    end

    # Allocating, view
    A = ArrayType{T}([MPI.Comm_rank(comm) + 1, 0])
    C = MPI.Gather(view(A, 1:1), root, comm)
    if isroot
        @test C isa ArrayType{T}
        @test C == ArrayType{T}(1:MPI.Comm_size(comm))
    end

    # Non Allocating
    A = ArrayType{T}(fill(T(rank+1), 4))
    C = ArrayType{T}(undef, 4sz)
    MPI.Gather!(A, C, length(A), root, comm)
    if isroot
        @test C == ArrayType{T}(repeat(1:sz, inner=4))
    end

    # In_place
    A = ArrayType{T}(fill(T(rank+1), 4))
    if isroot
        A = ArrayType{T}(fill(T(rank+1), 4*sz))
    end
    if root == MPI.Comm_rank(comm)
        MPI.Gather!(nothing, A, 4, root, comm)
    else
        MPI.Gather!(A, nothing, 4, root, comm)
    end
    if isroot
        @test A == ArrayType{T}(repeat(1:sz, inner=4))
    end

end

MPI.Finalize()
@test MPI.Finalized()


text:
function D3Trees.D3Tree(D::DESPOTAlphaTree; title="DESPOTAlpha Tree", kwargs...)
    lenb = D.b
    lenba = D.ba
    len = lenb + lenba
    children = Vector{Vector{Int}}(undef, len)
    text = Vector{String}(undef, len)
    tt = fill("", len)
    link_style = fill("", len)
    for b in 1:lenb
        as = retrieve_as(D, b)
        obs = b == 1 ? "<root>" : string(D.obs[D.as_map[as]][b-D.ba_children[D.parent[b]][1]+1])
        children[b] = D.children[b] .+ lenb
        text[b] = @sprintf("""
                           o:%s prob:%6.2f
                           l:%6.2f, u:%6.2f""",
                           obs,
                           D.obs_prob[b],
                           D.l[b],
                           D.u[b],
                          )
        tt[b] = """
                o: $(obs)
                prob: $(D.obs_prob[b])
                l: $(D.l[b])
                u: $(D.u[b])
                $(length(D.children[b])) children
                """
        link_width = 2.0
        link_style[b] = "stroke-width:$link_width"
        for ba in D.children[b]
            link_style[ba+lenb] = "stroke-width:$link_width"
        end

        push!(as, zero(eltype(as)))
        for ba in D.children[b]
            as[length(as)] = D.ba_action[ba]
            as_ba = D.as_map[as]
            children[ba+lenb] = D.ba_children[ba]
            text[ba+lenb] = @sprintf("""
                                     a:%s r:%6.2f |ϕ|:%2d
                                     l:%6.2f, u:%6.2f""",
                                     D.ba_action[ba], D.weights[b]' * D.r[as_ba], length(D.particles[as_ba]),
                                     D.ba_l[ba], D.ba_u[ba],
                                     )
            tt[ba+lenb] = """
                          a: $(D.ba_action[ba])
                          r: $(D.weights[b]' * D.r[as_ba])
                          |ϕ|:$(length(D.particles[as_ba]))
                          l: $(D.ba_l[ba])
                          u: $(D.ba_u[ba])
                          $(length(D.ba_children[ba])) children
                          """
        end

    end
    return D3Tree(children;
                  text=text,
                  tooltip=tt,
                  link_style=link_style,
                  title=title,
                  kwargs...
                 )
end

Base.show(io::IO, mime::MIME"text/html", D::DESPOTAlphaTree) = show(io, mime, D3Tree(D))
Base.show(io::IO, mime::MIME"text/plain", D::DESPOTAlphaTree) = show(io, mime, D3Tree(D))

function retrieve_as(D::DESPOTAlphaTree{S,A}, b::Int) where {S,A}
    as = Vector{A}(undef, D.Delta[b])
    while b != 1
        bp = b
        ba = D.parent[bp]
        b = D.ba_parent[ba]
        as[D.Delta[bp]] = D.ba_action[ba]
    end
    return as
end


text:
## from base/boot.jl:
#
# immutable ASCIIString <: DirectIndexString
#     data::Array{Uint8,1}
# end
#

## required core functionality ##

endof(s::ASCIIString) = length(s.data)
getindex(s::ASCIIString, i::Int) = (x=s.data[i]; x < 0x80 ? char(x) : '\ufffd')

## overload methods for efficiency ##

sizeof(s::ASCIIString) = sizeof(s.data)

getindex(s::ASCIIString, r::Vector) = ASCIIString(getindex(s.data,r))
getindex(s::ASCIIString, r::UnitRange{Int}) = ASCIIString(getindex(s.data,r))
getindex(s::ASCIIString, indx::AbstractVector{Int}) = ASCIIString(s.data[indx])
search(s::ASCIIString, c::Char, i::Integer) = c < 0x80 ? search(s.data,uint8(c),i) : 0
rsearch(s::ASCIIString, c::Char, i::Integer) = c < 0x80 ? rsearch(s.data,uint8(c),i) : 0

function string(c::ASCIIString...)
    if length(c) == 1
        return c[1]
    end
    n = 0
    for s in c
        n += length(s.data)
    end
    v = Array(Uint8,n)
    o = 1
    for s in c
        ls = length(s.data)
        unsafe_copy!(v, o, s.data, 1, ls)
        o += ls
    end
    ASCIIString(v)
end

function ucfirst(s::ASCIIString)
    if length(s) > 0 && 'a' <= s[1] <= 'z'
        t = ASCIIString(copy(s.data))
        t.data[1] -= 32
        return t
    end
    return s
end
function lcfirst(s::ASCIIString)
    if length(s) > 0 && 'A' <= s[1] <= 'Z'
        t = ASCIIString(copy(s.data))
        t.data[1] += 32
        return t
    end
    return s
end

function uppercase(s::ASCIIString)
    d = s.data
    for i = 1:length(d)
        if 'a' <= d[i] <= 'z'
            td = copy(d)
            for j = i:length(td)
                if 'a' <= td[j] <= 'z'
                    td[j] -= 32
                end
            end
            return ASCIIString(td)
        end
    end
    return s
end
function lowercase(s::ASCIIString)
    d = s.data
    for i = 1:length(d)
        if 'A' <= d[i] <= 'Z'
            td = copy(d)
            for j = i:length(td)
                if 'A' <= td[j] <= 'Z'
                    td[j] += 32
                end
            end
            return ASCIIString(td)
        end
    end
    return s
end

reverse(s::ASCIIString) = ASCIIString(reverse(s.data))

## outputing ASCII strings ##

print(io::IO, s::ASCIIString) = (write(io, s);nothing)
write(io::IO, s::ASCIIString) = write(io, s.data)

## transcoding to ASCII ##

ascii(x) = convert(ASCIIString, x)
convert(::Type{ASCIIString}, s::ASCIIString) = s
convert(::Type{ASCIIString}, s::UTF8String) = ascii(s.data)
convert(::Type{ASCIIString}, a::Array{Uint8,1}) = is_valid_ascii(a) ? ASCIIString(a) : error("invalid ASCII sequence")
function convert(::Type{ASCIIString}, a::Array{Uint8,1}, invalids_as::ASCIIString)
    l = length(a)
    idx = 1
    iscopy = false
    while idx <= l
        (a[idx] < 0x80) && (idx +=1; continue)
        !iscopy && (a = copy(a); iscopy = true)
        endn = idx
        while endn <= l
            (a[endn] < 0x80) && break
            endn += 1
        end
        (endn > idx) && (endn -= 1)
        splice!(a, idx:endn, invalids_as.data)
        l = length(a)
    end
    convert(ASCIIString, a)
end
convert(::Type{ASCIIString}, s::String) = ascii(bytestring(s))

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
L = 200
labels = rand(Bool, L);
scores = add_noise.(labels, 0.6)

B = BinaryMetrics(labels, scores)
println(B)

BinaryMetrics([0;1], [0.4;0.1])
@test_throws ErrorException BinaryMetrics([-1;1], [0.4;0.1])
@test_throws ErrorException BinaryMetrics([0;1;1], [0.4;0.1])

B = BinaryMetrics([0;0;1;1], [0.4;0.7;0.4;1])
@test tp(B,0.5) == 1
@test tn(B,0.5) == 1
@test fp(B,0.5) == 1
@test fn(B,0.5) == 1
@test cm(B,0.5) == ones(2,2)

L = 200
labels = rand(Bool, L);
scores = add_noise.(labels, 0.6)
B = BinaryMetrics(labels, scores)
OP = [0.2;0.5;0.7]
tps = tp(B, OP)
tns = tn(B, OP)
fps = fp(B, OP)
fns = fn(B, OP)
@test cm(B,OP) == (tps, fps, fns, tns)

@test isapprox(tpr(B, OP) , 1 .- fnr(B,OP); atol=1e-8)
@test isapprox(tpr(B, OP) , tps./(tps.+fns); atol=1e-8)
@test isapprox(tnr(B, OP) , 1 .- fpr(B,OP); atol=1e-8)
@test isapprox(tnr(B, OP) , tns./(tns.+fps); atol=1e-8)
@test isapprox(fnr(B, OP) , 1 .- tpr(B,OP); atol=1e-8)
@test isapprox(fnr(B, OP) , fns./(fns.+tps); atol=1e-8)
@test isapprox(fpr(B, OP) , 1 .- tnr(B,OP); atol=1e-8)
@test isapprox(fpr(B, OP) , fps./(fps.+tns); atol=1e-8)
@test isapprox(ppv(B, OP) , 1 .- fdr(B, OP); atol=1e-8)
@test isapprox(npv(B, OP) , 1 .- false_omission_rate(B, OP); atol=1e-8)
@test csi(B, OP) == tp(B,OP) ./ (tp(B,OP) + fn(B,OP) +fp(B,OP))
@test acc(B, OP) == (tp(B,OP)+tn(B,OP)) ./ ( tp(B,OP) + tn(B,OP)
                                            +fp(B,OP) + fn(B,OP) )
F1 = (2 .* tp(B,OP)) ./ (2 .* tp(B,OP) .+ fp(B,OP) .+ fn(B,OP))
@test isapprox(f1(B,OP), F1 ; atol=1e-8)
MCC = (tp(B,OP).*tn(B,OP).-fp(B,OP).*fn(B,OP))./ # <- determinant of CM
sqrt.((tp(B,OP).+fp(B,OP)).*(tp(B,OP).+fn(B,OP)).*(tn(B,OP).+fp(B,OP)).*(tn(B,OP).+fn(B,OP)))
@test mcc(B,OP) == MCC
@test eer(B) == (0.17,100)


text:
module ONCRPC

export make_rpc

@enum msg_type CALL=0 REPLY=1
@enum reply_stat MSG_ACCEPTED=0 MSG_REKECTED=1
@enum accept_stat SUCCESS=0 PROG_UNAVAIL=1 PROG_MISMATCH=2 PROC_UNAVAIL=3 GARBAGE_ARGS=4 SYSTEM_ERR=5
@enum reject_stat RPC_MISMATCH=0 AUTH_ERROR=1
@enum auth_stat AUTH_OK=0 AUTH_BADCRED=1 AUTH_REJECTEDCRED=2 AUTH_BADVERF=3 AUTH_REJECTEDVERF=4 AUTH_TOOWEAK=5 AUTH_INVALIDRESP=6 AUTH_FAILED=7
@enum auth_flavor AUTH_NONE=0 AUTH_SYS=1 AUTH_SHORT=2
struct Channel{T<:IO}
    io::T
end

function read_reply(c::Channel{TCPSocket})
    reply_buf = IOBuffer()
    while true
        fragment_header = ntoh(read(c.io, UInt32))
        is_last = (fragment_header & (UInt32(1) << 31)) != 0
        length = fragment_header & ~(UInt32(1) << 31)
        write(reply_buf, read(c.io, length))
        if is_last
            break
        end
    end
    seekstart(reply_buf)
end

const max_fragment_size = ~(UInt32(1) << 31)
function frame_and_send(c::Channel{TCPSocket}, data)
    data = take!(data)
    length = sizeof(data)
    offset = 1
    while length > max_fragment_size
        write(c.io, max_fragment_size)
        write(c.io, data[offset:(offset+max_fragment_size-1)])
        length -= max_fragment_size
    end
    buf = IOBuffer()
    framing_length = UInt32(length-offset+1) | (UInt32(1) << 31)
    write(buf, hton(framing_length))
    write(buf, data[offset:end])
    write(c.io, take!(buf))
end

function make_call(c::Channel, prog::UInt32, vers::UInt32, proc::UInt32, arg_data::Vector{UInt8}=UInt8[], xid::UInt32=UInt32(0x1))
    msg = IOBuffer()
    write(msg, hton(xid))
    write(msg, hton(UInt32(CALL)))
    write(msg, hton(UInt32(2)))       # rpcvers
    write(msg, hton(prog))
    write(msg, hton(vers))
    write(msg, hton(proc))
    # Authentication - only none for now
    write(msg, AUTH_NONE); write(msg, UInt32(0));
    write(msg, AUTH_NONE); write(msg, UInt32(0));
    write(msg, arg_data)
    frame_and_send(c, msg)
    read_reply(c)
end
make_call(c::Channel, prog::Integer, vers::Integer, proc::Integer, arg_data=UInt8[], xid::Integer=0x1) =
    make_call(c, UInt32(prog), UInt32(vers), UInt32(proc), arg_data, UInt32(xid))

function marshal_data(buf::IOBuffer, T::Type, arg)
    if T <: Integer
        # All small integers are promoted to 4 bytes
        if sizeof(T) <= sizeof(UInt32)
            write(buf, hton(arg % UInt32))
        else
            write(buf, hton(arg))
        end
    elseif T <: AbstractString || T == Vector{UInt8}
        size = UInt32(sizeof(arg))
        write(buf, hton(size))
        write(buf, Vector{UInt8}(arg))
        # Pad to 4 byte boundary
        while (size % 4) != 0
            write(buf, UInt8(0))
            size += 1
        end
    else
        for field in fieldnames(T)
            marshal_data(buf, fieldtype(T, field), getfield(arg, field))
        end
    end
    buf
end
marshal_data(t::Type, arg) = take!(marshal_data(IOBuffer(), t, arg))

function demarshal_data(T, reply_buf)
    if T <: Integer
        # All small integers are promoted to 4 bytes
        return if sizeof(T) <= sizeof(UInt32)
            ntoh(read(reply_buf, UInt32)) % T
        else
            ntoh(read(reply_buf, T))
        end
    elseif T <: AbstractString || T == Vector{UInt8}
        ssize = ntoh(read(reply_buf, UInt32))
        str = read(reply_buf, ssize)
        while (ssize % 4) != 0
            read(reply_buf, UInt8)
            ssize += 1
        end
        return T <: AbstractString ? String(str) : str
    else
        args = Any[]
        for field in fieldnames(T)
            push!(args, demarshal_data(fieldtype(T, field), reply_buf))
        end
        return T(args...)
    end
end

function parse_response(T::Type, reply_buf, xid)
    reply_xid = ntoh(read(reply_buf, UInt32))
    (xid == reply_xid) || error("Mismatched XID in reply")
    mt = msg_type(ntoh(read(reply_buf, Int32)))
    (mt == REPLY) || error("Not a REPLY")
    rt = reply_stat(ntoh(read(reply_buf, Int32)))
    (rt == MSG_ACCEPTED) || error("Message was rejected")
    verft = auth_flavor(ntoh(read(reply_buf, UInt32)))
    (verft == AUTH_NONE) || error("Authentication not yet supported")
    verfs = read(reply_buf, UInt32)
    as = accept_stat(ntoh(read(reply_buf, UInt32)))
    (as == SUCCESS) || error("Query failed")
    demarshal_data(T, reply_buf)
end

function make_rpc(prog, vers, proc, ret_type, param_type)
    function(c::Channel, arg)
        arg::param_type
        xid = rand(UInt32)
        reply_buf = make_call(c, prog, vers, proc, marshal_data(param_type, arg), xid)
        parse_response(ret_type, reply_buf, xid)
    end
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ a29c2ea0-6ff4-11eb-2954-372c6727bee2
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.Registry.update()

	Pkg.add("Plots")
	using Plots
end

# ╔═╡ 724acb50-6db3-11eb-07a8-7d048b8e4b3e
md"# Introduction to Elementary Cellular Automata"

# ╔═╡ 9464f530-6db3-11eb-0629-9f05df1979bd
md"
Cellular automata are fascination \"mathematical machines\" which, despite normally being associated with having relatively simple rules, can still produce highly complex output. This point is stressed over and over in Stephen Wolfram's *A New Kind of Science*. This is how Wikipedia defines them:
> A cellular automaton is a discrete model of computation studied in automata theory.
In this notebook we will explore some simple cellular automata, called the *Elementary Cellular Automata*.
"

# ╔═╡ f3f9f3b0-6db3-11eb-1446-cfe467e42fee
md"## Defining Rules"

# ╔═╡ fdc546b0-6db3-11eb-1a4c-514e32ae4a27
md"Put simply, a cellular automata takes some simple rule and applies it over and over on an infinite string of numbers. We will be taking a look at 1 dimensional cellular automata, which is to say there is one spatial dimension and one temporal dimension. This can easily be visualized in a 2D image unlike higher-dimensional automata. Chances are you have heard of Conway's *Game of Life*, which is an example of a two dimensional cellular automaton."

# ╔═╡ 87bd67d0-6db4-11eb-0931-f9a3ec94cd70
md"Let's define some rules as a dictionary of associations"

# ╔═╡ b3fdfd00-6db4-11eb-3fa3-79d5911c567d
rule = Dict(
	0b000 => 0,
	0b001 => 1,
	0b010 => 1,
	0b011 => 1,
	0b100 => 1,
	0b101 => 0,
	0b110 => 0,
	0b111 => 0
)

# ╔═╡ f21ba14e-6db4-11eb-382a-87bc07250f65
md"For our rule, each cell can have a value of either a $0$ or a $1$. To determine what the cell's next value will be, our rule looks at the value of its left and right neighbors, where the left neighbor is the most significant digit and the right neighbor is the least significant digit, and associates that configuration, which has $2^3=8$ different combinations, to a new value of either $1$ or $0$."

# ╔═╡ 7895f0e0-6db6-11eb-1d6e-b12880ebd5b9
md"## Executing Rules

To start out we can define a recursive function which will return a list of lists, where each list represents the state of the automaton at a particular timestep."

# ╔═╡ 0ab21310-6db6-11eb-00b5-1f2bd2203c79
function value_at(state, pos)
	if pos < 1 || pos > length(state)
		return 0
	end
	return state[pos]
end

# ╔═╡ 656f4030-6db5-11eb-0360-bd9ac3ac6b45
function automaton(initial_state, steps, rule)
	# Base case
	if steps <= 1
		return [initial_state]
	end

	# Otherwise...
	previous_states = automaton(initial_state, steps - 1, rule)
	ps = previous_states[end]
	next_state = [rule[value_at(ps, i-1) << 2 | value_at(ps, i) << 1 | value_at(ps, i+1)] for i = 1:length(initial_state)]

	[previous_states..., next_state]
end

# ╔═╡ 232031c0-6db6-11eb-2f48-81afe98ccac3
md"Despite not being the most efficient it will do for smaller step counts. Let's try it out for our rule for 50 steps"

# ╔═╡ 71d85310-6db6-11eb-0728-f34323017fe1
history = automaton([i == 51 ? 1 : 0 for i = 1:101], 50, rule)

# ╔═╡ aefbdff0-6db6-11eb-0441-a924d0f567c1
history_mat = [history[i][j] for i=1:length(history), j=1:length(history[1])]

# ╔═╡ ea3c54f0-6db6-11eb-1312-3daa47f3fa91
md"## Rendering the Results"

# ╔═╡ f018e190-6db6-11eb-3eef-c984453d36d9
md"Now that we have a matrix of values we can make use of the Plots package to make a plot displaying it."

# ╔═╡ 06c775f0-6db7-11eb-2603-4942353b201e
heatmap(1 .- history_mat, color = :greys, yflip = true, size=(640, 320))

# ╔═╡ 8bb7e6a0-6db7-11eb-3491-039358de17d7
md"Wow! Interpretting this plot is fairly easy. Black is a $1$ and white is a $0$, and time flows downwards. It's rather remarkable how starting from just a single $1$ and the rest $0$'s at $t=1$ the behavior of this cellular automata still appears quite complex! In fact, the center column of this particular rule has been shown to exhibit many of the qualities of a pseudo-random number generator.

Just for fun let's run this automaton for 500 steps just to see what happens if we let it run for a little longer."

# ╔═╡ dfe4d210-6db7-11eb-1425-732f8753c0f5
longer_history = automaton([i == 501 ? 1 : 0 for i = 1:1001], 500, rule)

# ╔═╡ f4bd6df0-6db7-11eb-15a0-c52c636db094
longer_history_mat = [longer_history[i][j] for i=1:length(longer_history), j=1:length(longer_history[1])]

# ╔═╡ f9fe56d0-6db7-11eb-1a02-6133a37e6426
heatmap(1 .- longer_history_mat, color = :greys, yflip = true, size=(640, 320))

# ╔═╡ 065d25f2-6db8-11eb-0ab4-95151bc3554f
md"It's a little tough to see the details since our image isn't high enough resolution to see each and every state value, but it is still fairly clear that the pseudo-random behavior does not settle down in the first 500 steps. However when viewing this many steps it begins to become clear that patterns, or rather regularities, begin to appear moving through space and time at the speed of causality on the left side."

# ╔═╡ 4117f030-6db8-11eb-0df5-4d606d3faf6b
md"## So what??

Beyond being an interesting thing to look at, what else are these automata useful for? The Wikipedia definition of what a cellular automata are mentioned that they are a **model of computation**. Put simply this means that fundamentally cellular automata are performing computations over time. In fact it has been shown that 1 particular rule using only nearest-neighbors and two colors in one dimension (just like our rule) is **Turing-complete**, which means it can do any computation that the computer you're viewing this on can!

I encourage you now to go back and experiment with some of the rule values on the right-hand side. Don't change the left since all 8 of the possible 3 bit combinations have to be present for the automata to be guarenteed to have a valid next state. Also, try using the rule below, which has been proven to be Turing-complete!"

# ╔═╡ caf3afb0-6db8-11eb-39de-852bb7609cb3
turing_complete_rule = Dict(
	0b000 => 0,
	0b001 => 1,
	0b010 => 1,
	0b011 => 1,
	0b100 => 0,
	0b101 => 1,
	0b110 => 1,
	0b111 => 0
)

# ╔═╡ Cell order:
# ╟─724acb50-6db3-11eb-07a8-7d048b8e4b3e
# ╟─9464f530-6db3-11eb-0629-9f05df1979bd
# ╟─f3f9f3b0-6db3-11eb-1446-cfe467e42fee
# ╟─fdc546b0-6db3-11eb-1a4c-514e32ae4a27
# ╟─87bd67d0-6db4-11eb-0931-f9a3ec94cd70
# ╠═b3fdfd00-6db4-11eb-3fa3-79d5911c567d
# ╟─f21ba14e-6db4-11eb-382a-87bc07250f65
# ╟─7895f0e0-6db6-11eb-1d6e-b12880ebd5b9
# ╠═0ab21310-6db6-11eb-00b5-1f2bd2203c79
# ╠═656f4030-6db5-11eb-0360-bd9ac3ac6b45
# ╟─232031c0-6db6-11eb-2f48-81afe98ccac3
# ╠═71d85310-6db6-11eb-0728-f34323017fe1
# ╠═aefbdff0-6db6-11eb-0441-a924d0f567c1
# ╟─ea3c54f0-6db6-11eb-1312-3daa47f3fa91
# ╟─f018e190-6db6-11eb-3eef-c984453d36d9
# ╠═a29c2ea0-6ff4-11eb-2954-372c6727bee2
# ╠═06c775f0-6db7-11eb-2603-4942353b201e
# ╟─8bb7e6a0-6db7-11eb-3491-039358de17d7
# ╠═dfe4d210-6db7-11eb-1425-732f8753c0f5
# ╠═f4bd6df0-6db7-11eb-15a0-c52c636db094
# ╠═f9fe56d0-6db7-11eb-1a02-6133a37e6426
# ╟─065d25f2-6db8-11eb-0ab4-95151bc3554f
# ╟─4117f030-6db8-11eb-0df5-4d606d3faf6b
# ╠═caf3afb0-6db8-11eb-39de-852bb7609cb3

max_stars_repo_name: ctrekker/pluto-static-render-tutorial

max_issues_repo_name: ctrekker/pluto-static-render-tutorial

max_forks_repo_name: ctrekker/pluto-static-render-tutorial


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
           Product1pBasis, Scal1pBasis
using Random: shuffle
using ACEbase.Testing: fdtest, print_tf

##

maxdeg = 5
r0 = 1.0
rcut = 3.0
maxorder = 3
Bsel = SimpleSparseBasis(maxorder, maxdeg)

trans = PolyTransform(1, r0)
Pk = ACE.scal1pbasis(:x, :k, maxdeg, trans, rcut)
RnYlm = ACE.Utils.RnYlm_1pbasis()

B1p = RnYlm * Pk
ACE.init1pspec!(B1p, Bsel)
length(B1p)

##

PosScalState{T} = ACE.State{(:rr, :x), Tuple{SVector{3, T}, T}}

Base.promote_rule(::Union{Type{S}, Type{PosScalState{S}}},
             ::Type{PosScalState{T}}) where {S, T} =
      PosScalState{promote_type(S, T)}

X = rand(PosScalState{Float64})
cfg = ACEConfig([ rand(PosScalState{Float64}) for _=1:10 ])

Rn = B1p.bases[1]
Ylm = B1p.bases[2]
Pk = B1p.bases[3]

##

@info("Checking correct differentiation w.r.t. known and unknown symbols")
for (B, sym, notsym) in zip( (Pk, Rn, Ylm),
                             (:x, :rr, :rr),
                             (:rr, :x, :x) )
   dB1 = evaluate_d(B, X)
   @which evaluate_d(B, X, sym)
   dB2 = evaluate_d(B, X, sym)
   dB3 = evaluate_d(B, X, notsym)
   dB4 = evaluate_d(B, X, :bob)
   print_tf(@test(dB1 == dB2))
   print_tf(@test(all(iszero, norm.(dB3))))
   print_tf(@test(all(iszero, norm.(dB4))))
end
println()

##

@info("Check how the product 1p basis handles a partial derivative")
dB = evaluate_d(B1p, X)
dB_x = evaluate_d(B1p, X, :x)
dB_rr = evaluate_d(B1p, X, :rr)
println(@test dB ≈ dB_x + dB_rr)

##

@info("Test partial derivatives for a symmetric basis")
basis = ACE.SymmetricBasis(ACE.Invariant(), B1p, Bsel)
evaluate(basis, cfg)
dB = evaluate_d(basis, cfg)
dB_x = evaluate_d(basis, cfg, :x)
dB_rr = evaluate_d(basis, cfg, :rr)
println(@test( dB ≈ dB_x + dB_rr ))

##

# try out a chainrule?
using Zygote
import ChainRules: rrule, NoTangent, ZeroTangent
using ACE: LinearACEModel, evaluate

c = randn(length(basis)) ./ (1:length(basis)).^2
model = LinearACEModel(basis, c)

# make up some features we can feed into the x variable.
function x_features(Rs)
   f(r) = exp(- r)
   Xi = [ sum(f(norm(Rs[i] - Rs[j])) for j = 1:length(Rs))
          for i = 1:length(Rs) ]
   Xs = [ ACE.State(rr = rr, x = x) for (rr, x) in zip(Rs, Xi) ]
   return ACEConfig(Xs)
end

function rrule(::typeof(x_features), Rs)
   f(r) = exp(- r)
   df(r) = - exp(-r)
   function x_features_pullback(dP, Rs)
      N = length(Rs)
      g = [ dP[i].rr for i = 1:N ]
      for i = 1:N, j = 1:N
         dxi = dP[i].x
         if j != i
            rr_ij = Rs[j] - Rs[i]
            r_ij = norm(rr_ij)
            g[j] += dxi * df(r_ij) * (rr_ij / r_ij)
            g[i] -= dxi * df(r_ij) * (rr_ij / r_ij)
         end
      end
      return NoTangent(), g
   end
   return x_features(Rs), dp -> x_features_pullback(dp, Rs)
end

eval_model(Rs) = ACE.val(evaluate( model, x_features(Rs) ))

Rs = 2.5 * randn(SVector{3, Float64}, 10)
eval_model(Rs)

##

Zygote.gradient(eval_model, Rs)[1]

# Us = randn(SVector{3, Float64}, length(Rs))

__floats(Rs) = collect(reinterpret(Float64, Rs))
__vecs(x) = collect(reinterpret(SVector{3, Float64}, x))
x0 = __floats(Rs)
F = x -> eval_model(__vecs(x))
dF = x -> Zygote.gradient(eval_model, __vecs(x))[1] |> __floats
println(@test all(ACEbase.Testing.fdtest(F, dF, x0)))


text:
using HSA
include("helpers.jl")
using FactCheck

# Check if the packet headers in the buffers have been
# initialized to "invalid" according to the spec
function check_packetheaders(q)
    for i = 1:q.size
        hdr = HSA.unsafe_convert(PacketHeader, convert(Ptr{Void}, q.base_address + i * 64))
        @fact hdr.typ --> HSA.PacketTypeInvalid "At Index $i"
        if hdr.typ != HSA.PacketTypeInvalid
            break
        end
    end
end

facts("A Queue") do
    rt = HSA.Runtime()
    agents = HSA.all_agents()

    @with_agents context("Can be created as software Queues") do
        a = get_testagent(agents)
        s = HSA.Signal(value = 0)
        r = Ref{Nullable{HSA.Region}}(nothing)
        q = get_softqueue(a,s,r)

        @fact isa(q, Queue) --> true
        @fact q.features & HSA.QueueFeatureAgentDispatch --> HSA.QueueFeatureAgentDispatch

        q_handle = q.doorbell_signal.handle.handle
        s_handle = s.handle.handle
        @fact q_handle --> s_handle "$q_handle/$s_handle"
        @fact q.region.value --> r.x.value

        check_packetheaders(q)
    end

    @with_agents context("Can be created as hardware Queues") do
        a = get_testagent(agents)
        a_minq = HSA.agent_info_queue_min_size(a)
        q = HSA.Queue(a, 0x04)

        @fact q.handle --> x -> x != C_NULL
        @fact q.is_active --> true
        @fact q.typ --> not(nothing)
        @fact q.features --> not((HSA.hsa_queue_feature_t)(0))
        @fact q.base_address --> not(UInt64(0))
        @fact isa(q.doorbell_signal, Signal) --> true
        @fact q.size --> greater_than_or_equal(0x04) "Requested/Min/Actual QueueSize: 4/$a_minq/$(q.size)"
        @fact q.id --> not(nothing)

        # try to see if the signal is actually valid
        HSA.store!(q.doorbell_signal, 1)

        check_packetheaders(q)
    end

    @with_agents context("Can be destroyed") do
        a = get_testagent(agents)
        q = HSA.Queue(a, 0x04)

        finalize(q)

        @fact q.handle --> C_NULL
        @fact q.is_active --> false
    end

    @with_agents context("Can have their error callback invoked") do
        local called_with_q, called_with_s

        a = get_testagent(agents)
        q = HSA.Queue(a, 0x04, error_callback = function(s,q_arg)
            called_with_q = q_arg
            called_with_s = s
        end)

        HSA.queue_err_cb(HSA.HSA_STATUS_ERROR, q.handle, C_NULL)

        @fact called_with_q --> q
        @fact called_with_s --> HSA.HSA_STATUS_ERROR
    end

    @with_agents context("Can have their indexes manipulated") do
        a = get_testagent(agents)
        q = HSA.Queue(a, 0x04)

        @fact HSA.load_read_index(q) --> 0
        @fact HSA.load_read_index(q, Val{HSA.Relaxed}) --> 0

        @fact HSA.load_write_index(q) --> 0
        HSA.store_write_index!(q, UInt64(1))
        @fact HSA.load_write_index(q) --> 1

        @fact HSA.add_write_index!(q, UInt64(1)) --> 1
        @fact HSA.load_write_index(q) --> 2
    end

    @with_agents context("Can be accessed by index") do
        a = get_testagent(agents)
        q = HSA.Queue(a, 0x04)
        p = HSA.AgentDispatchPacket(0x8000)

        q[q.size] = p
        @fact q[0] --> p
        @fact q[q.size] --> p
        q[1] = p
        @fact q[1] --> p
    end

    # currently fails with a segmentation fault when accessing the
    # doorbell signal because the runtime returns a broken software queue
    @with_agents @pending context("Can submit new packets via push") do
        a = get_testagent(agents)
        q = get_softqueue(a)

        p = HSA.AgentDispatchPacket(0x8000)

        push!(q,p)

        @fact HSA.load_write_index(q) --> 1
    end

    @with_agents context("Can be inactivated") do
        a = get_testagent(agents)
        q = HSA.Queue(a, 0x04)

        HSA.queue_inactivate(q)

        @fact q.is_active --> false
    end

    finalize(rt)
end


text:
#!/usr/bin/env julia

if length(ARGS) < 2
    println("Usage: ./royprefit.jl DETX HDF5FILE OUTFILE")
    exit(1)
end

using LinearAlgebra
using NeRCA
using HDF5
using ProgressMeter

function main()
    detx, filename, outfile = ARGS
    if isfile(outfile)
        println("Output file already exists, aborting...")
        exit(1)
    end

    outf = open(outfile, "w")
    write(outf, "group_id,dx,dy,dz,x,y,z,v,t0\n")

    @showprogress 1 for event in NeRCA.MCEventReader(filename, detx)
        hits = calibrate(event.hits, event.calib)
        sort!(hits, by=h->h.t)
        sort!(hits, by=h->h.dom_id)
        NeRCA.count_multiplicities!(hits)

        shits = filter(h->h.multiplicity.count >= 4, hits)
        doms = unique(map(h->h.dom_id, shits))
        track = NeRCA.prefit(shits)
        dir = Direction(normalize(track.dir))
        v = norm(track.dir)
        write(outf, "$(event.info.group_id),$(dir.x),$(dir.y),$(dir.z),$(track.pos.x),$(track.pos.y),$(track.pos.z),$(v),$(track.time)\n")
    end

    close(outf)
end

main()


text:
using PyPlot, Random, ProgressMeter

"""
Train 2-layer relu NN with GD on square loss
INPUT
 m is the number of neurons,
 alpha1 (resp. alpha2): scale of the initialization of the input (resp. output) weights
 ouly_output: train just the output layer (instead of both)
 sym_init: trick to make the output = 0 exactly at initialization
OUTPUT
 Ws: the whole training path
 loss: evolution of loss during training
"""
function twonet_square(X, Y, m, stepsize, niter;alpha1=1.0, alpha2=1.0, weight_decay=0.0, sym_init=false, only_output = false)
    (n,d) = size(X) # n samples in R^d
    # initialize
    W = randn(m, d+1)
    # input weights initialized uniformly on the unit sphere
    W[:,1:d] = W[:,1:d] .* alpha1 ./ sqrt.(sum(W[:,1:d].^2, dims=2))
    # output weights initialized by -1 or 1 with equal proba
    W[1:div(m,2),end] .= alpha2
    W[div(m,2)+1:end,end] .= -alpha2
    W = W
    if sym_init
        @assert iseven(m)
        W[div(m,2)+1:end,1:d]  .=  W[1:div(m,2),1:d]
    end
    Ws    = zeros(m, d+1, niter) # store optimization path
    loss  = zeros(niter) # loss is -log of the empirical risk

    @showprogress 1 "Training neural network..." for iter = 1:niter
        Ws[:,:,iter] = W
        act  =  max.( W[:,1:end-1] * X', 0.0) # activations
        out  =  (1/m) * sum( W[:,end] .* act , dims=1)[:] # predictions of the network
        gradR   = (out .- Y)/n  # size n
        grad_w1 = (W[:,end] .* float.(act .> 0) * ( X .* gradR  )) + weight_decay * W[:,1:end-1] # gradient for input weights
        grad_w2 = act * gradR + weight_decay * W[:,end]  # gradient for output weights
        if only_output
            grad = cat(zeros(m,d), grad_w2, dims=2) # size (m × d+1)
        else
            grad = cat(grad_w1, grad_w2, dims=2) # size (m × d+1)
        end
        loss[iter] = (1/2)*sum( ( out - Y).^2 )/n + (weight_decay/2)*sum(W.^2)/m
        W = W - stepsize * grad
    end

    return Ws, loss
end

X = [-1 1; -1/16 1; 1/16 1; 1 1]
Y = [-1; 1/2; -1/2; 1]
alpha1 = 50
alpha2s=[0 5 10 50]
m = 5000
WW = zeros(m,3,length(alpha2s))
for k = 1:length(alpha2s)
    alpha2 = alpha2s[k]
    stepsize = 0.003/(1+alpha2/20)
    niter = 10000
    Ws, loss = twonet_square(X, Y, m, stepsize, niter, alpha1=alpha1, alpha2=alpha2, sym_init=true)
    WW[:,:,k] .= Ws[:,:,end]
end
#loglog(loss)

cm = plt.get_cmap("viridis")
figure(figsize=[8,4])
for k = 1:length(alpha2s)
X_test = cat(range(-1.2,1.2,length=200),ones(200),dims=2 )
    Y_test = (1/m) * sum( WW[:,end,k] .* max.( WW[:,1:end-1,k] * X_test', 0.0) , dims=1)[:]
plot(X_test[:,1], Y_test, color=cm((k-1)/(length(alpha2s)-1)),label="r=$(alpha2s[k]/50)",lw=3)
end
legend()
plot(X[:,1],Y,"ko",ms=10)
xticks([-1,0,1]);xlabel(L"x")
yticks([-1,0,1]);ylabel(L"y")

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
"""
    EdgeFtLayer(v_in=>v_out, e_in=>e_out)

Edge Features Layers. This is used in [this article] to compute a state representation of the TSPTW.

[this article]: https://arxiv.org/abs/2006.01610

# Arguments
- `v_in`: the dimension of input features for vertices.
- `v_out`: the dimension of output features for vertices.
- `e_in`: the dimension of input features for edges.
- `e_out`: the dimension of output features for edges.
"""
struct EdgeFtLayer{T <: Real}
    W_a::AbstractMatrix{T}
    W_T::AbstractMatrix{T}
    b_T::AbstractVector{T}
    W_e::AbstractMatrix{T}
    W_ee::AbstractMatrix{T}
    σ
end

function EdgeFtLayer(vDim::Pair{<:Integer,<:Integer}, eDim::Pair{<:Integer,<:Integer}, σ=identity; init=Flux.glorot_uniform, T::DataType=Float32, bias::Bool=true)

    vIn, vOut = vDim
    eIn, eOut = eDim

    # Used to compute node features
    W_a = T.(init(vOut, 2 * vIn + eIn))
    W_T = T.(init(vOut, 2 * vIn + eIn))
    b_T = bias ? T.(init(vOut)) : zeros(T, vOut)

    # Used to compute edge features
    W_e = T.(init(eOut, vIn))
    W_ee = T.(init(eOut, eIn))

    EdgeFtLayer(W_a, W_T, b_T, W_e, W_ee, σ)
end

Flux.@functor EdgeFtLayer

function (g::EdgeFtLayer)(fg::FeaturedGraph)
    mask = Flux.unsqueeze(fg.graph, 1) # used to broadcast arrays to the right shape
    infMask = - Inf32 .* (mask .< 0.5) # used to zero missing edges in the softmax

    srcNodes = Flux.unsqueeze(fg.nf, 2) # Fx1xN
    dstNodes = Flux.unsqueeze(fg.nf, 3) # FxNx1
    inpFeatures = vcat(srcNodes .* mask, fg.ef .* mask, dstNodes .* mask) # (F+E+F)xNxN

    attention = Flux.softmax(g.σ.(g.W_a ⊠ inpFeatures) .+ infMask, dims=2) # F'xNxN
    nodeFeatures = (sum(attention .* (g.W_T ⊠ inpFeatures), dims=2) .+ g.b_T)[:,1,:] # F'xN
    edgeFeatures = (g.W_e ⊠ (srcNodes .+ dstNodes) + g.W_ee ⊠ fg.ef) .* mask # E'xNxN

    return FeaturedGraph(
        fg.graph,
        nodeFeatures,
        edgeFeatures,
        fg.gf,
        fg.directed
    )
end

function (g::EdgeFtLayer)(fgs::BatchedFeaturedGraph)
    mask = Flux.unsqueeze(fgs.graph, 1) # used to broadcast arrays to the right shape
    infMask = - Inf32 .* (mask .< 0.5) # used to zero missing edges in the softmax

    srcNodes = Flux.unsqueeze(fgs.nf, 2) # Fx1xNxB
    dstNodes = Flux.unsqueeze(fgs.nf, 3) # FxNx1xB
    inpFeatures = vcat(srcNodes .* mask, fgs.ef .* mask, dstNodes .* mask) # (F+E+F)xNxNxB

    nodeTargetSize = (:, size(inpFeatures)[3:end]...)
    edgeTargetSize = (:, size(inpFeatures)[2:end]...)
    inpFeaturesFlatten = flatten_batch(inpFeatures) # (F+E+F)xNx(N*B)
    infMaskFlatten = flatten_batch(infMask) # 1xNx(N*B)
    nodeSumFlatten = flatten_batch(srcNodes .+ dstNodes) # FxNx(N*B)
    edgeFeaturesFlatten = flatten_batch(fgs.ef) # ExNx(N*B)

    attention = Flux.softmax(g.σ.(g.W_a ⊠ inpFeaturesFlatten) .+ infMaskFlatten, dims=2) # F'xNx(N*B)
    nodeFeatures = reshape(sum(attention .* (g.W_T ⊠ inpFeaturesFlatten), dims=2), nodeTargetSize) .+ g.b_T # F'xNxB
    edgeFeatures = reshape(g.W_e ⊠ nodeSumFlatten + g.W_ee ⊠ edgeFeaturesFlatten, edgeTargetSize) #E'xNxNxB

    return BatchedFeaturedGraph(
        fgs.graph,
        nodeFeatures,
        edgeFeatures,
        fgs.gf
    )
end


text:
function LinearGaussian(μ::Float64, σ::Float64)
    α = 5.0
    x = rand(:x, Normal(μ, σ))
    y = rand(:y, Normal(α * x, 1.0))
    z = rand(:z, Normal(y, 5.0))
    return z
end

function LinearGaussianProposal()
    α = 10.0
    x = rand(:x, Normal(α * 3.0, 3.0))
    y = rand(:y, Normal(0.0, 1.0))
end

function OneSiteProposal()
    x = rand(:x, Normal(0.0, 1.0))
end

@testset "Importance sampling" begin
    z = 3.0
    observations = Jaynes.target([(:z,) =>  z])
    n_traces = 5

    @testset "Linear Gaussian model" begin
        ps, nw = Jaynes.importance_sampling(observations, n_traces, LinearGaussian, (0.0, 1.0))
        @test length(ps.calls) == n_traces
        @test length(ps.lws) == n_traces
        @test isapprox(sum(nw), 1., atol = 1e-9)
        @test !isnan(ps.lmle)
        for call in ps.calls
            @test (call[:z]) == z
        end
    end

    @testset "Linear Gaussian proposal" begin
        ps, nw = Jaynes.importance_sampling(observations, n_traces, LinearGaussian, (0.0, 1.0), LinearGaussianProposal, ())
        @test length(ps.calls) == n_traces
        @test length(ps.lws) == n_traces
        @test isapprox(sum(nw), 1., atol = 1e-9)
        @test !isnan(ps.lmle)
        for call in ps.calls
            @test (call[:z]) == z
        end
    end
end


text:
include("../src/solvers/descent.jl")
# include("../src/solvers/faye_timing.jl")
# include("../src/solvers/dynprog_timing.jl")
include("../src/solvers/earliest_timing.jl")
include("../src/solvers/lp_timing.jl")

using Random: randperm
using StatsBase: sample, ProbabilityWeights
include("../src/solvers/neighbour_operators.jl")

# ===========
# Préparartion des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:outdir, "/tmp") # eqivalent à  -d/tmp
Args.set(:itermax, 300)

# ===========
inst = instance_build_mini10()
@test inst.nb_planes == 10

# ===========
lg1("1. Création sol avec args :lp... ")
# Args.set(:timing_algo_solver, :earliest)
Args.set(:timing_algo_solver, :lp)

Args.set(:level, 0)
sol = Solution(inst)
Args.set(:level, TEST_LEVEL)

@test Args.get("timing_algo_solver") == :lp
@test sol.timing_algo_solver == :lp
@test isa(sol.solver, LpTimingSolver)

@test sol.cost == 25650.0 # pour LpTimingSolver
# @test sol.cost == 25910.0 # pour EarliestTimingSolver

ln1("fait.")

# ===========
lg1("2. Création sol avec args :lp mais construite avec :earliest ... ")
# Args.set(:timing_algo_solver, :earliest)
Args.set(:timing_algo_solver, :lp)

sol = Solution(inst, algo=:earliest)

@test Args.get("timing_algo_solver") == :lp
@test sol.timing_algo_solver == :earliest
@test isa(sol.solver, EarliestTimingSolver)

# @test sol.cost == 25650.0 # pour LpTimingSolver
@test sol.cost == 25910.0 # pour EarliestTimingSolver

ln1("fait.")

# ===========
lg1("3. Création sol avec args :earliest... ")
Args.set(:timing_algo_solver, :earliest)

sol = Solution(inst)

@test Args.get("timing_algo_solver") == :earliest
@test sol.timing_algo_solver == :earliest
@test isa(sol.solver, EarliestTimingSolver)

# @test sol.cost == 25650.0 # pour LpTimingSolver
@test sol.cost == 25910.0 # pour EarliestTimingSolver

ln1("fait.")

# ===========
lg1("4. Création d'une sol alp01 aléatoire avec algo :earliest... ")
sol = Solution(inst, algo=:earliest)
# initial_sort!(solver.cursol, presort=:shuffle)
initial_sort!(sol, presort=:shuffle)
# @test sol.solver == nothing
@test sol.timing_algo_solver == :earliest
@test isa(sol.solver, EarliestTimingSolver)
# copy!(solver.bestsol, solver.cursol)
ln1("(sol.cost=$(sol.cost)) fait.")

# ===========
lg1("5. Création DescentSolver nb_cons_reject_max=$(Args.get(:itermax)) ")
ln1("(cost=700<=1000.0 ?)... ")
sv = DescentSolver(inst, startsol=sol)
sv.do_save_bestsol = false
sv.nb_cons_reject_max = 1000
sv.durationmax = 1
# Args.set(:level, 3)

Args.set(:level, 0)
solve(sv, permutation_operator!, startsol=sol, nb_cons_reject_max=300, durationmax=3)
Args.set(:level, TEST_LEVEL)

@test sv.durationmax == 3
@test sv.nb_cons_reject_max == 300
# @show typeof(sv.bestsol.solver)
# @show sv.bestsol.solver
@test typeof(sv.bestsol.solver) == EarliestTimingSolver
@test sv.bestsol.cost <= 2500.0 # 2110
@test sv.bestsol.cost >= 2110.0
ln1("   => fait : trouvé $(sv.bestsol.cost) (au lieu de 2110*))")
ln1("      ", to_s(sv.bestsol), "\n")


text:
function main()
	N = 2_000_000
	res = 0
	for i = 1:N
		if isprime(i)
			res += i
		end
	end
	println(res)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
using Test
using TestSetExtensions
using Logging
using Dates
using DataFrames
using DataStructures
import InfrastructureSystems
import InfrastructureSystems: Deterministic, Probabilistic, Scenarios, Forecast
using PowerSystems
using PowerGraphics
using PlotlyJS
using PowerSimulations
using GLPK
using Weave
using TimeSeries
using DataStructures

const PG = PowerGraphics
const IS = InfrastructureSystems
const PSY = PowerSystems
const PSI = PowerSimulations
const LOG_FILE = "PowerGraphics-test.log"

const BASE_DIR = dirname(dirname(pathof(PowerGraphics)))
const TEST_DIR = joinpath(BASE_DIR, "test")
template_dir = joinpath(BASE_DIR, "report_templates")
const generic_template = joinpath(template_dir, "generic_report_template.jmd")

LOG_LEVELS = Dict(
    "Debug" => Logging.Debug,
    "Info" => Logging.Info,
    "Warn" => Logging.Warn,
    "Error" => Logging.Error,
)

macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = TEST_DIR
        if length(tests) == 0
            tests = readdir(rootfile)
            tests = filter(
                f ->
                    startswith(f, "test_") && endswith(f, ".jl") && f != basename(rootfile),
                tests,
            )
        else
            tests = map(f -> string(f, ".jl"), tests)
        end
        println()
        for test in tests
            print(splitext(test)[1], ": ")
            include(joinpath(TEST_DIR, test))
            println()
        end
    end
end

function get_logging_level(env_name::String, default)
    level = get(ENV, env_name, default)
    log_level = get(LOG_LEVELS, level, nothing)
    if log_level === nothing
        error("Invalid log level $level: Supported levels: $(values(LOG_LEVELS))")
    end

    return log_level
end

function run_tests()
    include(joinpath(BASE_DIR, "test", "test_data", "get_test_data.jl"))
    console_level = get_logging_level("PS_CONSOLE_LOG_LEVEL", "Error")
    console_logger = ConsoleLogger(stderr, console_level)
    file_level = get_logging_level("PS_LOG_LEVEL", "Info")

    IS.open_file_logger(LOG_FILE, file_level) do file_logger
        levels = (Logging.Info, Logging.Warn, Logging.Error)
        multi_logger =
            IS.MultiLogger([console_logger, file_logger], IS.LogEventTracker(levels))
        global_logger(multi_logger)

        # Testing Topological components of the schema
        @time @testset "Begin PowerGraphics tests" begin
            @includetests ARGS
        end

        # TODO: once all known error logs are fixed, add this test:
        #@test length(get_log_events(multi_logger.tracker, Logging.Error)) == 0

        @info IS.report_log_summary(multi_logger)
    end
end

logger = global_logger()

try
    run_tests()
finally
    # Guarantee that the global logger is reset.
    @info("removing test files")
    #rm(TEST_OUTPUTS, recursive = true)
    global_logger(logger)
    nothing
end


text:
"""
    Node

struct
"""
mutable struct Node{K,D}
    parent::Union{Node{K,D},Nothing}
    left::Union{Node{K,D},Nothing}
    right::Union{Node{K,D},Nothing}
    key::K
    bf::Int8
    data::D
end # Node

Node{K,D}(key, data, parent) where {K,D} =
    Node{K,D}(parent, nothing, nothing, key, Int8(0), data)
Node{K,D}(key, data) where {K,D} = Node{K,D}(key, data, nothing)

Node(key::K, data::D) where {K,D} = Node{K,D}(key, data)
Node(key::K, data::D, parent::Union{Node{K,D},Nothing}) where {K,D} =
    Node{K,D}(key, data, parent)

Base.setproperty!(x::Node{K,D}, f::Symbol, v) where {K,D} =
    setfield!(x, f, v)

Base.show(io::IO, ::MIME"text/plain", node::Node{K,D}) where {K,D} =
    print(io, "Node{$(K),$(D)}: $(node.key) -> $(node.data)")


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"supertype for clocks in `Simulate.jl`"
abstract type AbstractClock end

"supertype for events"
abstract type AbstractEvent end

"supertype for clock states"
abstract type ClockState end

"supertype for clock events"
abstract type ClockEvent end

"""
    Timing

Enumeration type for scheduling events and timed conditions:

- `at`: schedule an event at a given time,
- `after`: schedule an event a given time after current time,
- `every`: schedule an event every given time from now on,
- `before`: a timed condition is true before a given time,
- `until`: delay until t.
"""
@enum Timing at after every before until

"""
```
Fun(f::Function, arg...; kw...)
```
Store a function and its arguments for being called later as an event.

# Arguments, fields
- `f::Function`:  event function to be executed at event time,
- `arg...`: arguments to the event function,
- `kw...`: keyword arguments to the event function.

Arguments and keyword arguments can be 1) values or variables mixed with 2) symbols,
expressions or even other functions. In the 2nd cases they are evaluated at
event time before they are passed to the event function.

!!! note
    Composite types or variables given symbolically can change until they
    are evaluated later at event time.

# Examples
```jldoctest
julia> using Simulate

julia> f(a,b,c; d=4, e=5) = a+b+c+d+e       # if you define a function and ...
f (generic function with 1 method)

julia> sf = Fun(f, 10, 20, 30, d=14, e=15);  # store it as Fun

julia> sf.f(sf.arg...; sf.kw...)         # it can be executed later
89

julia> d = Dict(:a => 1, :b => 2);          # we set up a dictionary

julia> g(t) = t[:a] + t[:b]                 # and a function adding :a and :b
g (generic function with 1 method)

julia> g(d)                                 # our add function gives 3
3

julia> ff = Fun(g, d);              # we set up a Fun

julia> d[:a] = 10;                          # later somehow we change d

julia> ff.f(ff.arg...)                   # calling ff then gives a different result
12
```
"""
struct Fun
    f::Function
    arg::Union{Nothing, Tuple}
    kw::Union{Nothing, Base.Iterators.Pairs}

    Fun(f::Function, arg...; kw...) =
        new(f, ifelse(isempty(arg), nothing, arg), ifelse(isempty(kw), nothing, kw))
end

"An action is either an `Expr` or a `Fun` or a `Tuple` of them."
const Action = Union{Expr, Fun, Tuple}

"""
```
DiscreteEvent{T<:Action}
```
A discrete event is a function or an expression or a tuple of them to be
executed at an event time.

# Arguments, fields
- `ex::T`: a function or an expression or a tuple of them,
- `scope::Module`: evaluation scope,
- `t::Float64`: event time,
- `Δt::Float64`: repeat rate with for repeating events.
"""
struct DiscreteEvent{T<:Action} <: AbstractEvent
    ex::T
    scope::Module
    t::Float64
    Δt::Float64
end
DiscreteEvent(ex::T, scope::Module, t::Real, Δt::Real) where {T<:Action} =
    DiscreteEvent(ex, scope, float(t), float(Δt))

"""
```
DiscreteCond{S<:Action, T<:Action}
```
A condition to be evaluated repeatedly with expressions or functions
to be executed if conditions are met.

# Arguments, fields
- `cond::S`: a conditional function or an expression or a tuple of them
    (conditions must evaluate to `Bool`),
- `ex::T`: a function or an expression or a tuple of them to be executed
    if conditions are met,
- `scope::Module`: evaluation scope
"""
struct DiscreteCond{S<:Action, T<:Action} <: AbstractEvent
    cond::S
    ex::T
    scope::Module
end

"""
    Sample{T<:Union{Fun,Expr}}

A sampling function or expression is called at sampling time.

# Arguments, fields
- `ex::T`: expression or function to be called at sample time,
- `scope::Module`: evaluation scope.
"""
struct Sample{T<:Union{Fun,Expr}} <: AbstractEvent
    ex::T
    scope::Module
end

"""
    ClockException(ev::ClockEvent, value=nothing)

Define a ClockException, which can be thrown to processes.

# Arguments, fields
- `ev::ClockEvent`: delivers an event to the interrupted task
- `value=nothing`: deliver some other value
"""
struct ClockException <: Exception
  ev::ClockEvent
  value::Any
  ClockException(ev::ClockEvent, value=nothing) = new(ev, value)
end

"""
```
Prc( id, f::Function, arg...; kw...)
alias   SP( id, f::Function, arg...; kw...)
```
Prepare a function to run as a process in a simulation.

# Arguments, fields
- `id`: some unique identification for registration,
- `task::Union{Task,Nothing}`: a task structure,
- `clk::Union{AbstractClock,Nothing}`: clock where the process is registered,
- `state::ClockState`: process state,
- `f::Function`: a function `f(clk, arg...; kw...)`, must take `clk` as its
    first argument,
- `arg...`: further arguments to `f`
- `kw...`: keyword arguments to `f`

!!! note
    A function started as a Prc most often runs in a loop. It has to
    give back control by e.g. doing a `take!(input)` or by calling
    `delay!` etc., which will `yield` it. Otherwise it will starve everything else!

!!! warn
    That `f` nust take `clk` as first argument is a breaking change in v0.3!

# Examples
```jldoctest
julia> using Simulate

```
"""
mutable struct Prc
    id::Any
    task::Union{Task,Nothing}
    clk::Union{AbstractClock,Nothing}
    f::Function
    arg::Tuple
    kw::Base.Iterators.Pairs

    Prc( id, f::Function, arg...; kw...) = new(id, nothing, nothing, f, arg, kw)
end

"""
    Schedule()

A Schedule contains events, conditional events and sampling functions to be
executed or evaluated on the clock's time line.

# Fields
- `events::PriorityQueue{DiscreteEvent,Float64}`: scheduled events,
- `cevents::Array{DiscreteCond,1}`: conditional events to evaluate at each tick,
- `samples::Array{Sample,1}`: sampling expressions to evaluate at each tick,
"""
mutable struct Schedule
    events::PriorityQueue{DiscreteEvent,Float64}
    cevents::Array{DiscreteCond,1}
    samples::Array{Sample,1}

    Schedule() = new(PriorityQueue{DiscreteEvent,Float64}(), DiscreteCond[], Sample[])
end

"""
    ClockChannel

Provide a channel to an active clock or a real time clock.

# Fields
- `ref::Ref{Task}`: a pointer to an active clock,
- `ch::Channel`: a communication channel to an active clock,
- `id::Int`: the thread id of the active clock.
- `done::Bool`: flag indicating if the active clock has completed its cycle.
"""
mutable struct ClockChannel{T <: ClockEvent}
    ref::Ref{Task}
    forth::Channel{T}
    back::Channel{T}
    thread::Int
    done::Bool
end

"""
```
Clock(Δt::Number=0; t0::Number=0, unit::FreeUnits=NoUnits)
```
Create a new simulation clock.

# Arguments
- `Δt::Number=0`: time increment. If no Δt is given, the simulation doesn't tick,
    but jumps from event to event. Δt can be set later with `sample_time!`.
- `t0::Number=0`: start time for simulation
- `unit::FreeUnits=NoUnits`: clock time unit. Units can be set explicitely by
    setting e.g. `unit=minute` or implicitly by giving Δt as a time or else setting
    t0 to a time, e.g. `t0=60s`.

# Fields
- `id::Int`: thread on which the clock is running,
- `state::ClockState`: clock state,
- `time::Float64`: clock time,
- `unit::FreeUnits`: time unit,
- `end_time::Float64`: end time for simulation,
- `Δt::Float64`: sampling time, timestep between ticks,
- `ac::Vector{ClockChannel}`: active clocks running on parallel threads,
- `sc::Schedule`: the clock schedule (events, cond events and sampling),
- `processes::Dict{Any, Prc}`: registered `Prc`es,
- `tn::Float64`: next timestep,
- `tev::Float64`: next event time,
- `evcount::Int`: event counter,
- `scount::Int`: sample counter

# Examples

```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s, minute, hr

julia> c = Clock()                 # create a unitless clock (standard)
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c = Clock(1s, unit=minute)  # create a clock with units, does conversions automatically
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=0.0 minute, Δt=0.01667 minute, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c = Clock(1s)               # create a clock with implicit unit setting
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=0.0 s, Δt=1.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c = Clock(t0=60s)           # another example of implicit unit setting
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=60.0 s, Δt=0.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c = Clock(1s, t0=1hr)       # if given times with different units, Δt takes precedence
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=3600.0 s, Δt=1.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
mutable struct Clock <: AbstractClock
    id::Int
    state::ClockState
    time::Float64
    unit::FreeUnits
    Δt::Float64
    ac::Vector{ClockChannel}
    sc::Schedule
    processes::Dict{Any, Prc}
    tn::Float64
    tev::Float64
    end_time::Float64
    evcount::Int
    scount::Int

    function Clock(Δt::Number=0;
                   t0::Number=0, unit::FreeUnits=NoUnits)
        if 1unit isa Time
            Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
        elseif Δt isa Time
            unit = Unitful.unit(Δt)
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
            Δt = Δt.val
        elseif t0 isa Time
            unit = Unitful.unit(t0)
            t0 = t0.val
        else
            nothing
        end
        new(0, Undefined(), t0, unit, Δt, ClockChannel[], Schedule(),
            Dict{Any, Prc}(), t0 + Δt, t0, t0, 0, 0)
    end
end

"""
```
ActiveClock(clock::Clock, master::Ref{Clock},
            cmd::Channel{ClockEvent}, ans::Channel{ClockEvent})
```
A thread specific clock which can be operated via a channel.

# Fields
- `clock::Clock`: the thread specific clock,
- `master::Ref{Clock}`: a pointer to the master clock (on thread 1),
- `cmd::Channel{ClockEvent}`: the command channel from master,
- `ans::Channel{ClockEvent}`: the response channel to master.
- `id::Int`: the id in master's ac array,
- `thread::Int`: the thread, the active clock runs on.

!!! note
You should not setup an `ActiveClock` explicitly. Rather this is done
implicitly by [`fork!`](@ref)ing a [`Clock`](@ref) to other available threads
or directly with [`PClock`](@ref).
It then can be accessed via [`pclock`](@ref) as in the following example.

# Example
```jldoctest
julia> using Simulate

julia> clk = Clock()
Clock thread 1 (+ 0 ac): state=Simulate.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> fork!(clk)

julia> clk    #  ⬇ you got 3 parallel active clocks
Clock thread 1 (+ 3 ac): state=Simulate.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> clk = PClock()
Clock thread 1 (+ 3 ac): state=Simulate.Undefined(), t=0.0 , Δt=0.01 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> pclock(clk, 1)    # get access to the 1st active clock (on thread 2)
Active clock 1 on thrd 2: state=Simulate.Idle(), t=0.0 , Δt=0.01 , prc:0
   scheduled ev:0, cev:0, sampl:0
```
"""
mutable struct ActiveClock{T <: ClockEvent} <: AbstractClock
    clock::Clock
    master::Ref{Clock}
    forth::Channel{T}
    back ::Channel{T}
    id::Int
    thread::Int
end

"""
```
RTClock
```
Real time clocks use system time for time keeping and are controlled over
a channel. They are independent from each other and other clocks and run
asynchronously as tasks on arbitrary threads. Multiple real time clocks
can be setup with arbitrary frequencies to do different things.

# Fields
- `clock::Clock`:
- `cmd::Channel{T}`:
- `back::Channel{T}`:
- `id::Int`:
- `thread::Int`:
- `t0::Float64`:
"""
mutable struct RTClock{T <: ClockEvent} <: AbstractClock
    clock::Clock
    cmd::Channel{T}
    back::Channel{T}
    id::Int
    thread::Int
    t0::Float64
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractLearningExample

A learning example (e.g. points, tiles with labels).
"""
abstract type AbstractLearningExample end

"""
    Base.iterate(learnexample, state)

Generate a learning example `learnexample` from a given `state`.
"""
Base.iterate(::AbstractLearningExample, state=1) = @error "not implemented"

"""
    AbstractLabeledExample

A learning example for supervised learning tasks.
"""
abstract type AbstractLabeledExample <: AbstractLearningExample end

"""
    AbstractUnlabeledExample

A learning example for unsupervised learning tasks.
"""
abstract type AbstractUnlabeledExample <: AbstractLearningExample end

#------------------
# IMPLEMENTATIONS
#------------------
include("examples/point_example.jl")


text:
@testset "Utils - containedIn" begin
	@test containedin([5], [5])
	@test containedin([1,2,3], [1,2,3])
	@test containedin([1,2,3], [1,5,2,5,3,5])
end


text:
_isapprox(x::Real, y::Real) = _isapprox(promote(x, y)...)
_isapprox(x::T, y::T) where {T<:Real} = x == y
_isapprox(x::T, y::T) where {T<:AbstractFloat} = y < x+1024*eps(T) && x < y+1024*eps(T)
_isapprox(x::AbstractVector{S}, y::AbstractVector{T}) where {S<:Real,T<:Real} = _isapprox(promote(x, y)...)
_isapprox(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Real} = x == y
_isapprox(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat} = _isapproxzero(norm(x - y))
_isapproxzero(x::T) where {T<:Real} = _isapprox(x, zero(T))

function _parallel(x, y)
    @assert length(x) == length(y)
    cstchosen = false
    cst = zero(Base.promote_op(/, eltype(y), eltype(x)))
    for i in 1:length(x)
        if _isapproxzero(x[i])
            _isapproxzero(y[i]) || return false
        elseif _isapproxzero(y[i])
            _isapproxzero(x[i]) || return false
        else
            c = y[i] / x[i]
            if cstchosen
                _isapprox(cst, c) || return false
            else
                cst = c
            end
        end
    end
    return true
end

function inaffspace(r, s, L, par::Bool=true)
    x = Polyhedra.coord(r)
    y = Polyhedra.coord(s)
    for l in L
        z = Polyhedra.coord(l)
        # remove component
        x = x * dot(z, z) - z * dot(z, x)
        y = y * dot(z, z) - z * dot(z, y)
    end
    par ? _parallel(x, y) : _isapprox(r, s)
end

function inequality_fulltest(hr::HRepresentation, exp::HRepresentation)
    @test nhyperplanes(hr) == nhyperplanes(exp)
    @test nhalfspaces(hr) == nhalfspaces(exp)

    _lift(h::Polyhedra.HRepElement) = [h.β; h.a]
    aff = map(_lift, hyperplanes(hr))
    inaff(x, y) = inaffspace(x, y, aff)

    for h in hyperplanes(exp)
        hl = _lift(h)
        @test any(hp -> inaff(hl, _lift(hp)), hyperplanes(hr))
    end
    for h in halfspaces(exp)
        hl = _lift(h)
        @test any(hs -> inaff(hl, _lift(hs)), halfspaces(hr))
    end
end
function inequality_fulltest(h::HRepresentation, hrepargs...)
    inequality_fulltest(h, hrep(hrepargs...))
end
function inequality_fulltest(p::Polyhedron, hrepargs...)
    detecthlinearity!(p)
    removehredundancy!(p)
    inequality_fulltest(hrep(p), hrepargs...)
end

function generator_fulltest(vr::VRepresentation, exp::VRepresentation)
    @test npoints(vr) == npoints(exp)
    @test nlines(vr) == nlines(exp)
    @test nrays(vr) == nrays(exp)
    linspace = collect(lines(vr))
    inaff(x, y, args...) = inaffspace(x, y, linspace, args...)
    for l in lines(exp)
        @test any(ll -> inaff(l, ll), lines(vr))
    end
    for r in rays(exp)
        @test any(s -> inaff(r, s), rays(vr))
    end
    for p in points(exp)
        @test any(q -> inaff(p, q, false), points(vr))
    end
end
function generator_fulltest(v::VRepresentation, vrepargs...)
    generator_fulltest(v, vrep(vrepargs...))
end
function generator_fulltest(p::Polyhedron, vrepargs...)
    detectvlinearity!(p)
    removevredundancy!(p)
    generator_fulltest(vrep(p), vrepargs...)
end


text:
using Base.Iterators: Iterators, product
using ..Datasets: getsavefolder, Cube

"""
    savecube(cube,name::String)

Save a [`YAXArray`](@ref) to the folder `name` in the YAXArray working directory.
"""
function savecube(c, name::AbstractString;
  chunksize = Dict(),
  max_cache = 1e8,
  backend = :zarr,
  backendargs...
)
allax = caxes(c)
if !isa(chunksize,AbstractDict) && !isa(chunksize,NamedTuple)
  @warn "Chunksize must be provided as a Dict mapping axis names to chunk size in the future"
  chunksize = OrderedDict(axname(i[1])=>i[2] for i in zip(allax,chunksize))
end
  firstaxes = sort!([findAxis(String(k),allax) for k in keys(chunksize)])
  lastaxes = setdiff(1:length(allax),firstaxes)
  allax = allax[[firstaxes;lastaxes]]
  dl = cumprod(length.(caxes(c)))
  isplit = findfirst(i->i>max_cache/sizeof(eltype(c)),dl)
  isplit isa Nothing && (isplit=length(dl)+1)
  forcesingle = (isplit+1)<length(firstaxes)
  axin = allax[1:isplit-1]
  axn = axname.(axin)
  indims = InDims(axn...)
  outdims = OutDims(axn..., backend=backend,chunksize=chunksize, path = name; backendargs...)
  function cop(xout,xin)
    xout .= xin
  end
  if forcesingle
    nprocs()>1 && println("Forcing single core processing because of bad chunk size")
    o = mapCube(cop,c,indims=indims, outdims=outdims,ispar=false,max_cache=max_cache)
  else
    o = mapCube(cop,c,indims=indims, outdims=outdims,max_cache=max_cache)
  end
end

function loadcube(s)
  Cube(getsavefolder(s, true))
end

function rmcube(s)
  p = getsavefolder(s, true)
  if isfile(p) || isdir(p)
    rm(p, recursive=true)
  end
end

function writefun(xout,xin::AbstractArray{Union{Missing,T}},a,nd,cont_loop,filename) where T

  #TODO replace -9999.0 this with a default missing value
  x = map(ix->ismissing(ix) ? convert(T,-9999.0) : ix,xin)

  count_vec = fill(-1,nd)
  start_vec = fill(1,nd)

  used_syms = Symbol[]
  for (k,v) in cont_loop
    count_vec[k] = 1
    start_vec[k] = a[Symbol(v)][1]
    push!(used_syms,Symbol(v))
  end

  splinds = Iterators.filter(i->!in(i,used_syms),keys(a))
  vn = join(string.([a[s][2] for s in splinds]),"_")
  isempty(vn) && (vn="layer")
  ncwrite(x,filename,vn,start=start_vec,count=count_vec)
end

global const projection = Dict(
"grid_mapping_name" => "transverse_mercator",
"longitude_of_central_meridian" => -9. ,
"false_easting" => 500000. ,
"false_northing" => 0. ,
"latitude_of_projection_origin" => 0. ,
"scale_factor_at_central_meridian" => 0.9996,
"long_name" => "CRS definition",
"longitude_of_prime_meridian" => 0.,
"semi_major_axis" => 6378137.,
"inverse_flattening" => 298.257223563 ,
"spatial_ref" => "PROJCS[\"WGS 84 / UTM zone 29N\",GEOGCS[\"WGS 84\",
DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,
AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],
PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],
UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],
AUTHORITY[\"EPSG\",\"4326\"]],
PROJECTION[\"Transverse_Mercator\"],
PARAMETER[\"latitude_of_origin\",0],
PARAMETER[\"central_meridian\",-9],
PARAMETER[\"scale_factor\",0.9996],
PARAMETER[\"false_easting\",500000],
PARAMETER[\"false_northing\",0],
UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],
AXIS[\"Easting\",EAST],
AXIS[\"Northing\",NORTH],
AUTHORITY[\"EPSG\",\"32629\"]]",
"GeoTransform" => "722576.2320803495 20 0 4115483.464603715 0 -20 ")

global const epsg4326=Dict(
  "grid_mapping_name" => "latitude_longitude",
  "long_name" => "CRS definition",
  "longitude_of_prime_meridian" => 0.,
  "semi_major_axis" => 6378137.,
  "inverse_flattening" => 298.257223563,
  "spatial_ref" => "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
  "GeoTransform" => "-98.68712696894481 0.0001796630568239077 0 20.69179551612753 0 -0.0001796630568239077 "
)


text:
"""
	$(SIGNATURES)

Given a range `1:n`, return indices in "interval-halving" order.
The purpose is to determine the order in which ordered computations should be performed efficiently when surrounding values provide bounds for inner values.

# Example
```
bisecting_indices(5) == [1, 5, 3, 2, 4]
```
"""
function bisecting_indices(iLow :: I1, iHigh :: I1) where I1 <: Integer
    n = iHigh - iLow + 1;
    if n < 3
        return iLow : iHigh;
    end

    idxV = zeros(I1, n);
    idxV[1] = iLow;
    idxV[2] = iHigh;
    iMid = round(I1, (iLow + iHigh) / 2);
    idxV[3] = iMid;

    # Lower half
    nLow = iMid - iLow - 1;
    if nLow >= 1
        idxV[4 : (3 + nLow)] = bisect_interior(iLow, iMid);
    end

    # Upper half
    nHigh = iHigh - iMid - 1;
    if nHigh >= 1
        idxV[(4 + nLow) : n] = bisect_interior(iMid, iHigh);
    end
    return idxV
end

function bisect_interior(iLow :: I1, iHigh :: I1) where I1 <: Integer
    @assert iHigh > iLow + one(I1)
    idxV = zeros(I1, iHigh - iLow - 1);

    iMid = round(I1, (iLow + iHigh) / 2);
    idxV[1] = iMid;

    # Lower half
    nLow = iMid - iLow - 1;
    if nLow > 0
        idxV[2 : (1 + nLow)] .= bisect_interior(iLow, iMid);
    else
        nLow = 0;
    end

    nHigh = iHigh - iMid - 1;
    if nHigh > 0
        idxHigh1 = 2 + nLow;
        idxV[idxHigh1 : (idxHigh1 - 1 + nHigh)] .= bisect_interior(iMid, iHigh);
    end

    return idxV
end

"""
	$(SIGNATURES)

For each value in `valueV` (a weakly increasing `Vector`), find the corresponding indices in `gridV` (also an increasing `Vector{Integer}`).
If no match is found, set indices to 0.
Returns a `Vector` even if only one value is given.

## Example
```
julia> find_indices([2, 4], -2 : 2 : 10) == [3, 4]
true
```
"""
function find_indices(valueV, gridV;  notFoundError :: Bool = false)
    @argcheck eltype(valueV) <: Integer
    nGrid = length(gridV);
    idxV = zeros(Int, size(valueV));
    # Last value of `gridV` that was matched
    lastIdx = 1;
    for (j, value) in enumerate(valueV)
        # lastIdx = find_index(value, gridV; startIdx = lastIdx);
        # if !isnothing(lastIdx)
        #     idxV[j] = lastIdx;
        # else
            # # No grid point matches current value. Done.
            # break
        # end
        lastIdx = findnext(x -> x >= value, gridV, lastIdx);
        if isnothing(lastIdx)
            # No grid point matches current value. Done.
            break
        end
        if value == gridV[lastIdx]
            idxV[j] = lastIdx;
        end
    end

    if notFoundError  &&  any(idxV .== 0)
        error("Some values not matched.")
    end
    return idxV
end

"""
	$(SIGNATURES)

Find the next grid point starting with index `startIdx` that is equal to `value`.

# Example
```
julia> find_index(6.0, collect(1.0 : 10.0); startIdx = 3)
6
```
"""
function find_index(value :: T1, gridV;  startIdx = 1)  where T1 <: Number
    # First grid point (weakly) above current value
    idx = findnext(x -> x >= value, gridV, startIdx);
    if isnothing(idx)
        return nothing
    elseif isequal(value, gridV[idx])
        return idx
    else
        return nothing
    end
end

## --------  Comparisons
# Easier to write versions that check whether all elements of a vector are inside bounds
# Must be efficient!

"""
	$(SIGNATURES)

Check that all elements are above a scalar lower bound.
"""
all_at_least(xV, lb :: Real) = all(x -> x >= lb, xV);

all_at_least(xV :: AbstractArray{T1,N}, lbV :: AbstractArray{T2,N}) where {T1,T2,N} =
    all(j -> xV[j] >= lbV[j], eachindex(xV));

"""
	$(SIGNATURES)

Check that all elements are below a scalar upper bounds.
"""
all_at_most(xV, ub :: Real) = all(x -> x <= ub, xV);

all_at_most(xV :: AbstractArray{T1,N}, lbV :: AbstractArray{T2,N}) where {T1,T2,N} =
    all(j -> xV[j] <= lbV[j], eachindex(xV));

"""
	$(SIGNATURES)

Strictly check that all elements are above a lower bound.

# Example
```
julia> all_greater([1.0, 2.0], 0.9999; atol = 0.01)
true
```
"""
all_greater(xV, lb :: F1; atol = zero(F1)) where F1 <: Real =
    all(x -> x > lb - atol, xV);

all_greater(xV :: AbstractArray{<:Real}, yV :: AbstractArray{T}; atol = zero(T)) where T <: Real =
    all(map((x,y) -> x > y - atol,  xV, yV));

"""
	$(SIGNATURES)

Strictly check that all elements are below an upper bound.
"""
all_less(xV, ub :: F1; atol = zero(F1)) where F1 <: Real =
    all(x -> x < ub + atol, xV);

all_less(xV :: AbstractArray{<:Real}, yV :: AbstractArray{T}; atol = zero(T)) where T <: Real =
    all(map((x,y) -> x < y + atol,  xV, yV));

"""
	$(SIGNATURES)

Check that at least some element is above a lower bound.
"""
any_at_least(xV, lb) = any(x -> x >= lb, xV);

"""
	$(SIGNATURES)

Check that at least some element is below an upper bound.
"""
any_at_most(xV, ub) = any(x -> x <= ub, xV);

"""
	$(SIGNATURES)

Strictly check that at least some element is above a lower bound.
"""
any_greater(xV, lb :: F1; atol = zero(F1)) where F1 =
    any(x -> x > lb + atol, xV);

"""
	$(SIGNATURES)

Check that at least some element is below an upper bound.
"""
any_less(xV, ub :: F1; atol = zero(F1)) where F1 =
    any(x -> x < ub - atol, xV);

"""
	$(SIGNATURES)

Check if any array element is NaN.
"""
any_nan(xV) = any(x -> isnan(x), xV);

# module VectorLH

# using StatsBase

# export countbool, count_elem, count_indices, counts_from_fractions, counts_to_indices
# export scale_vector!, validate_vector

# # Count number of true elements in a vector
# function countbool(xV :: Array{Bool, 1})
#   count(z -> (z == true), xV)
# end

# """
# count_elem()

# Count the number of times each element occurs in a vector
# NaN is kept as a separate value unless `omitNan` is set

# OUT
#     outM
#         matrix with values and counts as columns
# """
# function count_elem(xV :: Array{T,1}, omitNan = false) where T <: Any
#     if omitNan
#         dictV = countmap(filter(y -> ~isnan(y), xV));
#     else
#         dictV = countmap(xV);
#     end

#     # Make dict into vector of counts
#     countM = hcat([[key, val] for (key, val) in dictV]...)';
#     # countV = zeros(T, size(xV));
#     # for (key, val) in dictV
#     #     countV[key] = val;
#     # end
#     return sortslices(countM, dims=1)
# end

# """
# count_indices()

# Count the number of times each element occurs in a vector
# Return vector of counts
# Assumes that values are integer indices
# """
# function count_indices(xV :: Array{T,1}) where T <: Integer
#     countM = count_elem(xV, false);
#     countV = zeros(Int64, maximum(countM[:,1]));
#     countV[countM[:,1]] = countM[:, 2];
#     return countV;
# end

# """
# Counts from fractions

# Given a total count and the fraction of time each value should occur
# compute the number of times each value occurs

# Fractions must sum to 1
# """
# function counts_from_fractions(fracV :: Vector{T1}, n :: Integer;
#     dbg :: Bool = false) where
#     T1 <: AbstractFloat

#     cumNV = round.(Int, cumsum(fracV) .* n);
#     countV = diff(vcat(0, cumNV));

#     if dbg
#         @assert sum(countV) == n
#         @assert length(countV) == length(fracV)
#         @assert all(fracV .>= 0.0)
#         @assert sum(fracV) ≈ 1.0
#     end
#     return countV :: Vector{Int}
# end

# """
# Counts to indices

# Given a vector of counts, return indices in a "stacked" vector
# First entry is (1, countV[1])

# OUT
#     Vector with start / end indices for each entry
#     Note: Cannot use Tuples b/c no empty Tuple{Int,Int} exists.
#     Empty when countV = 0
# """
# function counts_to_indices(countV :: Vector{T1}; dbg :: Bool = false) where
#     T1 <: Integer

#     n = length(countV);
#     if dbg
#         @assert n > 0
#         @assert all(countV .>= 0)
#     end

#     outV = Vector{Vector{Int}}(undef, n);
#     idxLast = 0;
#     for i1 = 1 : n
#         if countV[i1] == 0
#             outV[i1] = Int[];
#         else
#             outV[i1] = [idxLast + 1,  idxLast + countV[i1]];
#             idxLast += countV[i1];
#         end
#     end

#     if dbg
#         @assert idxLast == sum(countV)
#         @assert length(outV) == n
#     end
#     return outV
# end # function

# """
# ## Scale a vector

# Impose minimum, maximum values. Ensure that it sums to a fixed value

# If the sum is fixed, this currently does not guarantee that the min / max
# are satisfied. One would have to call the function iteratively to ensure this.
# """
# function scale_vector!(v :: Vector{T1};  xMin = [],  xMax = [],  vSum = [],
#     dbg :: Bool = false)  where T1 <: AbstractFloat

#     if dbg
#         # Check that scaling is feasible
#         n = length(v);
#         if !isempty(xMax)  &&  !isempty(xMin)
#             @assert xMax > xMin
#         end
#         if !isempty(xMin)  &&  !isempty(vSum)
#             @assert n * xMin < vSum
#         end
#         if !isempty(xMax)  &&  !isempty(vSum)
#             @assert n * xMax > vSum
#         end
#     end

#     if !isempty(xMin)
#         idxMinV = findall(v .< xMin);
#         v[idxMinV] .= xMin;
#     else
#         idxMinV = [];
#     end
#     if !isempty(xMax)
#         idxMaxV = findall(v .> xMax);
#         v[idxMaxV] .= xMax;
#     else
#         idxMaxV = [];
#     end
#     if !isempty(vSum)
#         scaledSum = sum(v);
#         canChangeV = trues(size(v));
#         if scaledSum > vSum
#             # Cannot change the values that hit the min
#             canChangeV[idxMinV] .= false;
#         else
#             canChangeV[idxMaxV] .= false;
#         end
#         scaledSum = sum(v[canChangeV]);
#         v[canChangeV] .*= (vSum - sum(v[.!canChangeV])) / sum(v[canChangeV]);

#         if dbg
#             @assert sum(v) ≈ vSum
#         end
#     end
#     return nothing
# end

# function validate_vector(v :: Vector{T1};  xMin = [],  xMax = [],
#     vSum = [])  where T1 <: AbstractFloat

#     isValid = true;
#     if !isempty(xMin)  &&  any(v .< xMin)
#         isValid = false;
#         vMin = minimum(v);
#         @warn "Values below minimum:  $vMin < $xMin"
#     end
#     if !isempty(xMax)  &&  any(v .> xMax)
#         isValid = false;
#         @warn "Values above maximum"
#     end
#     if !isempty(vSum)  &&  !(sum(v) ≈ vSum)
#         isValid = false;
#         @warn "Sum is wrong"
#     end
#     return isValid
# end

# end # module


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

using QuadGK

function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
    if method == "right"
        meth = f -> (lr -> begin l,r = lr; f(r) * (r-l) end)
    elseif method == "left"
        meth = f -> (lr -> begin l,r = lr; f(l) * (r-l) end)
    elseif method == "trapezoid"
        meth = f -> (lr -> begin l,r = lr; (1/2) * (f(l) + f(r)) * (r-l) end)
    elseif method == "simpsons"
        meth = f -> (lr -> begin l,r=lr; (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l) end)
    end

    xs = a .+ (0:n) * (b-a)/n

    sum(meth(f), zip(xs[1:end-1], xs[2:end]))
end

"""
    integration(f)

Show graphically the left Riemann approximation, the trapezoid approximation, and Simpson's approximation to the integral of `f` over [-1,1].

Assumes `f` is non-negative.
"""
function integration(f=nothing)
    if f == nothing
        f = x -> x^2*exp(-x/3)
    end

    a, b = -1, 1

function left_riemann_pts(n)
    xs = range(a, b, length=n+1)
    pts = Point2f0[(xs[1], 0)]
    for i in 1:n
        xᵢ, xᵢ₊₁ = xs[i], xs[i+1]
        fᵢ = f(xᵢ)
        push!(pts, (xᵢ, fᵢ))
        push!(pts, (xᵢ₊₁, fᵢ))
    end
    push!(pts, (xs[end], 0))
    pts
end

function trapezoid_pts(n)
    xs = range(a, b, length=n+1)
    pts = Point2f0[(xs[1], 0)]
    for i in 1:n
        xᵢ, xᵢ₊₁ = xs[i], xs[i+1]
        fᵢ = f(xᵢ)
        push!(pts, (xᵢ, f(xᵢ)))
    end
    push!(pts, (xs[end], f(xs[end])))
    push!(pts, (xs[end], 0))
    pts
end

function simpsons_pts(n)
    xs = range(a, b, length=n+1)
    pts = Point2f0[(xs[1], 0), (xs[1], f(xs[1]))]

    for i in 1:n
        xi, xi1 = xs[i], xs[i+1]
        m = xi/2 + xi1/2
        p = x -> f(xi)*(x-m)*(x - xi1)/(xi-m)/(xi-xi1) + f(m) * (x-xi)*(x-xi1)/(m-xi)/(m-xi1) + f(xi1) * (x-xi)*(x-m) / (xi1-xi) / (xi1-m)

        xs′ = range(xi, xi1, length=10)
        for j in 2:10
            x = xs′[j]
            push!(pts, (x, p(x)))
        end
    end
    push!(pts, (xs[end], 0))
    pts
end

# where we lay our scene:
scene, layout = layoutscene()
layout.halign = :left
layout.valign = :top

p1 = layout[1,1] = LScene(scene)
p2 = layout[1,2] = LScene(scene)
p3 = layout[1,3] = LScene(scene)

n = layout[2,1] = LSlider(scene, range=2:25, startvalue=2)
output = layout[2, 2:3] = LText(scene, "...")

lpts = lift(n.value) do n
    left_riemann_pts(n)
end

poly!(p1.scene, lpts, color=:gray75)
lines!(p1.scene, a..b, f, color=:black, strokewidth=10)
title(p1.scene, "Left Riemann")

tpts = lift(n.value) do n
    trapezoid_pts(n)
end
poly!(p2.scene, tpts, color=:gray75)
lines!(p2.scene, a..b, f, color=:black, strokewidth=10)
title(p2.scene, "Trapezoid")

spts = lift(n.value) do n
    simpsons_pts(n)
end
poly!(p3.scene, spts, color=:gray75)
lines!(p3.scene, a..b, f, color=:black, strokewidth=10)
title(p3.scene, "Simpson's")

on(n.value) do n
    actual,err = quadgk(f, -1, 1)
    lrr  = riemann(f, -1, 1, n, method="left")
    trap = riemann(f, -1, 1, n, method="trapezoid")
    simp = riemann(f, -1, 1, n, method="simpsons")

    Δleft = round(abs(lrr - actual), digits=8)
    Δtrap = round(abs(trap - actual), digits=8)
    Δsimp = round(abs(simp - actual), digits=8)

    txt = "Riemann: $Δleft, Trapezoid: $Δtrap, Simpson's: $Δsimp"
    output.text[] = txt
end
n.value[] = n.value[]

scene
end


text:
include("types/__init__.jl")
include("printing/__init__.jl")

include("utils/__init__.jl")
include("corpus/__init__.jl")

include("matchers/__init__.jl")
include("proving/__init__.jl")


text:
"""
This is a library of functions to help out in the Robotlab at LTH
The module imports the following files\n
include("utils.jl")\n
include("DH.jl")\n
include("kinematics.jl")\n
include("robotplot.jl")\n
include("Frames.jl")\n
include("read_log.jl")\n
include("csv2mat.jl")\n
include("posDepFric.jl")\n

The module includes a submodule, Frames, which is aimed at replacing the Nikon K600 software. It supports creation of frames, simple projections, fitting of planes, lines etc. and has a number of plotting options. It must be separately imported with `using Robotlib.Frames`

The module includes a submodule, Calibration, which includes a number of calibration routines. It must be separately imported with `using Robotlib.Calibration`

Usage:
fkine, ikine, jacobian = get_kinematic_functions("yumi")
data = orcalog2mat(pathopen, pathsave)
q = getData(\"robot_0.\*posRawAbs\", data, 1, removeNaN = false)

For YuMi, joint angles `q` must be converted to logical order using e.g. abb2logical!(q)
You must also consider the base transform of YuMi
"""
module Robotlib
include("utils.jl")
include("DH.jl")
include("kinematics.jl")
include("robotplot.jl")
include("Frames.jl")
include("read_log.jl")
include("csv2mat.jl")
include("posDepFric.jl")

module Calibration
using ..Robotlib
import Robotlib.xii
import Robotlib.Ai
include("calibLPOE.jl")
include("calibForce.jl")
export calibLPOE, calibLPOEdual, calibForce, calibPOE_offsets_from_points

end

export skewcoords, twistcoords, skew, skew4, expω, expξ, expξ2, expξ!, logT, logR, trinv, isrot, isse3, Rangle, conformize, DH2twistsPOE, DH2twistsLPOE, dh2Tn, toOrthoNormal, toOrthoNormal!, rpy2R, Quaternion, xyθ, smartDiff, R2rpy, traj2quat, centralDiff
export fkinePOE, fkineLPOE, ikinePOE, dh2Tn, jacobianPOE, jacobianPOEikine, jacobian, get_kinematic_functions
export trajplot, trajplot3, plot3smart
export DH, DH2400, DHYuMi, DH7600, DHtest, abb2logical!, logical2abb!, abb2logical, logical2abb
export csv2mat, orcalog2mat, getData, readmat
export frictionRBFN, getCenters

end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function annotatedash(pos, dashes, offset)
    @layer  begin
        translate(pos)
        @layer begin
            translate(0, 15)
            fontsize(3)
            sethue("black")
            x = 0
            for n in dashes
                text(string(n), O + (-offset + x + n/2, 0), halign=:center)
                x += n
            end
            fontsize(6)
            text("offset $offset", O + (-10, -10), halign=:right)
        end
        @layer begin
            translate(0, 2.5)
            setdash(dashes, offset)
            setline(0.5)
            sethue("grey60")
            line(O + (0, 6), O +  (3 * sum(dashes), 6), :stroke)
            sethue("red")
            setline(5)
            setdash(dashes, offset)
            line(O, O + (3 * sum(dashes), 0), :stroke)
        end
    end
end

function dash_tests(fname)
    pagewidth, pageheight = 800, 600
    Drawing(pagewidth, pageheight, fname)
    origin()
    background("ivory")
    setline(3)
    dashes = [30.0,  # ink
        10.0,  # skip
        5.0,  # ink
        15.0,  # skip
        20.0,  # ink
        4. # skip
        ]
    offset = -50.0
    sethue("black")
    text("Pattern is $(dashes)", boxtopcenter(BoundingBox() * 0.9),     halign=:center)
    sethue("blue")
    for y in -pageheight/2 + 50:30:pageheight/2 - 50
        randomhue()
        annotatedash(O + (-200, y), dashes, offset)
        offset += 10
    end
    @test finish() == true
end

fname = "dash-tests.svg"
dash_tests(fname)
println("...finished dash tests: output in $(fname)")


text:
"""
    minumum_operations(str_1::String, str_2::String)

The minimum operations algorithm explores how many string edits (`insert`, `remove`,
`replace`) are needed to convert str_1 to str_2.

# Arguments
- `str_1::String`: String 1 of length `m`
- `str_2::String`: String 2 of length `n`

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> str1 = "sunday"
julia> str2 = "saturday"
julia> ClassicAlgorithmsCollections.minumum_operations(str1,str2)
3
```
"""
function minumum_operations(str_1::String, str_2::String)
    m = length(str_1)
    n = length(str_2)
    if m == 1
        return n
    end
    if n == 1
        return m
    end

    if str_1[m] == str_2[n]
        return minumum_operations(str_1[1:m-1], str_2[1:n-1])
    end

    insert = minumum_operations(str_1[1:m], str_2[1:n-1])
    remove = minumum_operations(str_1[1:m-1], str_2[1:n])
    replace = minumum_operations(str_1[1:m-1], str_2[1:n-1])

    if insert < remove < replace || insert < replace < remove
        return insert
    elseif remove < insert < replace || remove < replace < insert
        return remove
    end
    return replace
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
"""
    polar(angle, radius)

Creates an `EChart` where region between points connected by a line, with polar coordinates.

## Methods
```julia
polar(angle::AbstractVector{<:Union{Missing, Real}}, radius::AbstractVector{<:Union{Missing, Real}})
polar(angle::AbstractVector{<:Union{Missing, Real}},radius::AbstractArray{<:Union{Missing, Real},2})
```

## Arguments
* `splitNumber::Int = 12` : number of splits shown in axis
* `showSymbol::Bool = false` : show the dots as well as the line?
* `legend::Bool` : display legend?
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `polar`, such as displaying a legend when two or more series present.

## Examples
```julia
t = 0:360
angle = [a / 180 * pi for a in t]
radius = [(sin(2t) * cos(2t)) for t in angle]
p = polar(angle, radius)
```
"""
function polar(angle::AbstractVector{<:Union{Missing, Real}},
               radius::AbstractVector{<:Union{Missing, Real}};
            splitNumber::Int = 12,
            showSymbol::Bool = false,
            legend::Bool = false,
            kwargs...)

    ec = newplot(kwargs, ec_charttype = "polar")

    ec.polar = [Polar()]
    ec.angleAxis = AngleAxis(_type = "value", splitNumber = splitNumber)
    ec.radiusAxis = RadiusAxis(min = 0)
    ec.series = [XYSeries(name = "Series 1", coordinateSystem = "polar", _type = "line", showSymbol = showSymbol, data = arrayofarray(radius,angle))]

    #Add legend if requested
    legend ? legend!(ec) : nothing

    return ec

end

function polar(angle::AbstractVector{<:Union{Missing, Real}},
               radius::AbstractArray{<:Union{Missing, Real},2};
            splitNumber::Int = 12,
            showSymbol::Bool = false,
            legend::Bool = true,
            kwargs...)

    #Instantiate polar graph with first series
    ec = polar(angle, radius[:,1]; splitNumber = splitNumber, showSymbol = showSymbol, kwargs...)

    #Append remaining series
	for i in 2:size(radius)[2]
		push!(ec.series, XYSeries(coordinateSystem = "polar", _type = "line", showSymbol = showSymbol, data = arrayofarray(radius[:,i],angle)))
	end

    #Add seriesnames for legend
    seriesnames!(ec)

    #Add legend if requested
    legend ? legend!(ec) : nothing

    return ec

end


text:
function overtime(img)
    img_result = zeros(size(img[:, :, 1]))

    for i = 1:length(img[:, 1, 1])
        for j = 1:length(img[1, :, 1])
            img_result[i, j] = sum(img[i, j, :])
        end
    end
    return img_result
end

function speed_cal(angle, line)
    band_speed = zeros(150, 2)
    current_point = 1
    band_num = 1
    for i = 2:length(line)-1
        if ((line[i-1] >= 0.02) && (line[i] <= 0.02)) ||
           ((line[i-1] <= 0.02) && (line[i] >= 0.02))
            dis = i - current_point
            band_speed[band_num, 2] = 34 * band_num
            band_speed[band_num, 1] = (i - 100) / 100 * angle
            current_point = i
            band_num += 1
        end
    end
    return band_speed
end

function smoothing(img, freq)
    designmethod = Butterworth(5)
    ff = digitalfilter(Lowpass(freq), designmethod)
    result = filtfilt(ff, img)
    return result
end

function fwhm(data)
    max_d = maximum(data) - data[end]
    mid_d = argmax(data)[1]
    half_value = max_d / 2 + data[end]
    sigma = zeros(1, 2)
    for i = 1:mid_d-1
        if data[mid_d-i] <= half_value && data[mid_d-i+1] >= half_value
            sigma[1] = mid_d - i
            break
        end
    end
    for i = mid_d:length(data)-1
        if data[i] <= half_value && data[i+1] >= half_value
            sigma[2] = i
            break
        end
    end
    return (sigma[2] - sigma[1])
end

function clearzero(data)
    for i = 1:length(data[:, 1])
        if data[end-i+1, 2] != 0.0
            return data[1:end-i+1, :]
            break
        end
    end
end

#print("Date: ")
#date=chomp(readline())
##date="20190315"
#
#cd(@sprintf("/mnt/jchang_005/jchang/%s/",date))
#files=readdir()
#filelist=filter(x->occursin(".imagine",x),files)
#@show filelist
#
#print("Filename: ")
#files=readdir()
#filename = chomp(readline())
#
#@printf("angle: ")
#A = chomp(readline())
#A = parse(Int,A)
#
#
#scene1 = Scene()
#scene2 = Scene()
#exp = load(@sprintf("%s_1.tif", filename))
##img1 = exp[:,:,t]
##img2 = convert(Array{N0f16}, img1)
#img3 = overtime(exp)
#idx = argmax(img3)
##scene1 = Makie.plot(smoothing(img2[idx,:],0.05))
#line_c = exp[idx[1],idx[2],:]
#line_s = smoothing(line_c,0.9)
#band_speed = speed_cal(A,line_c)
#band_pos = clearzero(band_speed)
#
#scene = Makie.lines(band_pos)
#
#time_pos = zeros(length(band_pos[:,1]),2)
#for i in 1:length(band_pos[:,1])
#        img3 = exp[:,:,Int(floor((band_pos[i,1]+A)*100/A))]
#        img4 = convert(Array{N0f16}, img3)
#        time_c = argmax(img4)[2]
#        time_pos[i,2] = std(img4[:,time_c])
#	time_pos[i,1] = band_pos[i,1]
#end
##time_pos = smoothing(time_pos,0.3)