text:
module GeckoScreen

using DataFrames
using CSV
using Query
using Statistics
using Dates

export readGeCKOlibrary, createScreenDF

"""    readGeCKOlibrary()

Creates a DataFrame of the GeCKO libraries. Created from CSV file of sequences belonging to "Human CRISPR Knockout Pooled LIbrary (GeCKO v2): can be found: https://www.addgene.org/pooled-library/zhang-human-gecko-v2/.

"""
function readGeCKOlibrary()

  ## These library files have 3 columns, 1 with the gene name, one with the UIDs, and the thirds has the specific sequences
  fields = ["","",""]
  collectedFields = ["","",""];
  for line in eachline("human_geckov2_library_a_09mar2015_v2.csv")
      fields = split(line, ",")
      collectedFields = hcat(collectedFields, fields)
  end

  genes = collectedFields[1,3:end];
  UIDs = collectedFields[2,3:end];
  seqs = collectedFields[3,3:end];

  ## Create DataFrame with first half of sequences
  af = DataFrame(gene=genes, UID=UIDs, seq=seqs);

  fieldsB = ["","",""]
  collectedFieldsB = ["","",""];
  for line in eachline("human_geckov2_library_b_09mar2015_v2.csv")
      fieldsB = split(line, ",")
      collectedFieldsB = hcat(collectedFieldsB, fieldsB)
  end

  genesB = collectedFieldsB[1,3:end];
  UIDsB = collectedFieldsB[2,3:end];
  seqsB = collectedFieldsB[3,3:end];

  ## Create DataFrame with second half of sequences
  bf = DataFrame(gene=genesB, UID=UIDsB, seq=seqsB);

  ## Combine all sequences into same DataFrame and make sure :seq is the correct type
  gmap = append!(af,bf)
  gmap[:seq] = convert(Array{Union{Missing, String},1}, gmap[:seq]) Array{AbstractString,1}

  return gmap

end

### create same function that will load libraries if given a file

struct RNAScreen
  drug::String
  day::Int
  replicate::Int
  path::String
end

function createScreenDF(files)

  for rnas in files
      df = CSV.read(rnas.path, header=[:rna, :count])
      df[:drug] = rnas.drug
      df[:day] = rnas.day
      df[:replicate] = rnas.replicate
      CSV.write(string("RNAScreens_RIGER_",Dates.today(),".csv"), df, append=true)
  end

  df = CSV.read(string("RNAScreens_RIGER_",Dates.today(),".csv"), header=[:rna, :count, :drug, :day, :replicate])

  return df

end


text:
using Luxor, Test

# test all points the same
pt1 = pt2 = pt3 = pt4 = O
flag, _ = intersectionlines(pt1, pt1, pt1, pt1)
@test flag == false

# test all points the same but not 0.0
pt1 = pt2 = pt3 = pt4 = Point(0.5, 0.75)
flag, _ = intersectionlines(pt1, pt1, pt1, pt1)
@test flag == false

#first pair of points the same
pt1 = pt2 = O
pt3 = Point(20, 20)
pt4 = Point(30, 20)
flag, _ = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

#last pair of points the same
pt1 = Point(20, 20)
pt2 = Point(30, 20)
pt3 = pt4 = O
flag, _ = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

#first point same as one of second pair
pt1 = Point(20, 20)
pt2 = Point(30, 20)
pt3 = Point(20, 20)
pt4 = Point(40, 50)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true
@test ip == Point(20.0, 20.0)

# first point same as one of second pair, crossing only
pt1 = Point(20, 20)
pt2 = Point(30, 20)
pt3 = Point(20, 20)
pt4 = Point(40, 50)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == true
@test ip == Point(20.0, 20.0)

#testing overlapping collinear horizontal lines
pt1 = Point(0, 0)
pt2 = Point(100, 0)
pt3 = Point(50, 0)
pt4 = Point(200, 0)
flag, _ = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

# testing overlapping collinear horizontal lines, with crossingonly
pt1 = Point(0, 0)
pt2 = Point(100, 0)
pt3 = Point(50, 0)
pt4 = Point(200, 0)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

# testing shared endpoints horizontal crossingonly = true
pt1 = Point(0, 0)
pt2 = Point(200, 0)
pt3 = Point(200, 0)
pt4 = Point(100, 40)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == true
@test ip == Point(200.0, 0.0)

# testing shared endpoints vertical
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 0)
pt4 = Point(0, 200)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true
@test ip == O

# testing shared endpoints vertical crossingonly = true
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 10)
pt4 = Point(0, 100)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == true
@test ip == Point(0.0, 100.0)

# testing nonoverlapping vertical lines
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 100)
pt4 = Point(0, 200)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true
@test ip == Point(0.0, 100.0)

# testing nonoverlapping vertical lines
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 200)
pt4 = Point(0, 300)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

# testing overlapping vertical lines
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 50)
pt4 = Point(0, 150)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

# testing overlapping vertical lines crossingonly true"
pt1 = Point(0, 0)
pt2 = Point(0, 100)
pt3 = Point(0, 50)
pt4 = Point(0, 150)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

# testing intersecting line segments not crossing"
pt1 = Point(0, 0)
pt2 = Point(150, -150)
pt3 = Point(150, 0)
pt4 = Point(0, -150)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true
@test ip == Point(75.0, -75.0)

# testing nearly vertical lines
pt1 = Point(0, 0)
pt2 = Point(0, -1500)
pt3 = Point(1, 0)
pt4 = Point(0.99, -1500)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true

# testing nearly vertical lines crossing only
pt1 = Point(0, 0)
pt2 = Point(0, -15000)
pt3 = Point(1, 0)
pt4 = Point(0.999999, -150000)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

# testing nearly vertical lines, they meet eventually
pt1 = Point(0, 0)
pt2 = Point(0, -15000)
pt3 = Point(1, 0)
pt4 = Point(0.999999, -150000)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test isapprox(ip.x, 0.0)
@test flag == true

# parallel vertical lines
pt1 = Point(0, 0)
pt2 = Point(0, -1500)
pt3 = Point(1, 0)
pt4 = Point(1, -1500)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

# parallel vertical lines crossingonly
pt1 = Point(0, 0)
pt2 = Point(0, -1500)
pt3 = Point(1, 0)
pt4 = Point(1, -1500)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

# parallel horizontal lines
pt1 = Point(-1000, 0)
pt2 = Point(1000,  0)
pt3 = Point(-1000, -10)
pt4 = Point(1000,  -10)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == false

# parallel horizontal lines crossingonly
pt1 = Point(-1000, 0)
pt2 = Point(1000,  0)
pt3 = Point(-1000, -10)
pt4 = Point(1000,  -10)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

# nearly parallel horizontal lines (they meet eventually)
pt1 = Point(-1000, 0)
pt2 = Point(1000,  0)
pt3 = Point(-1000, -10)
pt4 = Point(1000,  -9.99)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4)
@test flag == true
@test ip.x > 1000.0
@test iszero(ip.y)

# nearly parallel horizontal lines crossing only
pt1 = Point(-100, 0)
pt2 = Point(100,  0)
pt3 = Point(-100, -10)
pt4 = Point(100,  -9.99)
flag, ip = intersectionlines(pt1, pt2, pt3, pt4, crossingonly=true)
@test flag == false

ext: jl


text:
"""
    function addCleanInteraction_!(p::Program_,platform::String)

Generate the functions related with Cleaning Interaction Updates.
"""
function addCleanInteraction_!(p::Program_,platform::String)

    if "UpdateLocalInteraction" in keys(p.agent.declaredUpdates)

        #Construct cleanup function
        up = quote end

        for i in p.agent.declaredSymbols["Local"]

            code =  postwalk(x->@capture(x, c_.g_ += f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , p.agent.declaredUpdates["UpdateLocalInteraction"])
            code =  postwalk(x->@capture(x, c_.g_ -= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ *= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ /= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ \= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ÷= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ %= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ^= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ &= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ |= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ⊻= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ <<= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)

            if inexpr(code,:ARGS_)
                pos = findfirst(p.agent.declaredSymbols["Local"] .== i)
                push!(up.args,:(localV[ic1_,$pos]=0))
            end

        end

        for i in p.agent.declaredSymbols["Identity"]

            code =  postwalk(x->@capture(x, c_.g_ += f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , p.agent.declaredUpdates["UpdateLocalInteraction"])
            code =  postwalk(x->@capture(x, c_.g_ -= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ *= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ /= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ \= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ÷= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ %= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ^= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ &= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ |= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ⊻= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ <<= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)

            if inexpr(code,:ARGS_)
                pos = findfirst(p.agent.declaredSymbols["Identity"] .== i)
                push!(up.args,:(identityV[ic1_,$pos]=0))
            end

        end

        if up != quote end
            fclean = simpleFirstLoopWrapInFunction_(platform,:cleanLocalInteraction_!,up)
            push!(p.declareF.args,fclean)
        end
    end

    if "UpdateInteraction" in keys(p.agent.declaredUpdates)

        #Construct cleanup function
        up = quote end

        for i in p.agent.declaredSymbols["Local"]

            code =  postwalk(x->@capture(x, c_.g_ += f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , p.agent.declaredUpdates["UpdateInteraction"])
            code =  postwalk(x->@capture(x, c_.g_ -= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ *= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ /= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ \= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ÷= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ %= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ^= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ &= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ |= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ⊻= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ <<= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)

            if inexpr(code,:ARGS_)
                pos = findfirst(p.agent.declaredSymbols["Local"] .== i)
                push!(up.args,:(localV[ic1_,$pos]=0))
            end

        end

        for i in p.agent.declaredSymbols["Identity"]

            code =  postwalk(x->@capture(x, c_.g_ += f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , p.agent.declaredUpdates["UpdateInteraction"])
            code =  postwalk(x->@capture(x, c_.g_ -= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ *= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ /= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ \= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ÷= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ %= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ^= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ &= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ |= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ ⊻= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ >>= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)
            code =  postwalk(x->@capture(x, c_.g_ <<= f_) && c == i && g in INTERACTIONSYMBOLS ? :ARGS_ : x , code)

            if inexpr(code,:ARGS_)
                pos = findfirst(p.agent.declaredSymbols["Identity"] .== i)
                push!(up.args,:(identityV[ic1_,$pos]=0))
            end

        end

        if up != quote end
            fclean = simpleFirstLoopWrapInFunction_(platform,:cleanInteraction_!,up)
            push!(p.declareF.args,fclean)
        end
    end

    return nothing
end


text:
# Minimal sanity checks for MakieLayout
@testset "Layoutables constructors" begin
    scene, layout = layoutscene()
    ax = layout[1, 1] = Axis(scene)
    cb = layout[1, 2] = Colorbar(scene)
    gl2 = layout[2, :] = MakieLayout.GridLayout()
    bu = gl2[1, 1] = Button(scene)
    sl = gl2[1, 2] = Slider(scene)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = Legend(scene, [scat], ["scatter"])

    to = gl2[1, 4] = Toggle(scene)
    te = layout[0, :] = Label(scene, "A super title")
    me = layout[end + 1, :] = Menu(scene, options=["one", "two", "three"])
    tb = layout[end + 1, :] = Textbox(scene)
    is = layout[end + 1, :] = IntervalSlider(scene)
    @test true
end

@testset "deleting from axis" begin
    f = Figure()
    ax = Axis(f[1, 1])
    sc = scatter!(ax, randn(100, 2))
    li = lines!(ax, randn(100, 2))
    hm = heatmap!(ax, randn(20, 20))
    @test length(ax.scene.plots) == 3
    delete!(ax, sc)
    @test length(ax.scene.plots) == 2
    @test sc ∉ ax.scene.plots
    empty!(ax)
    @test isempty(ax.scene.plots)
end


text:
import Base.Broadcast: broadcasted

"""
    .+(x1::Var, x2::Var)
"""
function broadcasted(::typeof(+), x1::Var, x2::Var)
    ydata = x1.data .+ x2.data
    Var(ydata, ∇broadcasted!, (+,x1,x2))
end
broadcasted(::typeof(+), x1::Node, x2) = Node(broadcasted, (+,x1,x2))
broadcasted(::typeof(+), x1, x2::Node) = Node(broadcasted, (+,x1,x2))

function ∇broadcasted!(y::Var, ::typeof(+), x1::Var, x2::Var)
    isnothing(x1.grad) || ∇broadcast_plus!(y.grad, x1.grad)
    isnothing(x2.grad) || ∇broadcast_plus!(y.grad, x2.grad)
end

function ∇broadcast_plus!(gy::UniArray{T}, gx::UniArray{T}) where T
    if length(gy) == length(gx)
        addto!(gx, gy)
    else
        dims = ()
        for i = 1:ndims(gy)
            if size(gx,i) == 1 && size(gy,i) > 1
                dims = tuple(dims..., i)
            end
        end
        addto!(gx, sum(gy,dims=dims))
    end
end

"""
    .-(x1::Var, x2::Var)
"""
function broadcasted(::typeof(-), x1::Var, x2::Var)
    ydata = x1.data .- x2.data
    Var(ydata, ∇broadcasted!, (-,x1,x2))
end
broadcasted(::typeof(-), x1::Node, x2) = Node(broadcasted, (-,x1,x2))
broadcasted(::typeof(-), x1, x2::Node) = Node(broadcasted, (-,x1,x2))

function ∇broadcasted!(y::Var, ::typeof(-), x1::Var, x2::Var)
    isnothing(x1.grad) || ∇broadcast_plus!(y.grad, x1.grad)
    isnothing(x2.grad) || ∇broadcast_minus!(y.grad, x2.grad)
end

function ∇broadcast_minus!(gy::UniArray{T}, gx::UniArray{T}) where T
    if length(gy) == length(gx)
        axpy!(T(-1), gy, gx)
    else
        dims = ()
        for i = 1:ndims(gy)
            if size(gx,i) == 1 && size(gy,i) > 1
                dims = tuple(dims..., i)
            end
        end
        axpy!(T(-1), sum(gy,dims=dims), gx)
    end
end

"""
    .*(x1::Var, x2::Var)
"""
function broadcasted(::typeof(*), x1::Var, x2::Var)
    ydata = x1.data .* x2.data
    Var(ydata, ∇broadcasted!, (*,x1,x2))
end
broadcasted(::typeof(*), x1::Node, x2) = Node(broadcasted, (*,x1,x2))
broadcasted(::typeof(*), x1, x2::Node) = Node(broadcasted, (*,x1,x2))

function ∇broadcasted!(y::Var, ::typeof(*), x1::Var, x2::Var)
    isnothing(x1.grad) || ∇dottimes!(y.grad, x2.data, x1.grad)
    isnothing(x2.grad) || ∇dottimes!(y.grad, x1.data, x2.grad)
end

function ∇dottimes!(gy::UniArray{T}, x2::UniArray{T}, gx1::UniArray{T}) where T
    g = gy .* x2
    if length(gx1) == length(g)
        addto!(gx1, g)
    else
        dims = ()
        for i = 1:ndims(gy)
            if size(gx1,i) == 1 && size(g,i) > 1
                dims = tuple(dims..., i)
            end
        end
        addto!(gx1, sum(g,dims=dims))
    end
end


text:
############
# Starring #
############

@api_default function stargazers(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/stargazers"; options...)
    return map(Owner, results), page_data
end

@api_default function starred(api::GitHubAPI, user; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(user))/starred"; options...)
    return map(Repo, results), page_data
end

@api_default star(api::GitHubAPI, repo; options...) = gh_put(api, "/user/starred/$(name(repo))"; options...)

@api_default unstar(api::GitHubAPI, repo; options...) = gh_delete(api, "/user/starred/$(name(repo))"; options...)

############
# Watching #
############

@api_default function watchers(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/subscribers"; options...)
    return map(Owner, results), page_data
end

@api_default function watched(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/subscriptions"; options...)
    return map(Repo, results), page_data
end

@api_default watch(api::GitHubAPI, repo; options...) = gh_put(api, "/repos/$(name(repo))/subscription"; options...)

@api_default unwatch(api::GitHubAPI, repo; options...) = gh_delete(api, "/repos/$(name(repo))/subscription"; options...)


text:
module TestWithCollect

using Test
using LazyGroupBy

dataset = [
    "map(identity, _)" => g -> map.(identity, g),
    "map(iseven, _)" => g -> map.(iseven, g),
    "sum(_)" => g -> sum.(g),
    "sum(x -> 2x, _)" => g -> sum.(x -> 2x, g),
    "prod(_)" => g -> prod.(g),
    "prod(x -> 2 + x, _)" => g -> prod.(x -> 2 + x, g),
    "length(_)" => g -> length.(g),
    "any(isodd, _)" => g -> any.(isodd, g),
    "all(isodd, _)" => g -> all.(isodd, g),
    "minimum(_)" => g -> minimum.(g),
    "minimum(x -> 2x, _)" => g -> minimum.(x -> 2x, g),
    "maximum(_)" => g -> maximum.(g),
    "maximum(x -> 2x, _)" => g -> maximum.(x -> 2x, g),
    "extrema(_)" => g -> extrema.(g),
    "extrema(x -> 2x, _)" => g -> extrema.(x -> 2x, g),
]

@testset "$label" for (label, reducer) in dataset
    g = grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])
    d = collect.(g)
    @test reducer(g) == Dict(zip(keys(d), reducer(values(d))))
end

end  # module


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to create bilateral Head-Ries index (symmetric bilateral trade costs)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    head_ries_index(Z::Matrix, F::Matrix, θ, N::Integer, S::Integer)

The function computes bilateral trade costs according to Head-Ries, i.e. symmetric inter country-industry trade costs and
    zero intra country trade costs. Trade costs are computed for both intermediate input and final demand. The trade elasticity can
    take the form of a uniform country-industry value or origin country-industry trade elasticity.

# Arguments
- `Z::Matrix`: NS×NS, origin country-industry destination country-industry intermediate demand matrix Z.
- `F::Matrix`: NS×N, origin country-industry destination country final demand matrix F.
- `θ`: assumption for the effective trade elasticity, can either be supplied as a simple value or a S×N matrix.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `τ_HR_Z::Matrix{Float64}`: NS×NS, Head-Reis Index, symmetric bilateral country-industry intermediate input trade costs
        with intra-country trade assumed to have no trade costs.
- `τ_HR_F::Matrix{Float64}`: NS×N,  Head-Reis Index, symmetric bilateral country-industry final demand trade costs
    with intra-country trade assumed to have no trade costs.

# Examples
```julia-repl
julia> τ_HR_Z, τ_HR_F = head_ries_index(Z, F, θ, N, S)
```
"""

function head_ries_index(Z::Matrix, F::Matrix, θ, N::Integer, S::Integer)

    if length(θ) == 1
        θ = fill(θ, N*S) # NS×1, work with country-industry elasticities
        θ = reshape(θ, S, N) # S×N
    elseif size(θ) == (S, N)
    else
        println("The effective trade elasticity θ is not supplied in the right format: either a simple number or a S×N matrix.")
    end

    Z_HR = ifelse.(Z .== 0.0, 1e-18, Z) # otherwise we would have infinite trade costs
    τ_HR_Z = zeros(N*S, N*S) # NS×NS

    for ctry1 in 1:N
        for ctry2 in 1:N
            for ind1 in 1:S
                for ind2 in 1:S
                    if ctry1 == ctry2
                        τ_HR_Z[ind1+(ctry1-1)*S, ind2+(ctry2-1)*S] = 1.0
                    else
                        Z1122 = Z_HR[ind1+(ctry1-1)*S, ind2+(ctry2-1)*S]
                        Z2122 = Z_HR[ind1+(ctry2-1)*S, ind2+(ctry2-1)*S]
                        Z2112 = Z_HR[ind1+(ctry2-1)*S, ind2+(ctry1-1)*S]
                        Z1112 = Z_HR[ind1+(ctry1-1)*S, ind2+(ctry1-1)*S]
                        τ_temp = ((Z1122*Z2112)/(Z2122*Z1112))^(-1/(2*θ[ind1,ctry1]))
                        τ_HR_Z[ind1+(ctry1-1)*S, ind2+(ctry2-1)*S] = τ_temp
                        τ_HR_Z[ind1+(ctry2-1)*S, ind2+(ctry1-1)*S] = τ_temp
                    end
                end
            end
        end
    end

    F_HR = ifelse.(F .== 0.0, 1e-18, F)
    τ_HR_F = zeros(N*S, N) # NS×N

    for ctry1 in 1:N
        for ctry2 in 1:N
            for ind1 in 1:S
                if ctry1 == ctry2
                    τ_HR_F[ind1+(ctry1-1)*S, ctry2] = 1.0
                else
                    F12 = F_HR[ind1+(ctry1-1)*S, ctry2]
                    F22 = F_HR[ind1+(ctry2-1)*S, ctry2]
                    F21 = F_HR[ind1+(ctry2-1)*S, ctry1]
                    F11 = F_HR[ind1+(ctry1-1)*S, ctry1]
                    τ_temp = ((F12*F21)/(F11*F22))^(-1/(2*θ[ind1,ctry1]))
                    τ_HR_F[ind1+(ctry1-1)*S, ctry2] = τ_temp
                    τ_HR_F[ind1+(ctry2-1)*S, ctry1] = τ_temp
                end
            end
        end
    end

    return τ_HR_Z, τ_HR_F
end


text:
# ---
# title: 295. Find Median from Data Stream
# id: problem295
# author: zhwang
# date: 2022-02-26
# difficulty: Hard
# categories: Heap, Design
# link: <https://leetcode.com/problems/find-median-from-data-stream/description/>
# hidden: true
# ---
#
# Median is the middle value in an ordered integer list. If the size of the list
# is even, there is no middle value. So the median is the mean of the two middle
# value.
#
# For example,
#
# `[2,3,4]`, the median is `3`
#
# `[2,3]`, the median is `(2 + 3) / 2 = 2.5`
#
# Design a data structure that supports the following two operations:
#
#   * void addNum(int num) - Add a integer number from the data stream to the data structure.
#   * double findMedian() - Return the median of all elements so far.
#
#
#
# **Example:**
#
#
#
#     addNum(1)
#     addNum(2)
#     findMedian() -> 1.5
#     addNum(3)
#     findMedian() -> 2
#
#
#
#
# **Follow up:**
#
#   1. If all integer numbers from the stream are between 0 and 100, how would you optimize it?
#   2. If 99% of all integer numbers from the stream are between 0 and 100, how would you optimize it?
#
#
## @lc code=start
using LeetCode

using DataStructures

Base.@kwdef struct MedianFinder
    maxheap::BinaryMaxHeap{Int} = BinaryMaxHeap{Int}()
    minheap::BinaryMinHeap{Int} = BinaryMinHeap{Int}()
end
## function Base.show(io::IO, heap::MedianFinder)
##     return print(io, "$(heap.maxheap.valtree):max -- min:$(reverse(heap.minheap.valtree))")
## end

function add_num!(heap::MedianFinder, num::Int)
    ## maximal element of maxheap <= minimal element of minheap
    ## length of maxheap <= length of minheap
    hi, ha = heap.minheap, heap.maxheap
    if length(ha) < length(hi)
        ## e.g. [2]--[4,5] =6> [2]--[4,5,6] => [2,4]--[5,6]
        push!(hi, num)
        push!(ha, pop!(hi))
    else
        ## e.g. [4]--[5] =2> [2,4]--[5] => [2]--[4,5]
        push!(ha, num)
        push!(hi, pop!(ha))
    end
    return nothing
end

function find_median(heap::MedianFinder)::Float64
    hi, ha = heap.minheap, heap.maxheap
    return length(hi) == length(ha) ? (first(hi) + first(ha)) / 2 : first(hi)
end
## @lc code=end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antddivider

"""
    'feffery'_antddivider(;kwargs...)
    'feffery'_antddivider(children::Any;kwargs...)
    'feffery'_antddivider(children_maker::Function;kwargs...)

An AntdDivider component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `className` (String; optional)
- `direction` (a value equal to: 'horizontal', 'vertical'; optional)
- `fontFamily` (String; optional)
- `fontSize` (String; optional)
- `fontStyle` (String; optional)
- `fontWeight` (String; optional)
- `innerTextOrientation` (a value equal to: 'left', 'center', 'right'; optional)
- `isDashed` (Bool; optional)
- `lineColor` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function 'feffery'_antddivider(; kwargs...)
        available_props = Symbol[:children, :id, :className, :direction, :fontFamily, :fontSize, :fontStyle, :fontWeight, :innerTextOrientation, :isDashed, :lineColor, :loading_state]
        wild_props = Symbol[]
        return Component("'feffery'_antddivider", "AntdDivider", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antddivider(children::Any; kwargs...) = 'feffery'_antddivider(;kwargs..., children = children)
'feffery'_antddivider(children_maker::Function; kwargs...) = 'feffery'_antddivider(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
"""
    AbstractBridge

Subtype of [`MathOptInterface.Bridges.AbstractBridge`](@ref) for objective
bridges.
"""
abstract type AbstractBridge <: MOIB.AbstractBridge end

"""
    bridge_objective(BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
                     model::MOI.ModelLike,
                     func::MOI.AbstractScalarFunction)

Bridge the objective function `func` using bridge `BT` to `model` and returns
a bridge object of type `BT`. The bridge type `BT` should be a concrete type,
that is, all the type parameters of the bridge should be set. Use
[`concrete_bridge_type`](@ref) to obtain a concrete type for a given function
type.
"""
function bridge_objective(
    ::Type{<:AbstractBridge},
    ::MOI.ModelLike,
    func::MOI.AbstractScalarFunction,
)
    return throw(
        MOI.UnsupportedAttribute(MOI.ObjectiveFunction{typeof(func)}()),
    )
end

"""
    function MOI.set(model::MOI.ModelLike, ::MOI.ObjectiveSense,
                     bridge::MOI.Bridges.Objective.AbstractBridge,
                     sense::MOI.ObjectiveSense)

Informs `bridge` that the objective sense is changed to `sense`. If changing
the objective sense is not supported, the bridge should not implement this
method.
"""
function MOI.set(
    ::MOI.ModelLike,
    ::MOI.ObjectiveSense,
    bridge::AbstractBridge,
    ::MOI.OptimizationSense,
)
    return throw(
        ArgumentError(
            "Objective bridge of type `$(typeof(bridge))` does not support" *
            " modifying the objective sense. As a workaround, set the sense to" *
            " `MOI.FEASIBILITY_SENSE` to clear the objective function and" *
            " bridges.",
        ),
    )
end

"""
    function MOI.get(model::MOI.ModelLike, attr::MOI.ObjectiveFunction,
                     bridge::MOI.Bridges.Objective.AbstractBridge)

Return the objective function object bridged by `bridge` for the model
`model`.
"""
function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ObjectiveFunction,
    bridge::AbstractBridge,
)
    return throw(
        ArgumentError(
            "ObjectiveFunction bridge of type `$(typeof(bridge))` does not" *
            " support getting the objective function.",
        ),
    )
end

"""
    function MOI.delete(model::MOI.ModelLike, bridge::MOI.Bridges.Objective.AbstractBridge)

Delete any variable or constraint added by `bridge`.
"""
function MOI.delete(::MOI.ModelLike, bridge::AbstractBridge)
    return throw(
        ArgumentError(
            "`MOI.delete` not implemented for `ObjectiveFunction` bridges of type `$(typeof(bridge))`",
        ),
    )
end

"""
    supports_objective_function(
        BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
        F::Type{<:MOI.AbstractScalarFunction})::Bool

Return a `Bool` indicating whether the bridges of type `BT` support bridging
objective functions of type `F`.
"""
function supports_objective_function(
    ::Type{<:AbstractBridge},
    ::Type{<:MOI.AbstractScalarFunction},
)
    return false
end

"""
    added_constrained_variable_types(BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
                                     F::Type{<:MOI.AbstractScalarFunction})

Return a list of the types of constrained variables that bridges of type `BT`
add for bridging objective functions of type `F`. This fallbacks to
`added_constrained_variable_types(concrete_bridge_type(BT, F))`
so bridges should not implement this method.
```
"""
function MOIB.added_constrained_variable_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractScalarFunction},
)
    return MOIB.added_constrained_variable_types(concrete_bridge_type(BT, F))
end

"""
    added_constraint_types(BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
                           F::Type{<:MOI.AbstractScalarFunction})

Return a list of the types of constraints that bridges of type `BT` add
for bridging objective functions of type `F`. This fallbacks to
`added_constraint_types(concrete_bridge_type(BT, S))`
so bridges should not implement this method.
"""
function MOIB.added_constraint_types(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractScalarFunction},
)
    return MOIB.added_constraint_types(concrete_bridge_type(BT, F))
end

"""
    set_objective_function_type(BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
                                F::Type{<:MOI.AbstractScalarFunction})

Return the type of objective function that bridges of type `BT` set
for bridging objective functions of type `F`. This fallbacks to
`set_objective_function_type(concrete_bridge_type(BT, F))`
so bridges should not implement this method.
"""
function MOIB.set_objective_function_type(
    BT::Type{<:AbstractBridge},
    F::Type{<:MOI.AbstractScalarFunction},
)
    return MOIB.set_objective_function_type(concrete_bridge_type(BT, F))
end

"""
    concrete_bridge_type(BT::Type{<:MOI.Bridges.Objective.AbstractBridge},
                         F::Type{<:MOI.AbstractScalarFunction})::DataType

Return the concrete type of the bridge supporting objective functions of type
`F`. This function can only be called if `MOI.supports_objective_function(BT, F)`
is `true`.
"""
function concrete_bridge_type(
    bridge_type::DataType,
    ::Type{<:MOI.AbstractScalarFunction},
)
    return bridge_type
end

function concrete_bridge_type(
    b::MOIB.AbstractBridgeOptimizer,
    F::Type{<:MOI.AbstractScalarFunction},
)
    return concrete_bridge_type(MOIB.bridge_type(b, F), F)
end


text:
using PyCall: pyimport_conda, pycall
using RCall, Conda

function installpypackage()
	try
		pyimport_conda("sklearn", "scikit-learn")
	catch
		try
			Conda.add("scikit-learn")
		catch
			println("scikit-learn failed to install")
		end
	end
end

function installrpackage(package::AbstractString)
	try
		rcall(:library,package,"lib=.libPaths()")
		#rcall(:library,package,"lib=Sys.getenv('R_LIBS_USER')")
	catch
		try
			R"dir.create(path = Sys.getenv('R_LIBS_USER'), showWarnings = FALSE, recursive = TRUE)"
			R"install.packages($package,lib=Sys.getenv('R_LIBS_USER'),repos='https://cloud.r-project.org',type='binary')"
		catch xerror
			println(xerror)
			println("package "*package*" failed to install")
		end
	end
end

function installrml()
	#packages=["caret", "earth","mda","e1071","gam","randomForest","nnet","kernlab","grid","MASS","pls"]
	#packages=["caret", "e1071","gam","randomForest"]
	packages=["caret","e1071","randomForest"]
	for pk in packages
		installrpackage(pk)
	end
end

installrml()
installpypackage()

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
"""
    ContainerIndices{N, Ax}

A `DataType` for storing the indexing information of arbitrary Julia arrays
and JuMP containers. This is to be used to enable
[`vectorize`](@ref InfiniteOpt.Collections.vectorize) and
[`unvectorize`](@ref InfiniteOpt.Collections.unvectorize). `ContainerIndices`
should be constructed via [`indices`](@ref InfiniteOpt.Collections.indices).
"""
struct ContainerIndices{N, Ax}
    indices::CartesianIndices{N, NTuple{N, Base.OneTo{Int}}}
    axes::Ax
end

## Make constructors for different array types
# Array
function ContainerIndices(
    arr::Array{T, N}
    )::ContainerIndices{N, Nothing} where {T, N}
    return ContainerIndices(CartesianIndices(arr), nothing)
end

# DenseAxisArray
function ContainerIndices(
    arr::JuMPC.DenseAxisArray{T, N, Ax}
    )::ContainerIndices{N, Ax} where {T, N, Ax}
    return ContainerIndices(CartesianIndices(arr), axes(arr))
end

# SparseAxisArray
function ContainerIndices(
    arr::JuMPC.SparseAxisArray{T, N, K}
    )::ContainerIndices{1, Vector{K}} where {T, N, K}
    axs = sort!(collect(keys(arr.data)))
    return ContainerIndices(CartesianIndices(axs), axs)
end

# Extend Base.:(==)
function Base.:(==)(i1::ContainerIndices, i2::ContainerIndices)::Bool
    return i1.indices == i2.indices && i1.axes == i2.axes
end

# Extend Base.size
function Base.size(inds::ContainerIndices)
    return size(inds.indices) # This should be avoided for SparseAxisArrays
end

# Extend Base.length
function Base.length(inds::ContainerIndices)::Int
    return length(inds.indices)
end

"""
    indices(arr)::Union{ContainerIndices, Nothing}

Return the indices of `arr` to aid with calls made to
[`vectorize`](@ref InfiniteOpt.Collections.vectorize) and
[`unvectorize`](@ref InfiniteOpt.Collections.unvectorize).
"""
function indices(arr)::Nothing
    return
end

# Array
function indices(arr::Array{T, N})::ContainerIndices{N, Nothing} where {T, N}
    return ContainerIndices(arr)
end

# DenseAxisArray
function indices(
    arr::JuMPC.DenseAxisArray{T, N, Ax}
    )::ContainerIndices{N, Ax} where {T, N, Ax}
    return ContainerIndices(arr)
end

# SparseAxisArray
function indices(
    arr::JuMPC.SparseAxisArray{T, N, K}
    )::ContainerIndices{1, Vector{K}} where {T, N, K}
    return ContainerIndices(arr)
end

"""
    vectorize(arr, [inds::Union{Nothing, ContainerIndices}])::Vector

Return the vectorized version of `arr`. If `arr` is not an `AbstractArray` then
it is returned unchanged. Optionally, the indices `inds` can be given to enhance
the efficiency if they were already made extracted via the
[`indices`](@ref InfiniteOpt.Collections.indices) method.
"""
function vectorize(arr::T)::T where {T}
    return arr
end

# Single argument Vector
function vectorize(arr::Vector{T})::Vector{T} where {T}
    return arr
end

# Single argument Array/DenseAxisArray
function vectorize(arr::AbstractArray{T})::Vector{T} where {T}
    if isempty(arr)
        return T[]
    else
        return reduce(vcat, arr)
    end
end

# Single argument SparseAxisArray
function vectorize(arr::JuMPC.SparseAxisArray{T})::Vector{T} where {T}
    ks = sort!(collect(keys(arr.data)))
    return [arr[i] for i in ks]
end

# Double argument singleton
function vectorize(arr::T, inds::Nothing)::T where {T}
    return vectorize(arr)
end

# Double argument Vector/Array/DenseAxisArray
function vectorize(
    arr::AbstractArray{T},
    inds::ContainerIndices
    )::Vector{T} where {T}
    return vectorize(arr)
end

# Double argument SparseAxisArray
function vectorize(
    arr::JuMPC.SparseAxisArray{T},
    inds::ContainerIndices
    )::Vector{T} where {T}
    return [arr[i] for i in inds.axes]
end

"""
    unvectorize(vect, inds::Union{Nothing, ContainerIndices})

Return the unvectorized version of `vect` as dictated by `inds`. Here `vect`
comes from a call of [`vectorize`](@ref InfiniteOpt.Collections.vectorize) and
`inds` comes from a call of [`indices`](@ref InfiniteOpt.Collections.indices).
"""
function unvectorize(vect::T, inds::Nothing)::T where {T}
    return vect
end

# Vector
function unvectorize(
    vect::Vector{T},
    inds::ContainerIndices{1, Nothing}
    )::Vector{T} where {T}
    return vect
end

# Array
function unvectorize(
    vect::Vector{T},
    inds::ContainerIndices{N, Nothing}
    )::Array{T, N} where {T, N}
    return [vect[i] for i in LinearIndices(inds.indices)]
end

# DenseAxisArray (vector)
function unvectorize(vect::Vector, inds::ContainerIndices{1})
    return JuMPC.DenseAxisArray(vect, inds.axes...)
end

# DenseAxisArray (general)
function unvectorize(vect::Vector, inds::ContainerIndices)
    arr = [vect[i] for i in LinearIndices(inds.indices)]
    return JuMPC.DenseAxisArray(arr, inds.axes...)
end

# SparseAxisArray
function unvectorize(
    vect::Vector{T},
    inds::ContainerIndices{1, Vector{K}}
    )::JuMPC.SparseAxisArray{T, N, K} where {T, N, K <: NTuple{N}}
    dict = Dict(inds.axes[i] => vect[i] for i in eachindex(inds.axes))
    return JuMPC.SparseAxisArray(dict)
end


text:
using NiLang, NiLang.AD
using Test, Random
using ForwardDiff
using NiLangCore: default_constructor

@testset "basic instructions, ULogarithmic" begin
    x, y = default_constructor(ULogarithmic{Int}, 1),
    default_constructor(ULogarithmic{Int}, 2)
	@instr x *= y
    @test x == default_constructor(ULogarithmic{Int}, 3)
    @test y == default_constructor(ULogarithmic{Int}, 2)

	@test PlusEq(gaussian_log)(1.0, 2.0) == (1.0+log(1+exp(2.0)), 2.0)

    x, y, z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 2.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= y + z
	@test check_inv(MulEq(+), (x, y, z))
	@test x.log ≈ log(exp(7.0) * (exp(2.0) + exp(3.0)))
    x, y, z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= y - z
	@test x.log ≈ log(exp(7.0) * (exp(5.0) - exp(3.0)))

    x, y, z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= y^3.4
	@test x.log ≈ log(exp(5.0)^3.4 * exp(7.0))

    x, y, z = default_constructor(ULogarithmic{Float64}, 7.0),
    default_constructor(ULogarithmic{Float64}, 5.0),
    default_constructor(ULogarithmic{Float64}, 3.0)
	@instr x *= 3
	@test x.log ≈ log(exp(7.0) * 3)
end

@testset "error on += and -=" begin
    @i function f(x::ULogarithmic)
        x += ULogarithmic(3.0)
    end

    @test_throws MethodError f(ULogarithmic(2.0))
    @test_throws MethodError (~f)(ULogarithmic(2.0))
end


text:
using Base: PermutedDimsArrays

# PermutedDimsArrays' custom copyto! doesn't know how to deal with GPU arrays
function PermutedDimsArrays._copy!(dest::PermutedDimsArray{T,N,<:Any,<:Any,<:GPUArray}, src) where {T,N}
    dest .= src
    dest
end


text:
function part1(input)
    x = collect(readchomp(input))
    return count(==('('), x) - count(==(')'), x)
end

function part2(input)
    x = collect(readchomp(input))
    level = 0
    for i = 1:length(x)
        level += x[i] == '(' ? 1 : -1
        level < 0 && return i
    end
end


text:
# generic functions fot distribution fitting

function suffstats{D<:Distribution}(dt::Type{D}, xs...)
    argtypes = tuple(D, map(typeof, xs)...)
    error("suffstats is not implemented for $argtypes.")
end

fit_mle{D<:UnivariateDistribution}(dt::Type{D}, x::Array) = fit_mle(D, suffstats(D, x))
fit_mle{D<:UnivariateDistribution}(dt::Type{D}, x::Array, w::Array) = fit_mle(D, suffstats(D, x, w))

fit_mle{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix) = fit_mle(D, suffstats(D, x))
fit_mle{D<:MultivariateDistribution}(dt::Type{D}, x::Matrix, w::Array) = fit_mle(D, suffstats(D, x, w))

fit{D<:Distribution}(dt::Type{D}, x) = fit_mle(D, x)
fit{D<:Distribution}(dt::Type{D}, args...) = fit_mle(D, args...)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DdosSettings <: SwaggerModel
    ddosCustomPolicy::Any # spec type: Union{ Nothing, SubResource } # spec name: ddosCustomPolicy
    protectionCoverage::Any # spec type: Union{ Nothing, String } # spec name: protectionCoverage
    protectedIP::Any # spec type: Union{ Nothing, Bool } # spec name: protectedIP

    function DdosSettings(;ddosCustomPolicy=nothing, protectionCoverage=nothing, protectedIP=nothing)
        o = new()
        validate_property(DdosSettings, Symbol("ddosCustomPolicy"), ddosCustomPolicy)
        setfield!(o, Symbol("ddosCustomPolicy"), ddosCustomPolicy)
        validate_property(DdosSettings, Symbol("protectionCoverage"), protectionCoverage)
        setfield!(o, Symbol("protectionCoverage"), protectionCoverage)
        validate_property(DdosSettings, Symbol("protectedIP"), protectedIP)
        setfield!(o, Symbol("protectedIP"), protectedIP)
        o
    end
end # type DdosSettings

const _property_map_DdosSettings = Dict{Symbol,Symbol}(Symbol("ddosCustomPolicy")=>Symbol("ddosCustomPolicy"), Symbol("protectionCoverage")=>Symbol("protectionCoverage"), Symbol("protectedIP")=>Symbol("protectedIP"))
const _property_types_DdosSettings = Dict{Symbol,String}(Symbol("ddosCustomPolicy")=>"SubResource", Symbol("protectionCoverage")=>"String", Symbol("protectedIP")=>"Bool")
Base.propertynames(::Type{ DdosSettings }) = collect(keys(_property_map_DdosSettings))
Swagger.property_type(::Type{ DdosSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DdosSettings[name]))}
Swagger.field_name(::Type{ DdosSettings }, property_name::Symbol) =  _property_map_DdosSettings[property_name]

const _allowed_DdosSettings_protectionCoverage = ["Basic", "Standard"]

function check_required(o::DdosSettings)
    true
end

function validate_property(::Type{ DdosSettings }, name::Symbol, val)
    if name === Symbol("protectionCoverage")
        Swagger.validate_param(name, "DdosSettings", :enum, val, _allowed_DdosSettings_protectionCoverage)
    end
end


text:
import JuLIP
using FunctionWrappers,  BenchmarkTools, StaticPolynomials, StaticArrays

import NBodyIPs
using DynamicPolynomials: @polyvar
using JuLIP.Potentials: coscut, F64fun

##
@info "[1] Profile FunctionWrappers"
println("    1000 evaluations of coscut")
f(r) = coscut(r, 1.123, 3.456)
fwrap = F64fun(f)
r = 1.0 .+ 3*rand(1000)
test1(f, r) = (s=0.0; for n = 1:length(r); s += f(r[n]); end; s)
print(" Plain Function: "); @btime test1($f, $r)
print("FunctionWrapper: "); @btime test1($fwrap, $r)

@info "[2] Profile StaticPolynomials"

##
@info "[2.1] From the StaticPolynomials README"

function f10(x)
    f  = 48*x[1]^3 + 72*x[1]^2*x[2] + 72*x[1]^2*x[3] + 72*x[1]^2*x[4] + 72*x[1]^2*x[5] + 72*x[1]^2*x[7]
    f += 72*x[1]^2*x[8] + 72*x[1]*x[2]^2 + 144*x[1]*x[2]*x[4] + 144*x[1]*x[2]*x[7] + 72*x[1]*x[3]^2
    f += 144*x[1]*x[3]*x[5] + 144*x[1]*x[3]*x[8] + 72*x[1]*x[4]^2 + 144*x[1]*x[4]*x[7] + 72*x[1]*x[5]^2
    f += 144*x[1]*x[5]*x[8] + 72*x[1]*x[7]^2 + 72*x[1]*x[8]^2 + 48*x[2]^3 + 72*x[2]^2*x[3]
    f += 72*x[2]^2*x[4] + 72*x[2]^2*x[6] + 72*x[2]^2*x[7] + 72*x[2]^2*x[9] + 72*x[2]*x[3]^2
    f += 144*x[2]*x[3]*x[6] + 144*x[2]*x[3]*x[9] + 72*x[2]*x[4]^2 + 144*x[2]*x[4]*x[7] + 72*x[2]*x[6]^2
    f += 144*x[2]*x[6]*x[9] + 72*x[2]*x[7]^2 + 72*x[2]*x[9]^2 + 48*x[3]^3 + 72*x[3]^2*x[5]
    f += 72*x[3]^2*x[6] + 72*x[3]^2*x[8] + 72*x[3]^2*x[9] + 72*x[3]*x[5]^2 + 144*x[3]*x[5]*x[8]
    f += 72*x[3]*x[6]^2 + 144*x[3]*x[6]*x[9] + 72*x[3]*x[8]^2 + 72*x[3]*x[9]^2 + 48*x[4]^3
    f += 72*x[4]^2*x[5] + 72*x[4]^2*x[6] + 72*x[4]^2*x[7] + 72*x[4]^2*x[10] + 72*x[4]*x[5]^2
    f += 144*x[4]*x[5]*x[6] + 144*x[4]*x[5]*x[10] + 72*x[4]*x[6]^2 + 144*x[4]*x[6]*x[10] + 72*x[4]*x[7]^2
    f += 72*x[4]*x[10]^2 + 48*x[5]^3 + 72*x[5]^2*x[6] + 72*x[5]^2*x[8] + 72*x[5]^2*x[10]
    f += 72*x[5]*x[6]^2 + 144*x[5]*x[6]*x[10] + 72*x[5]*x[8]^2 + 72*x[5]*x[10]^2 + 48*x[6]^3
    f += 72*x[6]^2*x[9] + 72*x[6]^2*x[10] + 72*x[6]*x[9]^2 + 72*x[6]*x[10]^2 + 48*x[7]^3
    f += 72*x[7]^2*x[8] + 72*x[7]^2*x[9] + 72*x[7]^2*x[10] + 72*x[7]*x[8]^2 + 144*x[7]*x[8]*x[9]
    f += 144*x[7]*x[8]*x[10] + 72*x[7]*x[9]^2 + 144*x[7]*x[9]*x[10] + 72*x[7]*x[10]^2 + 48*x[8]^3
    f += 72*x[8]^2*x[9] + 72*x[8]^2*x[10] + 72*x[8]*x[9]^2 + 144*x[8]*x[9]*x[10] + 72*x[8]*x[10]^2
    f += 48*x[9]^3 + 72*x[9]^2*x[10] + 72*x[9]*x[10]^2 + 48*x[10]^3
    return f
end

@polyvar x[1:10]
p10 = Polynomial(f10(x))
y = @SVector rand(10)
test2(f, x) = (s=0.0; for n=1:1000; s+= f(x); end; s)
print("     Manual Implementation: "); @btime test2($f10, $y)
print("         StaticPolynomials: "); @btime test2($p10, $y)
print("StaticPolynomials.gradient: "); @btime StaticPolynomials.gradient($p10, $y)

##
@info "[2.2] A BL5 Invariant"
f10I(x) = (x[1]^2*x[2]*x[5]*x[3] + x[1]*x[2]^2*x[5]*x[3] + x[1]^2*x[2]*x[5]*x[6] + x[1]*x[2]*x[5]^2*x[6] + x[1]^2*x[2]*x[3]*x[6] + x[1]^2*x[5]*x[3]*x[6] + x[1]*x[2]*x[3]^2*x[6] + x[1]*x[5]*x[3]*x[6]^2 + x[1]*x[2]^2*x[5]*x[8] + x[1]*x[2]*x[5]^2*x[8] + x[1]*x[2]^2*x[3]*x[8] + x[2]^2*x[5]*x[3]*x[8] + x[1]*x[2]*x[3]^2*x[8] + x[1]*x[5]^2*x[6]*x[8] + x[2]*x[5]^2*x[6]*x[8] + x[1]*x[3]^2*x[6]*x[8] + x[2]*x[3]^2*x[6]*x[8] + x[1]*x[5]*x[6]^2*x[8] + x[1]*x[3]*x[6]^2*x[8] + x[5]*x[3]*x[6]^2*x[8] + x[2]*x[5]*x[3]*x[8]^2 + x[2]*x[5]*x[6]*x[8]^2 + x[2]*x[3]*x[6]*x[8]^2 + x[5]*x[3]*x[6]*x[8]^2 + x[1]^2*x[2]*x[5]*x[4] + x[1]*x[2]^2*x[5]*x[4] + x[1]^2*x[3]*x[6]*x[4] + x[1]*x[3]^2*x[6]*x[4] + x[2]^2*x[3]*x[8]*x[4] + x[2]*x[3]^2*x[8]*x[4] + x[1]^2*x[2]*x[5]*x[7] + x[1]*x[2]*x[5]^2*x[7] + x[1]^2*x[3]*x[6]*x[7] + x[1]*x[3]*x[6]^2*x[7] + x[5]^2*x[6]*x[8]*x[7] + x[5]*x[6]^2*x[8]*x[7] + x[1]^2*x[2]*x[4]*x[7] + x[1]^2*x[5]*x[4]*x[7] + x[1]^2*x[3]*x[4]*x[7] + x[1]^2*x[6]*x[4]*x[7] + x[1]*x[2]*x[4]^2*x[7] + x[1]*x[3]*x[4]^2*x[7] + x[1]*x[5]*x[4]*x[7]^2 + x[1]*x[6]*x[4]*x[7]^2 + x[1]*x[2]^2*x[5]*x[9] + x[1]*x[2]*x[5]^2*x[9] + x[2]^2*x[3]*x[8]*x[9] + x[5]^2*x[6]*x[8]*x[9] + x[2]*x[3]*x[8]^2*x[9] + x[5]*x[6]*x[8]^2*x[9] + x[1]*x[2]^2*x[4]*x[9] + x[2]^2*x[5]*x[4]*x[9] + x[2]^2*x[3]*x[4]*x[9] + x[2]^2*x[8]*x[4]*x[9] + x[1]*x[2]*x[4]^2*x[9] + x[2]*x[3]*x[4]^2*x[9] + x[1]*x[5]^2*x[7]*x[9] + x[2]*x[5]^2*x[7]*x[9] + x[5]^2*x[6]*x[7]*x[9] + x[5]^2*x[8]*x[7]*x[9] + x[1]*x[4]^2*x[7]*x[9] + x[2]*x[4]^2*x[7]*x[9] + x[1]*x[5]*x[7]^2*x[9] + x[5]*x[6]*x[7]^2*x[9] + x[1]*x[4]*x[7]^2*x[9] + x[5]*x[4]*x[7]^2*x[9] + x[2]*x[5]*x[4]*x[9]^2 + x[2]*x[8]*x[4]*x[9]^2 + x[2]*x[5]*x[7]*x[9]^2 + x[5]*x[8]*x[7]*x[9]^2 + x[2]*x[4]*x[7]*x[9]^2 + x[5]*x[4]*x[7]*x[9]^2 + x[1]*x[3]^2*x[6]*x[10] + x[1]*x[3]*x[6]^2*x[10] + x[2]*x[3]^2*x[8]*x[10] + x[5]*x[6]^2*x[8]*x[10] + x[2]*x[3]*x[8]^2*x[10] + x[5]*x[6]*x[8]^2*x[10] + x[1]*x[3]^2*x[4]*x[10] + x[2]*x[3]^2*x[4]*x[10] + x[3]^2*x[6]*x[4]*x[10] + x[3]^2*x[8]*x[4]*x[10] + x[1]*x[3]*x[4]^2*x[10] + x[2]*x[3]*x[4]^2*x[10] + x[1]*x[6]^2*x[7]*x[10] + x[5]*x[6]^2*x[7]*x[10] + x[3]*x[6]^2*x[7]*x[10] + x[6]^2*x[8]*x[7]*x[10] + x[1]*x[4]^2*x[7]*x[10] + x[3]*x[4]^2*x[7]*x[10] + x[1]*x[6]*x[7]^2*x[10] + x[5]*x[6]*x[7]^2*x[10] + x[1]*x[4]*x[7]^2*x[10] + x[6]*x[4]*x[7]^2*x[10] + x[2]*x[8]^2*x[9]*x[10] + x[5]*x[8]^2*x[9]*x[10] + x[3]*x[8]^2*x[9]*x[10] + x[6]*x[8]^2*x[9]*x[10] + x[2]*x[4]^2*x[9]*x[10] + x[3]*x[4]^2*x[9]*x[10] + x[5]*x[7]^2*x[9]*x[10] + x[6]*x[7]^2*x[9]*x[10] + x[2]*x[8]*x[9]^2*x[10] + x[5]*x[8]*x[9]^2*x[10] + x[2]*x[4]*x[9]^2*x[10] + x[8]*x[4]*x[9]^2*x[10] + x[5]*x[7]*x[9]^2*x[10] + x[8]*x[7]*x[9]^2*x[10] + x[3]*x[6]*x[4]*x[10]^2 + x[3]*x[8]*x[4]*x[10]^2 + x[3]*x[6]*x[7]*x[10]^2 + x[6]*x[8]*x[7]*x[10]^2 + x[3]*x[4]*x[7]*x[10]^2 + x[6]*x[4]*x[7]*x[10]^2 + x[3]*x[8]*x[9]*x[10]^2 + x[6]*x[8]*x[9]*x[10]^2 + x[3]*x[4]*x[9]*x[10]^2 + x[8]*x[4]*x[9]*x[10]^2 + x[6]*x[7]*x[9]*x[10]^2 + x[8]*x[7]*x[9]*x[10]^2)
p10I = Polynomial(f10I(x))
# print("     Manual Implementation: "); @btime test2($f10I, $y)
print("         StaticPolynomials: "); @btime test2($p10I, $y)
print("StaticPolynomials.gradient: "); @btime StaticPolynomials.gradient($p10I, $y)

##
@info "[2.3] A generic 4B total-degree 10 polynomial"
function sp_exps(A)
   exps = zeros(Int, (length(A), 12))
   for (i, a) in enumerate(A)
      exps[i,1:6] .= a[1:6]
      exps[i,7+a[end]] = 1
   end
   return exps
end
desc = NBodyIPs.BondLengthDesc("r->r", (:cos, 1.0, 2.0))
A = NBodyIPs.PolyBasis.gen_tuples(desc, 4, 10)
C = rand(Float64, length(A))
p4B = StaticPolynomials.Polynomial(C, collect(sp_exps(A)'))
x23 = @SVector rand(12)
print("         StaticPolynomials: "); @btime $p4B($x23)
print("StaticPolynomials.gradient: "); @btime StaticPolynomials.gradient($p4B, $x23)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# test format spec parsing

# default spec
@testset "Default spec" begin
    fs = FormatSpec("")
    @test fs.typ == 's'
    @test fs.fill == ' '
    @test fs.align == '<'
    @test fs.sign == '-'
    @test fs.width == -1
    @test fs.prec == -1
    @test fs.ipre == false
    @test fs.zpad == false
    @test fs.tsep == false
end

_contains(s, r) = @static VERSION < v"0.7.0-DEV" ? contains(s, r) : occursin(r, s)

@testset "Show" begin
    x = FormatSpec("#8,d")
    io = IOBuffer()
    show(io, x)
    str = String(take!(io))
    @test _contains(str, "width = 8")
end

@testset "Literal incorrect" begin
    @test_throws ErrorException FormatSpec("Z")
end

# more cases

@testset "FormatSpec(\"d\")" begin
    fs = FormatSpec("d")
    @test fs == FormatSpec('d')
    @test fs.align == '>'

    @test FormatSpec("8x") == FormatSpec('x'; width=8)
    @test FormatSpec("08b") == FormatSpec('b'; width=8, zpad=true)
    @test FormatSpec("12f") == FormatSpec('f'; width=12, prec=6)
    @test FormatSpec("12.7f") == FormatSpec('f'; width=12, prec=7)
    @test FormatSpec("+08o") == FormatSpec('o'; width=8, zpad=true, sign='+')

    @test FormatSpec("8") == FormatSpec('s'; width=8)
    @test FormatSpec(".6f") == FormatSpec('f'; prec=6)
    @test FormatSpec("<8d") == FormatSpec('d'; width=8, align='<')
    @test FormatSpec("#<8d") == FormatSpec('d'; width=8, fill='#', align='<')
    @test FormatSpec("⋆<8d") == FormatSpec('d'; width=8, fill='⋆', align='<')
    @test FormatSpec("#8,d") == FormatSpec('d'; width=8, ipre=true, tsep=true)
end

@testset "Format prefix" begin
    @test pyfmt("#b", 6) == "0b110"
    @test pyfmt("#o", 6) == "0o6"
    @test pyfmt("#x", 6) == "0x6"
end

@testset "Format string" begin
    @test pyfmt("", "abc") == "abc"
    @test pyfmt("", "αβγ") == "αβγ"
    @test pyfmt("s", "abc") == "abc"
    @test pyfmt("s", "αβγ") == "αβγ"
    @test pyfmt("2s", "abc") == "abc"
    @test pyfmt("2s", "αβγ") == "αβγ"
    @test pyfmt("5s", "abc") == "abc  "
    @test pyfmt("5s", "αβγ") == "αβγ  "
    @test pyfmt(">5s", "abc") == "  abc"
    @test pyfmt(">5s", "αβγ") == "  αβγ"
    @test pyfmt("*>5s", "abc") == "**abc"
    @test pyfmt("⋆>5s", "αβγ") == "⋆⋆αβγ"
    @test pyfmt("*<5s", "abc") == "abc**"
    @test pyfmt("⋆<5s", "αβγ") == "αβγ⋆⋆"
end

@testset "Format Char" begin
    @test pyfmt("", 'c') == "c"
    @test pyfmt("", 'γ') == "γ"
    @test pyfmt("c", 'c') == "c"
    @test pyfmt("c", 'γ') == "γ"
    @test pyfmt("3c", 'c') == "c  "
    @test pyfmt("3c", 'γ') == "γ  "
    @test pyfmt(">3c", 'c') == "  c"
    @test pyfmt(">3c", 'γ') == "  γ"
    @test pyfmt("*>3c", 'c') == "**c"
    @test pyfmt("⋆>3c", 'γ') == "⋆⋆γ"
    @test pyfmt("*<3c", 'c') == "c**"
    @test pyfmt("⋆<3c", 'γ') == "γ⋆⋆"
end

@testset "Format integer" begin

    @test pyfmt("", 1234) == "1234"
    @test pyfmt("d", 1234) == "1234"
    @test pyfmt("n", 1234) == "1234"
    @test pyfmt("x", 0x2ab) == "2ab"
    @test pyfmt("X", 0x2ab) == "2AB"
    @test pyfmt("o", 0o123) == "123"
    @test pyfmt("b", 0b1101) == "1101"

    @test pyfmt("d", 0) == "0"
    @test pyfmt("d", 9) == "9"
    @test pyfmt("d", 10) == "10"
    @test pyfmt("d", 99) == "99"
    @test pyfmt("d", 100) == "100"
    @test pyfmt("d", 1000) == "1000"

    @test pyfmt("06d", 123) == "000123"
    @test pyfmt("+6d", 123) == "  +123"
    @test pyfmt("+06d", 123) == "+00123"
    @test pyfmt(" d", 123) == " 123"
    @test pyfmt(" 6d", 123) == "   123"
    @test pyfmt("<6d", 123) == "123   "
    @test pyfmt(">6d", 123) == "   123"
    @test pyfmt("*<6d", 123) == "123***"
    @test pyfmt("⋆<6d", 123) == "123⋆⋆⋆"
    @test pyfmt("*>6d", 123) == "***123"
    @test pyfmt("⋆>6d", 123) == "⋆⋆⋆123"
    @test pyfmt("< 6d", 123) == " 123  "
    @test pyfmt("<+6d", 123) == "+123  "
    @test pyfmt("> 6d", 123) == "   123"
    @test pyfmt(">+6d", 123) == "  +123"

    @test pyfmt("+d", -123) == "-123"
    @test pyfmt("-d", -123) == "-123"
    @test pyfmt(" d", -123) == "-123"
    @test pyfmt("06d", -123) == "-00123"
    @test pyfmt("<6d", -123) == "-123  "
    @test pyfmt(">6d", -123) == "  -123"
end

@testset "Format floating point (f)" begin

    @test pyfmt("", 0.125) == "0.125"
    @test pyfmt("f", 0.0) == "0.000000"
    @test pyfmt("f", -0.0) == "-0.000000"
    @test pyfmt("f", 0.001) == "0.001000"
    @test pyfmt("f", 0.125) == "0.125000"
    @test pyfmt("f", 1.0/3) == "0.333333"
    @test pyfmt("f", 1.0/6) == "0.166667"
    @test pyfmt("f", -0.125) == "-0.125000"
    @test pyfmt("f", -1.0/3) == "-0.333333"
    @test pyfmt("f", -1.0/6) == "-0.166667"
    @test pyfmt("f", 1234.5678) == "1234.567800"
    @test pyfmt("8f", 1234.5678) == "1234.567800"

    @test pyfmt("8.2f", 8.376) == "    8.38"
    @test pyfmt("<8.2f", 8.376) == "8.38    "
    @test pyfmt(">8.2f", 8.376) == "    8.38"
    @test pyfmt("8.2f", -8.376) == "   -8.38"
    @test pyfmt("<8.2f", -8.376) == "-8.38   "
    @test pyfmt(">8.2f", -8.376) == "   -8.38"

    @test pyfmt("<08.2f", 8.376) == "00008.38"
    @test pyfmt(">08.2f", 8.376) == "00008.38"
    @test pyfmt("<08.2f", -8.376) == "-0008.38"
    @test pyfmt(">08.2f", -8.376) == "-0008.38"
    @test pyfmt("*<8.2f", 8.376) == "8.38****"
    @test pyfmt("⋆<8.2f", 8.376) == "8.38⋆⋆⋆⋆"
    @test pyfmt("*>8.2f", 8.376) == "****8.38"
    @test pyfmt("⋆>8.2f", 8.376) == "⋆⋆⋆⋆8.38"
    @test pyfmt("*<8.2f", -8.376) == "-8.38***"
    @test pyfmt("⋆<8.2f", -8.376) == "-8.38⋆⋆⋆"
    @test pyfmt("*>8.2f", -8.376) == "***-8.38"
    @test pyfmt("⋆>8.2f", -8.376) == "⋆⋆⋆-8.38"

    @test pyfmt(".2f", 0.999) == "1.00"
    @test pyfmt(".2f", 0.996) == "1.00"
    @test pyfmt("6.2f", 9.999) == " 10.00"
    # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
    @test (pyfmt(".2f", 0.995) == "1.00" || pyfmt(".2f", 0.995) == "0.99")
    @test pyfmt(".2f", 0.994) == "0.99"
end

@testset "Format floating point (e)" begin

    @test pyfmt("E", 0.0) == "0.000000E+00"
    @test pyfmt("e", 0.0) == "0.000000e+00"
    @test pyfmt("e", 0.001) == "1.000000e-03"
    @test pyfmt("e", 0.125) == "1.250000e-01"
    @test pyfmt("e", 100/3) == "3.333333e+01"
    @test pyfmt("e", -0.125) == "-1.250000e-01"
    @test pyfmt("e", -100/6) == "-1.666667e+01"
    @test pyfmt("e", 1234.5678) == "1.234568e+03"
    @test pyfmt("8e", 1234.5678) == "1.234568e+03"

    @test pyfmt("<12.2e", 13.89) == "1.39e+01    "
    @test pyfmt(">12.2e", 13.89) == "    1.39e+01"
    @test pyfmt("*<12.2e", 13.89) == "1.39e+01****"
    @test pyfmt("⋆<12.2e", 13.89) == "1.39e+01⋆⋆⋆⋆"
    @test pyfmt("*>12.2e", 13.89) == "****1.39e+01"
    @test pyfmt("⋆>12.2e", 13.89) == "⋆⋆⋆⋆1.39e+01"
    @test pyfmt("012.2e", 13.89) == "00001.39e+01"
    @test pyfmt("012.2e", -13.89) == "-0001.39e+01"
    @test pyfmt("+012.2e", 13.89) == "+0001.39e+01"

    @test pyfmt(".1e", 0.999) == "1.0e+00"
    @test pyfmt(".1e", 0.996) == "1.0e+00"
    # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
    @test (pyfmt(".1e", 0.995) == "1.0e+00" || pyfmt(".1e", 0.995) == "9.9e-01")
    @test pyfmt(".1e", 0.994) == "9.9e-01"
    @test pyfmt(".1e", 0.6) == "6.0e-01"
    @test pyfmt(".1e", 0.9) == "9.0e-01"

    @test pyfmt("10.2e", 1.2e100) == " 1.20e+100"
    @test pyfmt("11.2e", BigFloat("1.2e1000")) == " 1.20e+1000"
    @test pyfmt("11.2e", BigFloat("1.2e-1000")) == " 1.20e-1000"
    @test pyfmt("9.2e", 9.999e9) == " 1.00e+10"
    @test pyfmt("10.2e", 9.999e99) == " 1.00e+100"
    @test pyfmt("11.2e", BigFloat("9.999e999")) == " 1.00e+1000"
    @test pyfmt("10.2e", -9.999e-100) == " -1.00e-99"
end

@testset "Format special floating point value" begin

    @test pyfmt("f", NaN) == "NaN"
    @test pyfmt("e", NaN) == "NaN"
    @test pyfmt("f", NaN32) == "NaN"
    @test pyfmt("e", NaN32) == "NaN"

    @test pyfmt("f", Inf) == "Inf"
    @test pyfmt("e", Inf) == "Inf"
    @test pyfmt("f", Inf32) == "Inf"
    @test pyfmt("e", Inf32) == "Inf"

    @test pyfmt("f", -Inf) == "-Inf"
    @test pyfmt("e", -Inf) == "-Inf"
    @test pyfmt("f", -Inf32) == "-Inf"
    @test pyfmt("e", -Inf32) == "-Inf"

    @test pyfmt("<5f", Inf) == "Inf  "
    @test pyfmt(">5f", Inf) == "  Inf"
    @test pyfmt("*<5f", Inf) == "Inf**"
    @test pyfmt("⋆<5f", Inf) == "Inf⋆⋆"
    @test pyfmt("*>5f", Inf) == "**Inf"
    @test pyfmt("⋆>5f", Inf) == "⋆⋆Inf"
end

@testset "Format Symbols (S) for Irrationals" begin
    @test pyfmt(">10s", pi) == "         π"
    @test pyfmt("10s", pi) == "π         "
    @test pyfmt("3", MathConstants.eulergamma) == "γ  "
    @test pyfmt("10.2f", MathConstants.eulergamma) == "      0.58"
    @test pyfmt("<3s", MathConstants.e) == "ℯ  "
end

@testset "Format Symbols (S) for Irrationals" begin
    pyfmt("10", 3//4) == "3//4      "
    pyfmt(">10", 3//4) == "      3//4"
    pyfmt("10.1f", 3//4) == "       0.8"
end


text:
module MarkovAlgos

struct MarkovRule{F,T}
    patt::F
    repl::T
    term::Bool
end

isterminating(r::MarkovRule) = r.term
Base.show(io::IO, rule::MarkovRule) =
    print(io, rule.patt, " → ", isterminating(rule) ? "." : "", rule.repl)
function Base.convert(::Type{MarkovRule}, s::AbstractString)
    rmatch = match(r"^(.+)\s+->\s*(\.)?(.*)?$", s)
    if rmatch ≡ nothing || isempty(rmatch.captures)
        throw(ParseError("not valid rule: " * s))
    end
    patt, term, repl = rmatch.captures
    return MarkovRule(patt, repl ≢ nothing ? repl : "", term ≢ nothing)
end

function ruleset(file::Union{AbstractString,IO})
    ruleset = Vector{MarkovRule}(0)
    for line in eachline(file)
        ismatch(r"(^#|^\s*$)", line) || push!(ruleset, MarkovRule(line))
    end
    return ruleset
end

apply(text::AbstractString, rule::MarkovRule) = replace(text, rule.patt, rule.repl)
function apply(file::Union{AbstractString,IO}, ruleset::AbstractVector{<:MarkovRule})
    text = readstring(file)
    redo = !isempty(text)
    while redo
        matchrule = false
        for rule in ruleset
            if contains(text, rule.patt)
                matchrule = true
                text = apply(text, rule)
                redo = !isterminating(rule)
                break
            end
        end
        redo = redo && matchrule
    end
    return text
end

end  # module MarkovAlgos

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

"""
    PointProcess

A spatial point process.
"""
abstract type PointProcess end

"""
    ishomogeneous(p)

Tells whether or not the spatial point process `p` is homogeneous.
"""
ishomogeneous(p::PointProcess) = false

"""
    rand([rng], p, g, n=1; [algo])

Generate `n` realizations of spatial point process `p`
inside geometry `g`. Optionally specify sampling
algorithm `algo` and random number generator `rng`.
"""
Base.rand(rng::Random.AbstractRNG,
          p::PointProcess, g::Geometry, n::Int;
          algo=default_sampling_algorithm(p, g)) =
  [rand_single(rng, p, g, algo) for i in 1:n]

Base.rand(rng::Random.AbstractRNG,
          p::PointProcess, g::Geometry;
          algo=default_sampling_algorithm(p, g)) =
  rand_single(rng, p, g, algo)

Base.rand(p::PointProcess, g::Geometry, n::Int;
          algo=default_sampling_algorithm(p, g)) =
  rand(Random.GLOBAL_RNG, p, g, n; algo=algo)

Base.rand(p::PointProcess, g::Geometry;
          algo=default_sampling_algorithm(p, g)) =
  rand(Random.GLOBAL_RNG, p, g; algo=algo)

"""
    rand_single(rng, p, g, algo)

Generate a single realization of spatial point process
`p` inside geometry `g` with sampling `algo`.
"""
function rand_single end

"""
    default_sampling_algorithm(p, g)

Default sampling algorithm for spatial point process `p`
on geometry `g`.
"""
function default_sampling_algorithm end

# --------------------
# SAMPLING ALGORITHMS
# --------------------

struct ProductSampling end
struct ThinnedSampling end
struct DiscretizedSampling end
struct UnionSampling end

#-----------------
# IMPLEMENTATIONS
#-----------------

include("processes/binomial.jl")
include("processes/poisson.jl")
include("processes/union.jl")


text:
"""
    max_directory_ctime(prefix::String)

Takes the `stat()` of all files in a directory root, keeping the maximum ctime,
recursively.  Comparing just this value allows for quick directory change detection.
"""
function max_directory_ctime(prefix::String)
    max_time = 0.0
    for (root, dirs, files) in walkdir(prefix)
        for f in files
            max_time = max(max_time, lstat(joinpath(root, f)).ctime)
        end
    end
    return max_time
end

"""
    is_ecryptfs(path::AbstractString; verbose::Bool=false)

Checks to see if the given `path` (or any parent directory) is placed upon an
`ecryptfs` mount.  This is known not to work on current kernels, see this bug
for more details: https://bugzilla.kernel.org/show_bug.cgi?id=197603

This method returns whether it is encrypted or not, and what mountpoint it used
to make that decision.
"""
function is_ecryptfs(path::AbstractString; verbose::Bool=false)
    # Canonicalize `path` immediately, and if it's a directory, add a "/" so
    # as to be consistent with the rest of this function
    path = abspath(path)
    if isdir(path)
        path = abspath(path * "/")
    end

    if verbose
        @info("Checking to see if $path is encrypted...")
    end

    # Get a listing of the current mounts.  If we can't do this, just give up
    if !isfile("/proc/mounts")
        if verbose
            @info("Couldn't open /proc/mounts, returning...")
        end
        return false, path
    end
    mounts = String(read("/proc/mounts"))

    # Grab the fstype and the mountpoints
    mounts = [split(m)[2:3] for m in split(mounts, "\n") if !isempty(m)]

    # Canonicalize mountpoints now so as to dodge symlink difficulties
    mounts = [(abspath(m[1]*"/"), m[2]) for m in mounts]

    # Fast-path asking for a mountpoint directly (e.g. not a subdirectory)
    direct_path = [m[1] == path for m in mounts]
    local parent
    if any(direct_path)
        parent = mounts[findfirst(direct_path)]
    else
        # Find the longest prefix mount:
        parent_mounts = [m for m in mounts if startswith(path, m[1])]
        if isempty(parent_mounts)
            # This is weird; this means that we can't find any mountpoints that
            # hold the given path.  I've only ever seen this in `chroot`'ed scenarios.
            return false, path
        end
        parent = parent_mounts[argmax(map(m->length(m[1]), parent_mounts))]
    end

    # Return true if this mountpoint is an ecryptfs mount
    val = parent[2] == "ecryptfs"
    if verbose && val
        @info("  -> $path is encrypted from mountpoint $(parent[1])")
    end
    return val, parent[1]
end

"""
    uname()

On Linux systems, return the strings returned by the `uname()` function in libc.
"""
function uname()
    if !Sys.islinux()
        return String[]
    end

    # Get libc and handle to uname
    libcs = filter(x -> occursin("libc.so", x), dllist())
    if isempty(libcs)
        error("Could not find libc, unable to call uname()")
    end
    libc = dlopen(first(libcs))
    uname_hdl = dlsym(libc, :uname)

    # The uname struct can have wildly differing layouts; we take advantage
    # of the fact that it is just a bunch of NULL-terminated strings laid out
    # one after the other, and that it is (as best as I can tell) at maximum
    # around 1.5KB long.  We bump up to 2KB to be safe.
    uname_struct = zeros(UInt8, 2048)
    ccall(uname_hdl, Cint, (Ptr{UInt8},), uname_struct)

    # Parse out all the strings embedded within this struct
    strings = String[]
    idx = 1
    while idx < length(uname_struct)
        # Extract string
        new_string = unsafe_string(pointer(uname_struct, idx))
        push!(strings, new_string)
        idx += length(new_string) + 1

        # Skip trailing zeros
        while uname_struct[idx] == 0 && idx < length(uname_struct)
            idx += 1
        end
    end

    return strings
end

"""
    get_kernel_version(;verbose::Bool = false)

Use `uname()` to get the kernel version and parse it out as a `VersionNumber`,
returning `nothing` if parsing fails or this is not `Linux`.
"""
function get_kernel_version(;verbose::Bool = false)
    if !Sys.islinux()
        return nothing
    end

    uname_strings = try
        uname()
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end

        @warn("Unable to run `uname()` to check version number!")
        return nothing
    end

    # Some distributions tack extra stuff onto the version number.  We walk backwards
    # from the end, searching for the longest string that we can extract a VersionNumber
    # out of.  We choose a minimum length of 5, as all kernel version numbers will be at
    # least `X.Y.Z`.
    for end_idx in length(uname_strings[3]):-1:5
        try
            return VersionNumber(uname_strings[3][1:end_idx])
        catch e
            if isa(e, InterruptException)
                rethrow(e)
            end
        end
    end

    # We could never parse anything good out of it. :(
    if verbose
        @warn("Unablet to parse a VersionNumber out of uname output", uname_strings)
    end
    return nothing
end

"""
    getuid()

Wrapper around libc's `getuid()` function
"""
getuid() = ccall(:getuid, Cint, ())

_sudo_cmd = nothing
function sudo_cmd()
    global _sudo_cmd

    # Use cached value if we've already run this
    if _sudo_cmd !== nothing
        return _sudo_cmd
    end

    if getuid() == 0
        # If we're already root, don't use any kind of sudo program
        _sudo_cmd = [""]
    elseif success(`sudo -V`)
        # If `sudo` is available, use that
        _sudo_cmd = ["sudo"]
    else
        # Fall back to `su` if all else fails
        _sudo_cmd = ["su", "root", "-c"]
    end
    return _sudo_cmd
end


text:
"""
    plot_Rt_both_SES_groups(model::KenyaCoVSD.CoVAreaModel,forecastdate::Date)

Plot R(t) for both SES groups from a fitted `CoVAreaModel` object.
"""
function plot_Rt_both_SES_groups(model::KenyaCoVSD.CoVAreaModel,forecastdate::Date)
        june1day = (Date(2021,6,1) - Date(2020,2,24)).value

        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])
        qs = quantile(model.MCMC_results.chain)
        ms = mean(model.MCMC_results.chain)
        idx = findfirst(ms[:,1].== :P_eff)
        mean_Peff_perc = round(ms[:,2][idx]*100,digits = 1)
        lower_CI = round(qs[:,2][idx]*100,digits = 1)
        higher_CI = round(qs[:,end][idx]*100,digits = 1)

        Rt_fits = get_Rt(model,forecastdate)
        plt = plot(Rt_fits.Rt₁_fit.pred,
                lab = "Lower SES group, $(mean_Peff_perc)% of county ($(lower_CI)% - $(higher_CI)%)",
                ribbon = (Rt_fits.Rt₁_fit.lb,Rt_fits.Rt₁_fit.ub),
                xticks = (xticktimes[1:2:(end-2)],xticklabs[1:2:(end-2)]),
                title = "Effective reproductive number by SES group: $(model.areaname)",
                size = (700,500),dpi = 250,
                xlims = (-5,june1day),
                ylabel = "R(t)",
                legendfont = 10,titlefont = 15,xtickfontsize=10,ytickfontsize=13,guidefont = 18,
                left_margin = 10mm,right_margin = 7.5mm)

        plot!(plt,Rt_fits.Rt₂_fit.pred,
                ribbon = (Rt_fits.Rt₂_fit.lb,Rt_fits.Rt₂_fit.ub),
                lab = "Higher SES group")
        return plt
end

function plot_ct(ct_fitted)
        june1day = (Date(2021,6,1) - Date(2020,2,20)).value
        ct_for_plot = copy(ct_fitted)
        if june1day > length(ct_fitted)
                ct_for_plot = vcat(ct_fitted,fill(ct_fitted[end], june1day - length(ct_fitted)))
        end

        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])

        ct_plt = plot(ct_for_plot,
                color = :black,lw = 3,
                lab = "",
                xticks = (xticktimes,xticklabs),
                legend = :topright,
                title = "Nairobi one-group fitted contact rates",
                xlims = (-5,june1day),
                ylabel = "Contact rate (rel. to pre-pandemic baseline)",
                size = (1100,500),dpi = 250,
                legendfont = 13,titlefont = 20,tickfontsize=10,guidefont = 12,
                left_margin = 10mm,right_margin = 7.5mm)

                return ct_plt
end
function plot_PCR(fit,linelist_data_with_pos_neg,linelist_data)
        june1day = (Date(2021,6,1) - Date(2020,2,24)).value

        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])

        n = size(fit.pred.mean_PCR_forecast,1)
        n_1 = size(linelist_data_with_pos_neg.cases,1) - 14
        county_pcr_forecast = fit.pred.mean_PCR_forecast[:]
        std_county_pcr_forecast =fit.pred.std_PCR_forecast[:]
        county_pcr_forecast_mv_av = weekly_mv_av(county_pcr_forecast)
        kenya_pcr_forecast_mv_av_std = sqrt.(weekly_mv_av(std_county_pcr_forecast.^2))
        uprname = uppercase(fit.name)
        county_pos = sum(linelist_data_with_pos_neg.cases[:,linelist_data_with_pos_neg.areas .== uprname,
                :,1],dims = [2,3])[:]
        county_pos_mv_av = weekly_mv_av(county_pos)

        PCR_plt = scatter(county_pos[1:(end)],
                ms = 4,markerstrokewidth = 0,color = :grey,alpha = 0.5,
                xticks = (xticktimes,xticklabs),
                lab = "Daily cases: Kenyan linelist (used for fitting)",legend = :topleft,
                title = "$(fit.name) PCR test positives",
                xlims = (-5,june1day),
                ylabel = "Daily PCR-confirmed cases",
                size = (1100,500),dpi = 250,
                legendfont = 13,titlefont = 24,tickfontsize=10,guidefont = 18,
                left_margin = 10mm,right_margin = 7.5mm)

        plot!(PCR_plt,(1+3):(length(county_pos_mv_av)+3),county_pos_mv_av,
                color = :black,lw = 3,
                lab = "Daily cases: 7 day mv-av (used in fitting)")

        smooth_cases_lookahead = weekly_mv_av(sum(linelist_data.cases[(n_1-2):end,linelist_data.areas .== uprname],dims = 2)[:])

        plot!(PCR_plt,(length(county_pos_mv_av)+4):(length(county_pos_mv_av)+3+length(smooth_cases_lookahead)),smooth_cases_lookahead,
                color = :black,lw = 3,ls = :dot,
                lab = "Daily cases: 7 day mv-av (not used in fitting)")

        plot!(PCR_plt,4:(n-3),county_pcr_forecast_mv_av,ribbon = min.(3*kenya_pcr_forecast_mv_av_std,county_pcr_forecast_mv_av),
                color = :red, lw = 2,lab = "Model fit and forecast (7 day mv-av)")

        return PCR_plt

end

function plot_deaths(fit,linelist_data_with_pos_neg)
        june1day = (Date(2021,6,1) - Date(2020,2,24)).value

        n = size(fit.pred.mean_PCR_forecast,1)
        n_1 = size(linelist_data_with_pos_neg.cases,1) - 14

        uprname = uppercase(fit.name)
        county_deaths_forecast = fit.pred.mean_deaths[:]
        std_county_deaths_forecast = fit.pred.std_deaths[:]
        county_deaths_forecast_mv_av = weekly_mv_av(county_deaths_forecast)
        county_deaths_forecast_mv_av_std = sqrt.(weekly_mv_av(std_county_deaths_forecast.^2))
        county_deaths = sum(deaths_data.deaths[:,deaths_data.areas .== uprname],dims = 2)[:]
        county_deaths_mv_av = weekly_mv_av(county_deaths)
        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])

        deaths_plt = scatter(county_deaths,
                ms = 4,markerstrokewidth = 0,color = :grey,alpha = 0.5,
                xticks = (xticktimes,xticklabs),
                lab = "Daily deaths: Kenyan linelist",legend = :topleft,
                title = "$(fit.name) PCR-confirmed deaths",
                xlims = (-5,june1day),
                ylabel = "Daily PCR-confirmed deaths",
                size = (1100,500),dpi = 250,
                legendfont = 13,titlefont = 24,xtickfontsize=10,ytickfontsize=13,guidefont = 18,
                left_margin = 10mm,right_margin = 7.5mm)
        plot!(deaths_plt,(1+3):(length(county_deaths_mv_av)+3),county_deaths_mv_av,
                color = :black,lw = 3,
                lab = "Daily deaths: 7 day mv-av")
        plot!(deaths_plt,4:(n-3),county_deaths_forecast_mv_av,ribbon = 3.0.*county_deaths_forecast_mv_av_std,
                color = :green, lw = 5, ls = :dot,lab = "Model fit and forecast (7 day mv-av)")

        plot!(deaths_plt,(1+3):(length(county_deaths_mv_av)+3),cumsum(county_deaths_mv_av),
                xticks = (xticktimes[4:4:end],xticklabs[4:4:end]),
                xlims = (-5,june1day),
                color = :black,lw = 3,
                grid = nothing,
                inset = (1,bbox(0.35, -0.25, 0.25, 0.25, :center)),
                lab="",
                subplot = 2,
                title = "Cumulative confirmed deaths",
                bg_inside = nothing)
        plot!(deaths_plt,4:(n-3),cumsum(county_deaths_forecast_mv_av),
                xlims = (-5,june1day),color = :green, lw = 5, ls = :dot,lab = "",subplot=2)
        return deaths_plt
end

function plot_pop_exposure(fit,serological_data,serology_data,N_kenya)
        june1day = (Date(2021,6,1) - Date(2020,2,24)).value

        n = size(fit.pred.mean_PCR_forecast,1)
        n_1 = size(linelist_data_with_pos_neg.cases,1) - 14

        uprname = uppercase(fit.name)

        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])

        zeropadtomonday = dayofweek(Date(2020,2,20)) - 1
        county_sero_pos_rnd1_2 = vcat(zeros(zeropadtomonday),
                                sum(serological_data.serodata[:,serological_data.areas .== uprname,:,1],dims = [2,3])[:])
        county_sero_neg_rnd1_2 = vcat(zeros(zeropadtomonday),
                                sum(serological_data.serodata[:,serological_data.areas .== uprname,:,2],dims = [2,3])[:])
        county_weekly_sero_pos_rnd1_2 = [sum(grp) for grp in Iterators.partition(county_sero_pos_rnd1_2,7)]
        county_weekly_sero_total_rnd1_2 = [sum(grp) for grp in Iterators.partition(county_sero_pos_rnd1_2.+county_sero_neg_rnd1_2,7)]

        county_sero_pos_rnd3 = vcat(zeros(zeropadtomonday),
                                sum(serology_data.sero[:,serology_data.areas .== uprname,:,1],dims = [2,3])[:])
        county_sero_neg_rnd3 = vcat(zeros(zeropadtomonday),
                                sum(serology_data.sero[:,serology_data.areas .== uprname,:,2],dims = [2,3])[:])
        county_weekly_sero_pos_rnd3 = [sum(grp) for grp in Iterators.partition(county_sero_pos_rnd3,7)]
        county_weekly_sero_total_rnd3 = [sum(grp) for grp in Iterators.partition(county_sero_pos_rnd3.+county_sero_neg_rnd3,7)]

        #Jeffery intervals
        seroidxs = county_weekly_sero_total_rnd3 .> 0
        # rnd_1_2idx =
        uerr = [invlogcdf(Beta(pos + 0.5,county_weekly_sero_total_rnd3[k] - pos + 0.5),log(0.975)) - pos/county_weekly_sero_total_rnd3[k] for (k,pos) in enumerate(county_weekly_sero_pos_rnd3) ]
        lerr = [pos/county_weekly_sero_total_rnd3[k] - invlogcdf(Beta(pos + 0.5,county_weekly_sero_total_rnd3[k] - pos + 0.5),log(0.025)) for (k,pos) in enumerate(county_weekly_sero_pos_rnd3) ]

        xs_mondays = [-3 + (k-1)*7 for k = 1:length(county_weekly_sero_total_rnd3)]
        rnd1_2_idxs = xs_mondays .< 300
        rnd3_idxs = .~rnd1_2_idxs

        county_sero_pos_nw = fit.pred.mean_serocoverted₁ .+ fit.pred.mean_serocoverted₂
        N = sum(N_kenya[:,fit.name])
        county_serology_forecast_nw = county_sero_pos_nw./N
        std_kenya_serology_forecast = sqrt.(((1.0/sum(N_kenya[:,fit.name]))^2 ).*(fit.pred.std_serocoverted₁.^2 .+ fit.pred.std_serocoverted₂.^2))
        county_infections_forecast = cumsum(fit.pred.mean_incidence₁[1:524] .+ fit.pred.mean_incidence₂[1:524])
        std_county_infections_forecast = sqrt.(cumsum((fit.pred.std_incidence₁[1:524] .+ fit.pred.std_incidence₂[1:524]).^2))

        plt_sero = scatter(xs_mondays[seroidxs.*rnd1_2_idxs],county_weekly_sero_pos_rnd3[seroidxs.*rnd1_2_idxs]./county_weekly_sero_total_rnd3[seroidxs.*rnd1_2_idxs],
                lab = "Weekly KNBTS: rounds 1 and 2 (used in fitting)",
                legend = :topleft,
                yerr = (lerr[seroidxs.*rnd1_2_idxs],uerr[seroidxs.*rnd1_2_idxs]),
                xticks = (xticktimes,xticklabs),
                title = "$(fit.name) population exposure",
                size = (1100,500),dpi = 250,
                xlims = xlims = (-5,june1day), ylims = (-0.025,1),
                ylabel = "Proportion of population",
                legendfont = 10,titlefont = 24,xtickfontsize=10,ytickfontsize=13,guidefont = 18,
                left_margin = 10mm,right_margin = 7.5mm)
        scatter!(plt_sero,xs_mondays[seroidxs.*rnd3_idxs],county_weekly_sero_pos_rnd3[seroidxs.*rnd3_idxs]./county_weekly_sero_total_rnd3[seroidxs.*rnd3_idxs],
                yerr = (lerr[seroidxs.*rnd3_idxs],uerr[seroidxs.*rnd3_idxs]),
                lab = "Weekly KNBTS: round 3 (not used in fitting)")
        plot!(plt_sero,county_serology_forecast_nw,lw = 2,color = :green,
                ribbon = min.(3*std_kenya_serology_forecast,county_serology_forecast_nw),
                lab = "Model fit: seroposivity")

        plot!(plt_sero,county_infections_forecast./N,
                ribbon = min.(9*std_county_infections_forecast./N,county_infections_forecast./N),
                lab = "Model fit: Overall population exposure",
                color = :red)

        return plt_sero
end

function plot_group_incidence(fit)
        june1day = (Date(2021,6,1) - Date(2020,2,24)).value

        xticktimes = [((Date(2020,2,1) + Month(k))- Date(2020,2,24)).value for k = 1:18 ]
        xticklabs = [monthname(k)[1:3]*"/20" for k = 3:12]
        xticklabs = vcat(xticklabs,[monthname(k)[1:3]*"/21" for k = 1:8])
        county_group1_incidence = fit.pred.mean_incidence₁
        std_county_group1_incidence = fit.pred.std_incidence₁
        county_group2_incidence = fit.pred.mean_incidence₂
        std_county_group2_incidence = fit.pred.std_incidence₂

        deleteat!(county_group1_incidence,county_group1_incidence.==0)
        deleteat!(std_county_group1_incidence,std_county_group1_incidence.==0)
        deleteat!(county_group2_incidence,county_group2_incidence.==0)
        deleteat!(std_county_group2_incidence,std_county_group2_incidence.==0)

        if length(std_county_group2_incidence) < length(county_group2_incidence)
                m = length(county_group2_incidence) - length(std_county_group2_incidence)
                std_county_group2_incidence = [std_county_group2_incidence;fill(0.0,m)]
        end

         if length(std_county_group2_incidence) > length(county_group2_incidence)
                m = length(county_group2_incidence) - length(std_county_group2_incidence)
                county_group2_incidence = [county_group2_incidence;fill(0.0,abs(m))]
        end

        plt_inc = plot(county_group1_incidence./1e5,
                lab = "Daily incidence: Lower SES",
                legend = :topleft,
                ribbon = min.(3*std_county_group1_incidence./1e5,county_group1_incidence./1e5),
                xticks = (xticktimes,xticklabs),
                title = "$(fit.name) transmission rates by SES group",
                size = (1100,500),dpi = 250,
                xlims = (-5,june1day),
                ylabel = "Daily infections (100,000s)",
                legendfont = 13,titlefont = 24,xtickfontsize=10,ytickfontsize=13,guidefont = 18,
                left_margin = 10mm,right_margin = 7.5mm)

        plot!(plt_inc,county_group2_incidence./1e5,
                lab = "Daily incidence: Higher SES",
                legend = :topleft,
                ribbon = min.(3*std_county_group2_incidence./1e5,county_group2_incidence./1e5) )
        return plt_inc
end


text:
module GaussianMixtureRegressionDrivers

using AutomotiveDrivingModels
using Distributions
using Compat
using PyCall

import AutomotiveDrivingModels: get_name, action_context, reset_hidden_state!, observe!

const SKLEARN_MIXTURE = PyCall.PyNULL()
function __init__()
    copy!(SKLEARN_MIXTURE, pyimport("sklearn.mixture"))
end

export
    GaussianMixtureRegressionDriver,
    GMR,

    n_learned_components,
    nsuffstats,

    GMRTrainParams,
    calc_bic_score

include(Pkg.dir("AutoDrivers", "src", "GaussianMixtureRegressionDrivers", "gaussian_mixture_regression.jl"))
include(Pkg.dir("AutoDrivers", "src", "GaussianMixtureRegressionDrivers", "gaussian_mixture_regression_drivers.jl"))
include(Pkg.dir("AutoDrivers", "src", "GaussianMixtureRegressionDrivers", "learning.jl"))
include(Pkg.dir("AutoDrivers", "src", "GaussianMixtureRegressionDrivers", "io.jl"))

end


text:
@testset "unimodal" begin
    #= Tests a histogram in which the bin directly following the maximum has an
    intensity of 1,therefore making the distance from the second bin to
    the line always the longest
    =#
    edges = 0:255
    samples = collect(255:-1:0)
    samples[2] = 1
    t = find_threshold(UnimodalRosin(),samples,edges)
    @test t == 1

    # Tests using a unimodal test image from the TestImages package
    img = testimage("moonsurface")
    edges, samples = build_histogram(img, 256)
    t = find_threshold(UnimodalRosin(),samples[1:end], edges)
    @test t == 0.5530790786724538

    #= Tests a histogram that resembles a unimodal image, using a
    surge function. Also tests the scenario in which there are no bins
    with an intensity of 0. In this case, the algorithm should use the
    last bin as the minimum.
    =#
    edges = 0:255
    samples = map(x->surge(x, A = 20 ,k=0.04),edges)
    t = find_threshold(UnimodalRosin(),samples,edges)
    @test t == 112

    edges = 0:200
    samples = collect(200:-1:0)
    samples[2] = 1
    t = find_threshold(UnimodalRosin(),samples,edges)
    @test t == 1
end


text:
# runner.jl
#
# Launch a series of runs to generate some performance logs

using DataFrames, Dates, OrderedCollections

    #CM_dir = joinpath(dirname(@__DIR__), "..", "ClimateMachine.jl")

function main(args::Array{String})
    df1 = parse_sbatch_log(args[1])
    sort!(df1, order(:name))
    df2 = parse_sbatch_log(args[2])
    sort!(df2, order(:name))
    @assert all(df1.name .== df2.name)
    dfd = df1.tottime ./ df2.tottime
    dfn = df1.name
    println(DataFrame(name = dfn, df1 = df1.tottime, df2 = df2.tottime, diff = dfd))
    return (dfd, dfn)
end