text:
using StringDistances, Unicode, Test

@testset "Modifiers" begin
	# Partial
	@test Partial(QGram(2))("martha", "marhta") == 6
	@test Partial(QGram(2))("martha", missing) === missing
	@test Partial(Levenshtein())("martha", "marhta") == 2
	@test Partial(RatcliffObershelp())("martha", "marhta") ≈ 0.16666666 atol = 1e-5
	@test Partial(RatcliffObershelp())("martha", "marhtaXXX") ≈ 0.16666666 atol = 1e-5
	@test Partial(RatcliffObershelp())("martha", missing) === missing

	# TokenSort
	@test TokenSort(QGram(2))("martha", "marhta") == 6
	@test TokenSort(QGram(2))("martha", missing) === missing
	@test TokenSort(Levenshtein())("martha", "marhta") == 2
	@test TokenSort(RatcliffObershelp())("martha", "marhta") ≈ 0.16666666 atol = 1e-5

	# TokenSet
	@test TokenSet(QGram(2))("martha", "marhta") == 6
	@test TokenSet(QGram(2))("martha", missing) === missing
	@test TokenSet(Levenshtein())("martha", "marhta") == 2
	@test TokenSet(RatcliffObershelp())("martha", "marhta") ≈ 0.16666666 atol = 1e-5

	# TokenMax
	@test TokenMax(QGram(2))("martha", "marhta") ≈ 0.6
	@test TokenMax(QGram(2))("martha", missing) === missing
	@test TokenMax(Levenshtein())("martha", "marhta") ≈ 1/3
	@test TokenMax(RatcliffObershelp())("martha", "marhta") ≈ 0.16666666 atol = 1e-5
end

@testset "Compare" begin
	# Qgram
	@test compare("", "abc", QGram(1)) ≈ 0.0 atol = 1e-4
	@test compare("abc", "cba", QGram(1)) ≈ 1.0 atol = 1e-4
	@test compare("abc", "ccc", QGram(1)) ≈ 1/3 atol = 1e-4
	compare("aüa", "aua", TokenMax(QGram(2)))
	@test compare("", "abc", Jaccard(2)) ≈ 0.0 atol = 1e-4
	@test compare("martha", "martha", Jaccard(2)) ≈ 1.0 atol = 1e-4
	@test compare("martha", "martha", Jaccard(2)) ≈ 1.0 atol = 1e-4
	@test compare("aa", "aa ", Partial(Jaccard(2))) ≈ 1.0
	@test compare("martha", "martha", Cosine(2)) ≈ 1.0 atol = 1e-4
	@test compare("martha", "martha", Overlap(2)) ≈ 1.0 atol = 1e-4
	@test compare("martha", "martha", SorensenDice(2)) ≈ 1.0 atol = 1e-4

	# Jaro
	@test compare("aüa", "aua", Hamming()) ≈ 2/3
	@test compare("aüa", "aua", Jaro()) ≈ 0.77777777 atol = 1e-4
	@test compare("New York Yankees",  "", Partial(Jaro())) ≈ 0.0

	# JaroWinkler
	@test compare("martha", "marhta", JaroWinkler()) ≈ 0.9611 atol = 1e-4
	@test compare("dwayne", "duane", JaroWinkler()) ≈ 0.84 atol = 1e-4
	@test compare("dixon", "dicksonx", JaroWinkler()) ≈ 0.81333 atol = 1e-4
	@test compare("william", "williams", JaroWinkler()) ≈ 0.975 atol = 1e-4
	@test compare("", "foo", JaroWinkler()) ≈ 0.0 atol = 1e-4
	@test compare("a", "a", JaroWinkler()) ≈ 1.0 atol = 1e-4
	@test compare("abc", "xyz", JaroWinkler()) ≈ 0.0 atol = 1e-4

	#Levenshtein
	compare("aüa", "aua", Levenshtein())
	@test compare("ok", missing, Levenshtein()) === missing
	compare("aüa", "aua", OptimalStringAlignement())
	@test StringDistances.Normalized(Partial(OptimalStringAlignement()))("ab", "cde") == 1.0
	@test compare("ab", "de", Partial(OptimalStringAlignement())) == 0

	# RatcliffObershelp
	@test compare("New York Mets vs Atlanta Braves", "", RatcliffObershelp())  ≈ 0.0
	@test round(Int, 100 * compare("为人子女者要堂堂正正做人，千万不可作奸犯科，致使父母蒙羞", "此前稍早些时候中国商务部发布消息称，中美经贸高级别磋商双方牵头人通话，中方就美拟9月1日加征关税进行了严正交涉。", RatcliffObershelp())) == 5
	compare("aüa", "aua", TokenMax(RatcliffObershelp()))
	@test compare("New York Yankees",  "Yankees", Partial(RatcliffObershelp())) ≈ 1.0
	@test compare("New York Yankees",  "", Partial(RatcliffObershelp())) ≈ 0.0
	#@test compare("mariners vs angels", "los angeles angels at seattle mariners", Partial(RatcliffObershelp())) ≈ 0.444444444444
	@test compare("HSINCHUANG", "SINJHUAN", Partial(RatcliffObershelp())) ≈ 0.875
	@test compare("HSINCHUANG", "LSINJHUANG DISTRIC", Partial(RatcliffObershelp())) ≈ 0.8
	@test compare("HSINCHUANG", "SINJHUANG DISTRICT", Partial(RatcliffObershelp())) ≈ 0.8
	@test compare("HSINCHUANG",  "SINJHUANG", Partial(RatcliffObershelp())) ≈ 0.8888888888888
	@test compare("New York Mets vs Atlanta Braves", "Atlanta Braves vs New York Mets", TokenSort(RatcliffObershelp())) ≈ 1.0
	@test compare(graphemes("New York Mets vs Atlanta Braves"), graphemes("Atlanta Braves vs New York Mets"), Partial(RatcliffObershelp()))  ≈ compare("New York Mets vs Atlanta Braves", "Atlanta Braves vs New York Mets", Partial(RatcliffObershelp()))
	@test compare("mariners vs angels", "los angeles angels of anaheim at seattle mariners", TokenSet(RatcliffObershelp())) ≈ 1.0 - 0.09090909090909094
	@test compare("New York Mets vs Atlanta Braves", "", TokenSort(RatcliffObershelp()))  ≈ 0.0
	@test compare("mariners vs angels", "", TokenSet(RatcliffObershelp())) ≈ 0.0
	@test compare("mariners vs angels", "los angeles angels at seattle mariners", TokenSet(Partial(RatcliffObershelp()))) ≈ 1.0
	@test compare("mariners", "mariner", TokenMax(RatcliffObershelp())) ≈ 0.933333333333333
	@test compare("为人子女者要堂堂正正做人，千万不可作奸犯科，致使父母蒙羞", "此前稍早些时候中国商务部发布消息称，中美经贸高级别磋商双方牵头人通话，中方就美拟9月1日加征关税进行了严正交涉。", RatcliffObershelp()) ≈ 5 / 100 atol = 1e-2
	@test compare("为人子女者要堂堂正正做人，千万不可作奸犯科，致使父母蒙羞", "此前稍早些时候中国商务部发布消息称，中美经贸高级别磋商双方牵头人通话，中方就美拟9月1日加征关税进行了严正交涉。", Partial(RatcliffObershelp())) ≈ 7 / 100 atol = 1e-2
	@test compare("mariners", "mariner are playing tomorrow", TokenMax(RatcliffObershelp())) ≈ 79 / 100 atol = 1e-2
	@test compare("mariners", "mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow", Partial(RatcliffObershelp())) ≈ 88 / 100 atol = 1e-2
	@test compare("mariners", "mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow", TokenSort(RatcliffObershelp())) ≈ 11 / 100 atol = 1e-2
	@test compare("mariners", "are mariner playing tomorrow", RatcliffObershelp()) ≈ 39 / 100 atol = 1e-2
	@test compare("mariners", "are mariner playing tomorrow", Partial(RatcliffObershelp())) ≈ 88 / 100 atol = 1e-2
	@test compare("mariners", "mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow", TokenSet(RatcliffObershelp())) ≈ 39 / 100 atol = 1e-2
	@test compare("mariners", "mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow", TokenSet(Partial(RatcliffObershelp()))) ≈ 88 / 100 atol = 1e-2
	# not exactly the same because tokenmax has uses the max of rounded tokenset etc
	@test compare("mariners", "mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow mariner are playing tomorrow", TokenMax(RatcliffObershelp())) ≈ 78.75 / 100 atol = 1e-2

	# check min
	strings = [
	("martha", "marhta"),
	("dwayne", "duane") ,
	("dixon", "dicksonx"),
	("william", "williams"),
	("", "foo"),
	("a", "a"),
	("abc", "xyz"),
	("abc", "ccc"),
	("kitten", "sitting"),
	("saturday", "sunday"),
	("hi, my name is", "my name is"),
	("alborgów", "amoniak"),
	("cape sand recycling ", "edith ann graham"),
	( "jellyifhs", "jellyfish"),
	("ifhs", "fish"),
	("leia", "leela"),
	]
	for dist in (Levenshtein, OptimalStringAlignement)
		for i in eachindex(strings)
			if compare(strings[i]..., dist()) <  1 / 3
				@test compare(strings[i]..., dist() ; min_score = 1/ 3) ≈ 0.0
			else
				@test compare(strings[i]..., dist() ; min_score = 1/ 3) ≈ compare(strings[i]..., dist())
			end
		end
	end
end

@testset "Find*" begin
	# findnearest
	@test findnearest("New York", ["NewYork", "Newark", "San Francisco"], Levenshtein()) == ("NewYork", 1)
	@test findnearest("New York", ["San Francisco", "NewYork", "Newark"], Levenshtein()) == ("NewYork", 2)
	@test findnearest("New York", ["Newark", "San Francisco", "NewYork"], Levenshtein()) == ("NewYork", 3)
	@test findnearest("New York", ["NewYork", "Newark", "San Francisco"], Jaro()) == ("NewYork", 1)
	@test findnearest("New York", ["NewYork", "Newark", "San Francisco"], QGram(2)) == ("NewYork", 1)
	@test findnearest("New York", ["Newark", "San Francisco", "NewYork"], QGram(2)) == ("NewYork", 3)

	# findall
	@test findall("New York", ["NewYork", "Newark", "San Francisco"], Levenshtein()) == [1]
	@test findall("New York", ["NewYork", "Newark", "San Francisco"], Jaro()) == [1, 2]
	@test findall("New York", ["NewYork", "Newark", "San Francisco"], Jaro(); min_score = 0.99) == Int[]
	@test findall("New York", ["NewYork", "Newark", "San Francisco"], QGram(2); min_score = 0.99) == Int[]

	if VERSION >= v"1.2.0"
		@test findnearest("New York", skipmissing(["NewYork", "Newark", missing]), Levenshtein()) == ("NewYork", 1)
		@test findnearest("New York", skipmissing(Union{AbstractString, Missing}[missing, missing]), Levenshtein()) == (nothing, nothing)
		@test findall("New York", skipmissing(["NewYork", "Newark", missing]), Levenshtein()) == [1]
		@test findall("New York", skipmissing(Union{AbstractString, Missing}[missing, missing]), Levenshtein()) == []
	end
end


text:
const affine_index = [1]  # which affine vector index to use

reset_affine_index() = affine_index[1] = 1

@doc """
An affine quantity for affine arithmetic.
The usual way to create an affine quantity is from an interval `X`, via `FullAffine(X)`.
"""
struct FullAffine{T<:AbstractFloat}
    c::T   # mid-point
    γ::Vector{T}  # error terms
end

function show(io::IO, C::FullAffine)
    print(io, "⟨", C.c, "; ", C.γ, "⟩")
end

==(C::FullAffine, D::FullAffine) = C.c == D.c && C.γ == D.γ

@doc """
    FullAffine(X::Interval)

Construct a new `FullAffine` quantity from an `Interval`.
"""
function FullAffine(X::Interval)
    c = mid(X)
    r = radius(X)

    index = affine_index[1]
    affine_index[1] += 1

    γ = zeros(index)
    γ[end] = r

    return FullAffine(c, γ)
end

# conversion of numerical constant to affine:
FullAffine(c::Real) = FullAffine(c, Float64[])

range(C::FullAffine) = C.c + sum(abs.(C.γ))*(-1..1)
range(X::Interval) = X

# morally:
# +(C::FullAffine, D::Afine) = FullAffine(C.c + D.c, C.γ + D.γ)

for op in (:+, :-)
    @eval function $op(C::FullAffine, D::FullAffine)
        k = length(C.γ)
        l = length(D.γ)

        # account for unequal lengths:
        m = min(k, l)
        common_γ = $op(C.γ[1:m], D.γ[1:m])

        if l >= k
            γ = [ common_γ; D.γ[m+1:l] ]
        else
            γ = [ common_γ; C.γ[m+1:k] ]
        end

        FullAffine($op(C.c, D.c), γ)
    end
end

+(C::FullAffine, α::Real) = FullAffine(C.c + α, C.γ)
+(α::Real, C::FullAffine) = C + α

-(C::FullAffine, α::Real) = FullAffine(C.c - α, C.γ)
-(α::Real, C::FullAffine) = FullAffine(α - C.c, [-x for x in C.γ])

function *(C::FullAffine, D::FullAffine)

    c = C.c
    d = D.c

    k = length(C.γ)
    l = length(D.γ)

    # account for unequal lengths:
    m = min(k, l)
    common_γ = d * C.γ[1:m] + c * D.γ[1:m]

    error_bound = sum(abs.(C.γ)) * sum(abs.(D.γ))

    if l >= k
        γ = [ common_γ; c*D.γ[m+1:l]; error_bound ]
    else
        γ = [ common_γ; d*C.γ[m+1:k]; error_bound ]
    end

    # |γ|₁ = sum(abs.(γ))

    FullAffine(C.c * D.c, γ)
end

*(α::Real, C::FullAffine) = FullAffine(α*C.c, α*C.γ)
*(C::FullAffine, α::Real) = α * C

^(x::FullAffine, n::Integer) = Base.power_by_squaring(x, n)

Base.literal_pow(::typeof(^), x::FullAffine{T}, ::Val{p}) where {T,p} = x^p

eltype(C::FullAffine{T}) where T = T
zero(C::FullAffine{T}) where T = FullAffine(zero(T))
zero(::Type{FullAffine{T}}) where T = FullAffine(zero(T))

one(C::FullAffine{T}) where T = FullAffine(one(T))
one(::Type{FullAffine{T}}) where T = FullAffine(one(T))

#one(C::FullAffine)

ext: jl

max_stars_repo_name: dpsanders/AffineArithmetic

max_issues_repo_name: dpsanders/AffineArithmetic

max_forks_repo_name: dpsanders/AffineArithmetic


text:
using Dates
include("filtering.jl");
##
odc = ODC(pokes)
odc[!,:ClearDuration] = [r.Reward ? r.PokeDuration > 0.7 : r.PokeDuration < 0.4 for r in eachrow(odc)]
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]

filter!(r -> r.Treatment in list &&
    r.ClearDuration &&
    r.Trial < 61 &&
    r.MouseID != "pc7",# ["pc1","pc2","pc3","pc4","pc5","pc6","pc8","pc9","pc10"],
    odc)

@df odc density(:PokeDuration,group = :Reward,
    linecolor = :auto,  color_palette = [:black,:blue],
    xticks = 0:0.25:4)

odc[odc.Treatment .== "PreVehicle",:Treatment] .= "Control"
odc[odc.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in odc[odc.Treatment .== "Saline",:Stim]]
odc[odc.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
odc[odc.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in odc[odc.Treatment .== "SB242084_opt",:Stim]]
## ODC from last reward
df1 =  filter(r -> !ismissing(r.ODC) &&
    !r.Reward &&
    !ismissing(r.PokeFromLastRew) &&
    -5 < r.PokeFromLastRew < 5 &&
    r.Treatment in ["Control"],
    odc)
# open_html_table(df1[1:200,[:Trial,:PokeFromTrial,:PokeFromLastRew]])
df2 = combine(groupby(df1,:Treatment)) do dd
    summarize(dd,:Poke_Hierarchy,:ODC)
end
Drug_colors!(df2)
@df df2 scatter(:Xaxis, :Mean,
    yerror = :SEM,
    group = :Treatment,
    color = :color)
## ODC Poke From Leaving
phases = ["Citalopram","Optogenetic","SB242084","Altanserin"]
allignment = :PokeFromLeaving
ps1 = plot_odc(odc, phases,allignment; limit = 5)
filename1 = string(allignment) * ".pdf"
savefig(joinpath(figs_loc,"ODC",filename1))
ps2 = plot_wilcoxon_odc(odc, phases,allignment; limit = 5)
filename2 = string(allignment) * "Test.pdf"
savefig(joinpath(figs_loc,"ODC",filename2))
##  ODC Poke Hierarchy
allignment = :Poke_Hierarchy
ps1 = plot_odc(odc, phases,allignment; limit = 5)
filename1 = string(allignment) * ".pdf"
savefig(joinpath(figs_loc,"ODC",filename1))
ps2 = plot_wilcoxon_odc(odc, phases,allignment; limit = 5)
filename2 = string(allignment) * "Test.pdf"
savefig(joinpath(figs_loc,"ODC",filename2))
##  ODC Poke From Last Reward
allignment = :PokeFromLastRew
ps1 = plot_odc(odc, phases,allignment; limit = 5)
filename1 = string(allignment) * ".pdf"
savefig(joinpath(figs_loc,"ODC",filename1))
ps2 = plot_wilcoxon_odc(odc, phases,allignment; limit = 5)
filename2 = string(allignment) * "Test.pdf"
savefig(joinpath(figs_loc,"ODC",filename2))
## Delta before-after last reward ODC in control data
df1 = filter(r -> r.Treatment == "Control" &&
    !r.Reward &&
    !ismissing(r.PokeFromLastRew) &&
    r.PokeFromLastRew in [-1,1] &&
    !ismissing(r.ODC), odc)
df2 = combine(:ODC => mean => :ODC, groupby(df1,[:Phase,:MouseID, :PokeFromLastRew]))
df3 = unstack(df2,:PokeFromLastRew,:ODC)
df3[!,:Delta_odc] = df3[:,3] .- df3[:,4]
df4 = dropmissing(df3)
df5 = wilcoxon(df4, Symbol(-1), Symbol(1))
plot_wilcoxon(df5)
filename1 = "DeltaODC_allcontrol.pdf"
savefig(joinpath(figs_loc,"ODC",filename1))
gd = groupby(df4,:Phase)
df6 = combine(gd) do dd
    wilcoxon(dd, Symbol(-1), Symbol(1))
end
plot_wilcoxon(df6)
filename2 = "DeltaODC_separatecontrol.pdf"
savefig(joinpath(figs_loc,"ODC",filename2))
##
odc = ODC(pokes)
odc[!,:ClearDuration] = [r.Reward ? r.PokeDuration > 0.7 : r.PokeDuration < 0.4 for r in eachrow(odc)]
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]

filter!(r -> r.Treatment in list &&
    r.ClearDuration &&
    r.Trial < 61 &&
    r.MouseID != "pc7",# ["pc1","pc2","pc3","pc4","pc5","pc6","pc8","pc9","pc10"],
    odc)

odc[odc.Treatment .== "PreVehicle",:Treatment] .= "Control"
odc[odc.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in odc[odc.Treatment .== "Saline",:Stim]]
odc[odc.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
odc[odc.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in odc[odc.Treatment .== "SB242084_opt",:Stim]]

union(odc.Treatment)

df1 = filter(r -> !r.Reward &&
    r.ExpDay > 3 &&
    !ismissing(r.PokeFromLastRew) &&
    0 < r.PokeFromLastRew &&
    !ismissing(r.ODC) &&
    r.Treatment in ["Control","Citalopram", "Optogenetic", "Altanserin","SB242084"] &&
    !in(r.Phase,["Methysergide", "Way_100135","SB242084_opt"]) ,odc)

categorical!(df1,[:MouseID,:Treatment,:Phase])

levels!(df1.Treatment, ["Control",
    "Citalopram",
    "Optogenetic",
    "Altanserin",
    "SB242084"])

mdl1 = @formula(ODC ~ 1 + PokeFromLastRew + (1|MouseID))
mdl2 = @formula(ODC ~ 1 + PokeFromLastRew * Treatment + (1|MouseID))
fm1 = fit!(LinearMixedModel(mdl1,df1))
Manipulation_effect = fit!(LinearMixedModel(mdl2,df1))
AICc_test = SwitchingAnalysis.AICcTest(Manipulation_effect,fm1)
##
df2 = combine(groupby(df1, :Treatment)) do dd
    summarize(dd,:PokeFromLastRew,:ODC)
end
Drug_colors!(df2)
##
TRM = ["Control", "Citalopram"]
@df filter(r -> r.Treatment  in TRM,df2) plot(:Xaxis, :Mean, ribbon = :SEM,
    fillalpha = 0.3, linewidth = 3, xlims = (1,12), ylims = (0.5,1),
    group = :Treatment, linecolor = :color, color = :color, xlabel = "Poke from last reward", ylabel = "ODC")
filename = TRM[2] * "ODCRegression.pdf"
savefig(joinpath(figs_loc,"ODC","Regressions",filename))
##
TRM = ["Control", "Optogenetic"]
@df filter(r -> r.Treatment  in TRM,df2) plot(:Xaxis, :Mean, ribbon = :SEM,
    fillalpha = 0.3, linewidth = 3, xlims = (1,12), ylims = (0.5,1),
    group = :Treatment, linecolor = :color, color = :color, xlabel = "Poke from last reward", ylabel = "ODC")
filename = TRM[2] * "ODCRegression.pdf"
savefig(joinpath(figs_loc,"ODC","Regressions",filename))
##
TRM = ["Control", "Altanserin"]
@df filter(r -> r.Treatment  in TRM,df2) plot(:Xaxis, :Mean, ribbon = :SEM,
    fillalpha = 0.3, linewidth = 3, xlims = (1,12), ylims = (0.5,1),
    group = :Treatment, linecolor = :color, color = :color, xlabel = "Poke from last reward", ylabel = "ODC")
filename = TRM[2] * "ODCRegression.pdf"
savefig(joinpath(figs_loc,"ODC","Regressions",filename))
##
TRM = ["Control", "SB242084"]
@df filter(r -> r.Treatment  in TRM,df2) plot(:Xaxis, :Mean, ribbon = :SEM,
    fillalpha = 0.3, linewidth = 3, xlims = (1,12), ylims = (0.5,1),
    group = :Treatment, linecolor = :color, color = :color, xlabel = "Poke from last reward", ylabel = "ODC")
filename = TRM[2] * "ODCRegression.pdf"
savefig(joinpath(figs_loc,"ODC","Regressions",filename))
## Logistic regression
df1[!,:Leave] = [r == 0 for r in df1.PokeFromLeaving]
# full model
m_full = @formula(Leave ~ ODC * PokeFromLastRew * Treatment  + (1|MouseID));
f_full = GeneralizedLinearMixedModel(m_full, df1, Bernoulli());
# simplified model without significant simple factors
m_a1 = @formula(Leave ~ ODC + PokeFromLastRew + ODC & PokeFromLastRew + ODC & Treatment & PokeFromLastRew  + (1|MouseID));
f_a1 = GeneralizedLinearMixedModel(m_a1, df1, Bernoulli());
AICcTest(f_full, f_a1)


text:
using DrWatson
@quickactivate
using HierarchicalAD
using Flux, CUDA
using BSON, FileIO
using ValueHistories, DataFrames
using Plots
using ArgParse
using IterTools

using HierarchicalAD: construct_conv_classifier, train_conv_classifier!, create_input_data
using HierarchicalAD: auc_val, classifier_score

arg_table = ArgParseSettings(;autofix_names=true)
@add_arg_table arg_table begin
	"--modeldir"
		default = "models/contaminated_training/run2/"
		help = "dir where the models are saved"
	"--model-id"
		default = nothing
		help = "id of the FVLAE model whose architecture is to be used for construction fo the classifier"
	"--model-ind"
		default = nothing
		help = "index of the FVLAE model whose architecture is to be used for construction fo the classifier"
end
args = parse_args(arg_table)
@unpack modeldir, model_id, model_ind = args

# get the data
X, y = HierarchicalAD.load_shapes2D();

# set device id
CUDA.device!(2)

# setup paths
modeldir = datadir(modeldir)
plotpath = datadir("plots/baseline_comparison/conv_classifier")
outpath = datadir("baseline_comparison/conv_classifier")
mkpath(outpath)
mkpath(plotpath)

# what models are we gonna use
inmfs = readdir(modeldir)
if !isnothing(model_ind)
	mfs = [inmfs[Meta.parse(model_ind)]]
elseif !isnothing(model_id)
	mfs = filter(x->occursin(model_id, x), mfs)
else
	mfs = inmfs # iterate over all
end

mfs = ["activation=relu_channels=8-8-16-32-64_data=shapes2D_gamma=50.0_lambda=0.0_latent_dim=4_model=fvlae_model_id=20210817174044364_xdist=bernoulli.bson"]

# some params
nns = floor.(Int, 10.0 .^ range(1, 5, length=10))
#nns = floor.(Int, 10.0 .^ range(1, 2, length=2))
nas = floor.(Int, 10.0 .^ range(1, 4, length=4))
#nas = floor.(Int, 10.0 .^ range(1, 2, length=2))
ntst = 25000
label_vec = [
	y[!,:shape] .== 2, y[!,:shape] .!= 2,
	y[!,:scale] .< 0.75, y[!,:scale] .> 0.75,
	y[!,:posX] .> 0.5, y[!,:posY] .> 0.5,
	y[!,:normalized_orientation] .< 0.1,
	.!((y[!,:shape] .== 1) .| ((y[!,:posX] .> 0.5) .& (y[!,:posY] .> 0.5)))
	]
label_desc = [
	"ovals normal", "ovals anomalous",
	"small normal", "large normal",
	"right normal", "bottom normal",
	"rotated anomalous",
	"squares or right bottom anomalous"
	]
conv_classifier_params = (lr = 0.001f0, λ = 0f0, batchsize = 64, nepochs = 1000, patience = 1, verb = true)

function train_and_evaluate_classifier(mf)
	model_id = split(split(mf, "model_id=")[2], "_")[1]
	autoencoder_data = load(joinpath(modeldir, mf))

	for seed = 1:10
		aucs_per_split = []
		scores_per_split = []
		tst_ys_per_split = []
		for (labels, desc) in zip(label_vec, label_desc)
			aucs_per_samples = []
			scores_per_samples = []
			tst_ys_per_samples = []
			for (nn, na) in product(nns, nas)
				@info "Training classifier using $(desc), seed=$seed."
				@info "# normal training/validation = $nn, # anomalous training/validation = $na, # test = $(ntst*2)"
				try
					global tr_nX, (tr_X, tr_y), (val_X, val_y), (tst_X, tst_y) = create_input_data(X, y, labels, nn, nn, na, na,
					    ntst; seed=seed);
				catch e
					if isa(e, BoundsError)
						@info "Not enough data, skipping training..."
						push!(aucs_per_samples, NaN)
						push!(scores_per_samples, [NaN])
						push!(tst_ys_per_samples, [NaN])
						continue
					else
						rethrow(e)
					end
				end
				classifier = construct_conv_classifier(
					autoencoder_data[:experiment_args].kernelsizes,
					autoencoder_data[:experiment_args].activation,
					autoencoder_data[:experiment_args].stride,
					autoencoder_data[:experiment_args].channels,
					size(X))

				classifier, history, opt = train_conv_classifier!(classifier, tr_X, tr_y, val_X, val_y; conv_classifier_params...);
				scores = cpu(classifier_score(classifier, tst_X, 128))
				auc = auc_val(tst_y, scores)
				@info "Finished training, test AUC = $(auc)"
				push!(aucs_per_samples, auc)
				push!(scores_per_samples, scores)
				push!(tst_ys_per_samples, tst_y)
			end
			push!(aucs_per_split, aucs_per_samples)
			push!(scores_per_split, scores_per_samples)
			push!(tst_ys_per_split, tst_ys_per_samples)
		end
		savef = joinpath(outpath, "classifier_output_model_id=$(model_id)_seed=$(seed).bson")
		outdata = Dict(
			:nns => nns, :nas => nas, :ntst => ntst,
			:model_id => model_id,
			:conv_classifier_params => conv_classifier_params,
			:aucs => aucs_per_split,
			:scores => scores_per_split,
			:tst_ys => tst_ys_per_split,
			:label_desc => label_desc
			)
		save(savef, outdata)
		@info "Saved to $savef"
	end
end

for mf in mfs
	train_and_evaluate_classifier(mf)
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
using PkgBenchmark
include("pprinthelper.jl")

if length(ARGS) == 2
    group_target = PkgBenchmark.readresults(ARGS[1])
    group_baseline = PkgBenchmark.readresults(ARGS[2])
else
    res_dir = joinpath(@__DIR__, "results")
    last_result = maximum(filter(x -> x != ".gitignore", readdir(res_dir)))
    baselines = filter(x -> (x != ".gitignore") & occursin("baseline", x), readdir(res_dir))
    last_baseline = isempty(baselines) ? maximum(filter(x -> (x != ".gitignore") & (x != last_result), readdir(res_dir))) : maximum(baselines)
    group_target = PkgBenchmark.readresults(joinpath(res_dir, last_result))
    group_baseline = PkgBenchmark.readresults(joinpath(res_dir, last_baseline))
end
judgement = judge(group_target, group_baseline)

displayresult(judgement)

printnewsection("Target result")
displayresult(group_target)

printnewsection("Baseline result")
displayresult(group_baseline)


text:
module Epochs

using EarthOrientation: getΔUT1
using LeapSeconds: offset_tai_utc

import Base: -, +, <, ==, isapprox, isless, show
import Dates
import Dates: format, parse

import ..AstroDates:
    Date,
    DateTime,
    Time,
    calendar,
    date,
    day,
    dayofyear,
    fractionofday,
    hour,
    j2000,
    julian,
    julian_twopart,
    millisecond,
    minute,
    month,
    second,
    secondinday,
    time,
    year,
    yearmonthday

export Epoch,
    CCSDS_EPOCH,
    FIFTIES_EPOCH,
    FUTURE_INFINITY,
    GALILEO_EPOCH,
    GPS_EPOCH,
    J2000_EPOCH,
    JULIAN_EPOCH,
    MODIFIED_JULIAN_EPOCH,
    PAST_INFINITY,
    UNIX_EPOCH,
    date,
    day,
    dayofyear,
    fractionofday,
    hour,
    j2000,
    julian,
    julian_twopart,
    millisecond,
    minute,
    modified_julian,
    month,
    now,
    second,
    secondinday,
    time,
    year,
    yearmonthday

using ..TimeScales
using ..AstroDates
using ..Periods

const J2000_TO_JULIAN = 2.451545e6days
const J2000_TO_MJD = 51544.5days

@inline function two_sum(a, b)
    hi = a + b
    a1 = hi - b
    b1 = hi - a1
    lo = (a - a1) + (b - b1)
    hi, lo
end

struct Epoch{S, T} <: Dates.AbstractDateTime
    epoch::Int64
    offset::T
    ts_offset::T
    function Epoch{S}(epoch::Int64, offset::T, ts_offset::T) where {S, T}
        new{S::TimeScale, T}(epoch, offset, ts_offset)
    end
end

function Epoch{S}(epoch::Int64, offset, ts_offset, Δt) where S
    sum, residual = two_sum(offset, Δt)

    if !isfinite(sum)
        offset′ = sum
        epoch′ = sum < 0 ? typemin(Int64) : typemax(Int64)
    else
        dl = floor(Int64, sum)
        offset′ = (sum - dl) + residual
        epoch′ = epoch + dl
    end

    Epoch{S}(epoch′, offset′, ts_offset)
end

"""
    Epoch{S}(jd1::T, jd2::T=zero(T); origin=:j2000) where {S, T<:Period}

Construct an `Epoch` with time scale `S` from a Julian date
(optionally split into `jd1` and `jd2`). `origin` determines the
variant of Julian date that is used. Possible values are:

- `:j2000`: J2000 Julian date, starts at 2000-01-01T12:00
- `:julian`: Julian date, starts at -4712-01-01T12:00
- `:modified_julian`: Modified Julian date, starts at 1858-11-17T00:00

### Examples ###

```jldoctest
julia> Epoch{UTC}(0.0days, 0.5days)
2000-01-02T00:00:00.000 UTC

julia> Epoch{UTC}(2.451545e6days, origin=:julian)
2000-01-01T12:00:00.000 UTC
```
"""
function Epoch{S}(jd1::T, jd2::T=zero(T), args...; origin=:j2000) where {S, T<:Period}
    jd1 = float(value(days(jd1)))
    jd2 = float(value(days(jd2)))
    if jd2 > jd1
        jd1, jd2 = jd2, jd1
    end

    if origin == :j2000
        # pass
    elseif origin == :julian
        jd1 -= value(J2000_TO_JULIAN)
    elseif origin == :modified_julian
        jd1 -= value(J2000_TO_MJD)
    else
        throw(ArgumentError("Unknown Julian epoch: $origin"))
    end

    jd1 *= SECONDS_PER_DAY
    jd2 *= SECONDS_PER_DAY

    sum, residual = two_sum(jd1, jd2)
    epoch = floor(Int64, sum)
    offset = (sum - epoch) + residual

    ftype = float(eltype(T))
    tai = Epoch{TAI}(epoch, ftype(offset), zero(ftype))
    ts_offset = tai_offset(S, tai, args...)
    ep = Epoch{TAI}(tai, -ts_offset)
    Epoch{S}(ep.epoch, ep.offset, ts_offset)
end

Epoch{S}(ep::Epoch{S}, Δt) where {S} = Epoch{S}(ep.epoch, ep.offset, ep.ts_offset, Δt)

function j2000(ep::Epoch, tai_offset)
    (ep.offset + tai_offset + ep.epoch) / SECONDS_PER_DAY * days
end
julian(ep::Epoch, tai_offset) = j2000(ep, tai_offset) + J2000_TO_JULIAN
modified_julian(ep::Epoch, tai_offset) = j2000(ep, tai_offset) + J2000_TO_MJD

function julian_twopart(ep::Epoch, tai_offset)
    jd = value(julian(ep, tai_offset))
    jd1 = trunc(jd)
    jd2 = jd - jd1
    (jd1 * days, jd2 * days)
end

"""
    j2000(ep)

Returns the J2000 Julian date for epoch `ep`.

### Example ###

```jldoctest
julia> j2000(UTCEpoch(2000, 1, 1, 12))
0.0 days
```
"""
j2000(ep::Epoch) = j2000(ep, ep.ts_offset)

"""
    julian(ep)

Returns the Julian Date for epoch `ep`.

### Example ###

```jldoctest
julia> julian(UTCEpoch(2000, 1, 1, 12))
2.451545e6 days
```
"""
julian(ep::Epoch) = julian(ep, ep.ts_offset)

"""
    modified_julian(ep)

Returns the Modified Julian Date for epoch `ep`.

### Example ###

```jldoctest
julia> modified_julian(UTCEpoch(2000, 1, 1, 12))
51544.5 days
```
"""
modified_julian(ep::Epoch) = modified_julian(ep, ep.ts_offset)

"""
    julian_twopart(ep)

Returns the two-part Julian date for epoch `ep`, which is a tuple consisting
of the Julian day number and the fraction of the day.

### Example ###

```jldoctest
julia> julian_twopart(UTCEpoch(2000, 1, 2))
(2.451545e6 days, 0.5 days)
```
"""
julian_twopart(ep::Epoch) = julian_twopart(ep, ep.ts_offset)

"""
    j2000(scale, ep)

Returns the J2000 Julian date for epoch `ep` within a specific time `scale`.

### Example ###

```jldoctest
julia> j2000(TAI, TTEpoch(2000, 1, 1, 12, 0, 32.184))
0.0 days
```
"""
j2000(scale, ep::Epoch) = j2000(ep, tai_offset(scale, ep))

"""
    julian(scale, ep)

Returns the Julian Date for epoch `ep` within a specific time `scale`.

### Example ###

```jldoctest
julia> julian(TAI, TTEpoch(2000, 1, 1, 12, 0, 32.184))
2.451545e6 days
```
"""
julian(scale, ep::Epoch) = julian(ep, tai_offset(scale, ep))

"""
    modified_julian(scale, ep)

Returns the Modified Julian Date for epoch `ep` within a specific time `scale`.

### Example ###

```jldoctest
julia> modified_julian(TAI, TTEpoch(2000, 1, 1, 12, 0, 32.184))
51544.5 days
```
"""
modified_julian(scale, ep::Epoch) = modified_julian(ep, tai_offset(scale, ep))

"""
    julian_twopart(scale, ep)

Returns the two-part Julian date for epoch `ep` within a specific time `scale`,
which is a tuple consisting of the Julian day number and the fraction of the day.

### Example ###

```jldoctest
julia> julian_twopart(TAI, TTEpoch(2000, 1, 1, 12, 0, 32.184))
(2.451545e6 days, 0.0 days)
```
"""
julian_twopart(scale, ep::Epoch) = julian_twopart(ep, tai_offset(scale, ep))

include("offsets.jl")
include("accessors.jl")

show(io::IO, ep::Epoch{S}) where {S} = print(io, DateTime(ep), " ", S)

function Epoch{S}(date::Date, time::Time, args...) where S
    seconds = second(Float64, time)
    ts_offset = tai_offset(S, date, time, args...)

    sum, residual = two_sum(seconds, ts_offset)
    dl = floor(Int64, sum)

    offset = (sum - dl) + residual
    epoch  = Int64(60) * ((j2000(date) * Int64(24) + hour(time)) * Int64(60)
                          + minute(time) - Int64(720)) + dl
    from_tai = tai_offset(S, Epoch{TAI}(epoch, offset, 0.0), args...)
    Epoch{S}(epoch, offset, from_tai)
end

Dates.default_format(::Type{Epoch}) = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss ttt")

"""
    Epoch(str[, format])

Construct an `Epoch` from a string `str`. Optionally a `format` definition can
be passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
object or as a string. In addition to the character codes supported by `DateFormat` the character
code `D` is supported which is parsed as "day of year" (see the example below) and the character
code `t` which is parsed as the time scale.  The default format is `yyyy-mm-ddTHH:MM:SS.sss ttt`.

**Note:** Please be aware that this constructor requires that the time scale is part of `str`, e.g.
`2018-02-06T00:00 UTC`. Otherwise use an explicit constructor, e.g. `Epoch{UTC}`.

### Example ###

```jldoctest
julia> Epoch("2018-02-06T20:45:00.0 UTC")
2018-02-06T20:45:00.000 UTC

julia> Epoch("2018-037T00:00 UTC", "yyyy-DDDTHH:MM ttt")
2018-02-06T00:00:00.000 UTC
```
"""
Epoch(str::AbstractString, format::Dates.DateFormat=Dates.default_format(Epoch)) = parse(Epoch, str, format)

Epoch(str::AbstractString, format::AbstractString) = Epoch(str, Dates.DateFormat(format))

Dates.default_format(::Type{Epoch{S}}) where {S} = Dates.ISODateTimeFormat

"""
    Epoch{S}(str[, format]) where S

Construct an `Epoch` with time scale `S` from a string `str`. Optionally a `format` definition can
be passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
object or as a string. In addition to the character codes supported by `DateFormat` the code `D` can
be used which is parsed as "day of year" (see the example below).  The default format is
`yyyy-mm-ddTHH:MM:SS.sss`.

### Example ###

```jldoctest
julia> Epoch{UTC}("2018-02-06T20:45:00.0")
2018-02-06T20:45:00.000 UTC

julia> Epoch{UTC}("February 6, 2018", "U d, y")
2018-02-06T00:00:00.000 UTC

julia> Epoch{UTC}("2018-037T00:00", "yyyy-DDDTHH:MM")
2018-02-06T00:00:00.000 UTC
```
"""
Epoch{S}(str::AbstractString,
         format::Dates.DateFormat=Dates.default_format(Epoch{S})) where {S} = parse(Epoch{S}, str, format)

Epoch{S}(str::AbstractString, format::AbstractString) where {S} = Epoch{S}(str, Dates.DateFormat(format))

Epoch{S}(d::Date) where {S} = Epoch{S}(d, AstroDates.H00)

Epoch{S}(dt::DateTime) where {S} = Epoch{S}(date(dt), time(dt))
Epoch{S}(dt::Dates.DateTime) where {S} = Epoch{S}(DateTime(dt))

"""
    now()

Get the current date and time as a `UTCEpoch`.
"""
now() = UTCEpoch(Dates.now())

"""
    Epoch{S}(year, month, day, hour=0, minute=0, second=0.0) where S

Construct an `Epoch` with time scale `S` from date and time components.

### Example ###

```jldoctest
julia> Epoch{UTC}(2018, 2, 6, 20, 45, 0.0)
2018-02-06T20:45:00.000 UTC

julia> Epoch{UTC}(2018, 2, 6)
2018-02-06T00:00:00.000 UTC
```
"""
function Epoch{S}(year::Int, month::Int, day::Int, hour::Int=0,
                  minute::Int=0, second::Float64=0.0, args...) where S
    Epoch{S}(Date(year, month, day), Time(hour, minute, second), args...)
end

function Epoch{S}(year::Int64, month::Int64, day::Int64, dayofyear::Int64,
                  hour::Int64, minute::Int64, second::Int64, milliseconds::Int64) where S
    if dayofyear != 0
        date = Date(year, dayofyear)
    else
        date = Date(year, month, day)
    end
    Epoch{S}(date, Time(hour, minute, second + 1e-3milliseconds))
end

function Epoch(year::Int64, month::Int64, day::Int64, dayofyear::Int64,
               hour::Int64, minute::Int64, second::Int64, milliseconds::Int64,
               scale::S) where S<:TimeScale
    if scale === TimeScales.NotATimeScale()
        throw(ArgumentError("Could not parse the provided string as an `Epoch`." *
                            " No time scale was provided."))
    end

    if dayofyear != 0
        date = Date(year, dayofyear)
    else
        date = Date(year, month, day)
    end
    Epoch{scale}(date, Time(hour, minute, second + 1e-3milliseconds))
end

"""
    Epoch{S}(Δtai, ep::TAIEpoch) where S

Convert `ep`, a `TAIEpoch`, to an `Epoch` with time scale `S` by overriding
the offset between `S2` and `TAI` with `Δtai`.

### Examples ###

```jldoctest
julia> ep = TAIEpoch(2000,1,1)
2000-01-01T00:00:00.000 TAI

julia> TTEpoch(32.184, ep)
2000-01-01T00:00:32.184 TT
```
"""
function Epoch{S}(Δtai, ep::Epoch{TAI}) where S
    Epoch{S}(ep.epoch, ep.offset, Δtai)
end

"""
    Epoch{S2}(ep::Epoch{S1}) where {S1, S2}

Convert `ep`, an `Epoch` with time scale `S1`, to an `Epoch` with time
scale `S2`.

### Examples ###

```jldoctest
julia> ep = TTEpoch(2000,1,1)
2000-01-01T00:00:00.000 TT

julia> TAIEpoch(ep)
1999-12-31T23:59:27.816 TAI
```
"""
function Epoch{S2}(ep::Epoch{S1}, args...) where {S1, S2}
    Epoch{S2}(ep.epoch, ep.offset, tai_offset(S2, ep, args...))
end

Epoch{TAI}(ep::Epoch) = Epoch{TAI}(ep.epoch, ep.offset, 0.0)

Epoch{S, T}(ep::Epoch{S, T}) where {S, T} = ep

function isapprox(a::T, b::T; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps()) where T <: Epoch
    sum_a, residual_a = two_sum(a.ts_offset, a.offset)
    Δep_a = floor(Int64, sum_a)
    epoch_a = a.epoch + Δep_a
    offset_a = (sum_a - Δep_a) + residual_a

    sum_b, residual_b = two_sum(b.ts_offset, b.offset)
    Δep_b = floor(Int64, sum_b)
    epoch_b = b.epoch + Δep_b
    offset_b = (sum_b - Δep_b) + residual_b

    epoch_a == epoch_b && isapprox(offset_a, offset_b; atol=atol, rtol=rtol)
end

function isapprox(a::Epoch, b::Epoch; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps())
    a.epoch == b.epoch && isapprox(a.offset, b.offset; atol=atol, rtol=rtol)
end

function ==(a::Epoch, b::Epoch)
    a.epoch == b.epoch && a.offset == b.offset && a.ts_offset == b.ts_offset
end

<(ep1::Epoch, ep2::Epoch) = value(ep1 - ep2) < 0.0
isless(ep1::Epoch, ep2::Epoch) = isless(value(ep1 - ep2), 0.0)

+(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, value(seconds(p)))
-(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, -value(seconds(p)))

"""
    -(a::Epoch, b::Epoch)

Return the duration between epoch `a` and epoch `b`.

### Examples ###

```jldoctest```
julia> UTCEpoch(2018, 2, 6, 20, 45, 20.0) - UTCEpoch(2018, 2, 6, 20, 45, 0.0)
20.0 seconds
```
"""
function -(a::T, b::T) where T <: Epoch
    sum_a, residual_a = two_sum(a.ts_offset, a.offset)
    Δep_a = floor(Int64, sum_a)
    epoch_a = a.epoch + Δep_a
    offset_a = (sum_a - Δep_a) + residual_a

    sum_b, residual_b = two_sum(b.ts_offset, b.offset)
    Δep_b = floor(Int64, sum_b)
    epoch_b = b.epoch + Δep_b
    offset_b = (sum_b - Δep_b) + residual_b

    ((epoch_a - epoch_b) + (offset_a - offset_b)) * seconds
end

-(a::Epoch, b::Epoch) = ((a.epoch - b.epoch) + (a.offset - b.offset)) * seconds

# Generate aliases for all defined time scales so we can use
# e.g. `TTEpoch` instead of `Epoch{TT}`
for scale in TimeScales.ACRONYMS
    epoch = Symbol(scale, "Epoch")
    name = string(epoch)
    @eval begin
        const $epoch = Epoch{$scale}
        export $epoch

        """
            $($name)(str[, format])

        Construct a $($name) from a string `str`. Optionally a `format` definition can be
        passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
        object or as a string. In addition to the character codes supported by `DateFormat`
        the code `D` is supported which is parsed as "day of year" (see the example below).
        The default format is `yyyy-mm-ddTHH:MM:SS.sss`.

        ### Example ###

        ```jldoctest
        julia> $($name)("2018-02-06T20:45:00.0")
        2018-02-06T20:45:00.000 $($scale)

        julia> $($name)("February 6, 2018", "U d, y")
        2018-02-06T00:00:00.000 $($scale)

        julia> $($name)("2018-37T00:00", "yyyy-DDDTHH:MM")
        2018-02-06T00:00:00.000 $($scale)
        ```
        """
        $epoch(::AbstractString)

        """
            $($name)(jd1::T, jd2::T=zero(T); origin=:j2000) where T<:Period

        Construct a $($name) from a Julian date (optionally split into
        `jd1` and `jd2`). `origin` determines the variant of Julian
        date that is used. Possible values are:

        - `:j2000`: J2000 Julian date, starts at 2000-01-01T12:00
        - `:julian`: Julian date, starts at -4712-01-01T12:00
        - `:modified_julian`: Modified Julian date, starts at 1858-11-17T00:00

        ### Examples ###

        ```jldoctest
        julia> $($name)(0.0days, 0.5days)
        2000-01-02T00:00:00.000 $($scale)

        julia> $($name)(2.451545e6days, origin=:julian)
        2000-01-01T12:00:00.000 $($scale)
        ```
        """
        $epoch(::Number, ::Number)

        """
            $($name)(year, month, day, hour=0, minute=0, second=0.0)

        Construct a $($name) from date and time components.

        ### Example ###

        ```jldoctest
        julia> $($name)(2018, 2, 6, 20, 45, 0.0)
        2018-02-06T20:45:00.000 $($scale)

        julia> $($name)(2018, 2, 6)
        2018-02-06T00:00:00.000 $($scale)
        ```
        """
        $epoch(::Int, ::Int, ::Int)
    end
end

include("leapseconds.jl")
include("range.jl")

const JULIAN_EPOCH = TTEpoch(AstroDates.JULIAN_EPOCH, AstroDates.H12)
const J2000_EPOCH = TTEpoch(AstroDates.J2000_EPOCH, AstroDates.H12)
const MODIFIED_JULIAN_EPOCH = TTEpoch(AstroDates.MODIFIED_JULIAN_EPOCH, AstroDates.H00)
const FIFTIES_EPOCH = TTEpoch(AstroDates.FIFTIES_EPOCH, AstroDates.H00)
const CCSDS_EPOCH = TTEpoch(AstroDates.CCSDS_EPOCH, AstroDates.H00)
const GALILEO_EPOCH = TTEpoch(AstroDates.GALILEO_EPOCH, AstroDates.H00)
const GPS_EPOCH = TTEpoch(AstroDates.GPS_EPOCH, AstroDates.H00)
const UNIX_EPOCH = TAIEpoch(AstroDates.UNIX_EPOCH, Time(0, 0, 10.0))

const PAST_INFINITY = TAIEpoch(UNIX_EPOCH, -Inf)
const FUTURE_INFINITY = TAIEpoch(UNIX_EPOCH, Inf)

const EPOCH_77 = TAIEpoch(1977, 1, 1)

function Dates.validargs(::Type{Epoch}, y::Int64, m::Int64, d::Int64,
                         h::Int64, mi::Int64, s::Int64, ms::Int64, ts::S) where S<:TimeScale
    err = Dates.validargs(Dates.DateTime, y, m, d, h, mi, s, ms)
    err !== nothing || return err
    return Dates.argerror()
end

abstract type DayOfYearToken end

@inline function Dates.tryparsenext(d::Dates.DatePart{'D'}, str, i, len, locale)
    next = Dates.tryparsenext_base10(str, i, len, 1, 3)
    next === nothing && return nothing
    val, i = next
    (val >= 1 && val <= 366) || throw(ArgumentError("Day number must be within 1 and 366."))
    return val, i
end

function Dates.format(io, d::Dates.DatePart{'D'}, ep)
    print(io, dayofyear(ep))
end

function Dates.format(io, d::Dates.DatePart{'t'}, ep)
    print(io, timescale(ep))
end

end


text:
"""
    QC(from, to; geno = 0.1, maf = 0.01, hwe = 0.0001, mind = 0.1)
---
Quality control of a genotype dataset in `plink` format.

# Example
    QC("pathfrom/from", "pathto/to")
This function will do a quality control of genotypes in file `from.bed`, `from.bim`,
`from.fam` in path `pathfrom`.  It will save the results in `pathto`, with stem file name as `to`.  By default, the function will remove:
- loci who has more than 10% genotypes missing
- loci of maf ≤ 0.01
- loci of P-value ≤ 0.0001
- ID who has ≥ 10% genotypes missing
"""
function QC(from, to; geno = 0.1, maf = 0.01, hwe = 0.0001, mind = 0.1)
    title("Quality control of genotype set $from")
    item("Parameter check")
    if !isfile(from * ".bed")
        warning("$from files don't exist")
        return
    end
    pathto, target = splitdir(to)
    if length(pathto) == 0
        warning("Target file is in the current folder")
    else
        isdir(pathto) || mkdir(pathto)
    end

    item("Filter ID and SNP")
    _ = read(`$plink  --sheep
			--bfile $from
			--geno $geno
			--maf $maf
			--hwe $hwe
			--mind $mind
			--make-bed
			--out $to`,
                String)
    done()

    item("Summary")
    obtain_set(file) = begin
        tv = String[]
        for line in eachline(file)
            str = split(line)[2]
            push!(tv, str)
        end
        Set(tv)
    end
    fsnp, tsnp = obtain_set("$from.bim"), obtain_set("$to.bim")
    diff = setdiff(fsnp, tsnp)
    if length(diff) > 0
        println(join(diff, ' '))
        msg = "Above " * string(length(diff)) * " SNP were removed"
        warning(msg)
    else
        message("No SNP were removed")
    end
    fid, tid = obtain_set("$from.fam"), obtain_set("$to.fam")
    diff = setdiff(fid, tid)
    if length(diff) > 0
        println(join(diff, ' '))
        msg = "Above " * string(length(diff)) * " ID were removed"
        warning(msg)
    else
        message("No ID were removed")
    end
    done()
end


text:
using Base.Test
const NPROCS = isempty(ARGS) ? 4 : parse(Int, ARGS[1])

function do_mul()
    println("creating initial data...")
    a1 = reshape([1:10;], 1, 10)
    a2 = reshape([1:6;], 6, 1)
    a3 = a2*a1
    a4 = a3'

    println("making blocks...")
    mb = MatOpBlock(a3, a4, :*, NPROCS)
    blk = Block(mb)
    println("running matrix operation...")
    t1 = @elapsed (result = op(blk))

    println("calculating locally...")
    t2 = @elapsed (tr = a3*a4)

    println("verifying...")
    @test_approx_eq tr result

    println("times: local:$t2 parallel:$t1")
end

function do_mul_large()
    println("creating initial data...")
    a1 = reshape([[1:5000;];[1:5000;]], 1, 10000)
    a2 = reshape([[1:3000;];[1:3000;]], 6000, 1)
    a3 = a2*a1
    a4 = a3'

    println("making blocks...")
    mb = MatOpBlock(a3, a4, :*, NPROCS)
    blk = Block(mb)
    println("running matrix operation...")
    t1 = @elapsed (result = op(blk))

    println("calculating locally...")
    t2 = @elapsed (tr = a3*a4)

    println("verifying...")
    @test_approx_eq tr result

    println("times: local:$t2 parallel:$t1")
end

function do_mul_test()
    a1 = rand(10000, 2000)
    a2 = rand(2000, 1)

    mb = MatOpBlock(a1, a2, :*, NPROCS)
    blk = Block(mb)
    println("running matrix operation...")
    t1 = @elapsed (result = op(blk))
    println("times: parallel:$t1")
end

function do_rand_mul()
    println("creating initial data...")
    a1 = RandomMatrix(Float64, (24000, 24000), 0)
    a2 = RandomMatrix(Float64, (24000, 1), 0)

    println("making blocks...")
    mb = MatOpBlock(a1, a2, :*, NPROCS)
    blk = Block(mb)
    println("running matrix operation...")
    t1 = @elapsed (result = op(blk))
    #t2 = 0
    #rmprocs(workers())

    println("calculating locally...")
    la1 = convert(DenseMatrix, mb.mb1)
    la2 = convert(DenseMatrix, mb.mb2)
    t2 = @elapsed (tr = la1*la2)

    println("verifying...")
    @test_approx_eq tr result

    println("times: local:$t2 parallel:$t1")
end

println("adding $NPROCS more processors...")
addprocs(NPROCS)
println("\tnworkers: $(nworkers())")
println("loading packages...")
using Blocks
using Blocks.MatOp

do_mul()
#do_mul_large()
#do_mul_test()
#do_rand_mul()


text:
"""
Parameters for NFW Profile, includes useful quantities defined by parameters
"""
type NFWParameters{G <: AbstractFloat}
    a::G
    d::G
    e::G
    Ec::G
    rlim::G
    b::G
    q::G
    Jb::G
    vmax::G
    rmax::G

    #secondary parameters
    rs::G
    Φs::G
    xlim::G
    Φlim::G

    #adjusted parameter values
    adjEc::G
    adjJb::G
end

function NFWParameters(a::G, d::G, e::G, Ec::G, rlim::G, b::G, q::G, Jb::G, vmax::G, rmax::G) where G <: AbstractFloat

    rs = rmax / 2.16
    Φs = (vmax / 0.465)^2
    xlim = rlim / rs
    adjEc = Φs * Ec
    adjJb = Jb * rs * vmax / 0.465

    p = NFWParameters(a, d, e, Ec, rlim, b, q, Jb, vmax, rmax, rs, Φs, xlim, 0.0, adjEc, adjJb)
    p.Φlim = gravitational_potential(xlim, p)

    return p
end


text:
using Test, JSON

@testset "./src/LocalPolytope/adjacency_decomposition.jl" begin

using BellScenario

test_dir = "./test/integration/files/"

@testset "adjacent_facets" begin
    @testset "41-2-14 polytope" begin
        scenario = LocalSignaling(4,4,2)
        vertices = LocalPolytope.vertices(scenario)

        F = [1,0,0,0,1,0,0,0,1,-1,-1,-1,1]
        adj_facets = LocalPolytope.adjacent_facets(vertices, F, dir=test_dir)

        @test length(adj_facets) == 48

        adj_games = map(f -> convert(BellGame, f, scenario, rep="normalized"), adj_facets)
        canonical_games = unique(map(g -> LocalPolytope.generator_facet(g, scenario), adj_games))

        @test canonical_games == [
            [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
            [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
            [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0]
        ]
    end
end

@testset "adjacency_decomposition" begin
    @testset  "41-2-14 polytope skip positivity" begin
        scenario = LocalSignaling(4,4,2)
        vertices = LocalPolytope.vertices(scenario)

        BG = BellGame([1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1], 2)

        skip = [BellGame([1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],1)]
        dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, skip_games=skip, dir=test_dir)

        # positvity included, but is skipped in computation
        @test collect(keys(dict)) == [
            [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
            [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1],
            [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0],
            [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1],
        ]

        @test dict[[1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0]]["skipped"]

        @test all(d -> haskey(d[2], "norm_facet"), dict)
        @test all(d -> haskey(d[2], "generator_facet"), dict)
    end

    @testset  "41-2-14 polytope no skips" begin
        scenario = LocalSignaling(4,4,2)
        vertices = LocalPolytope.vertices(scenario)

        BG = BellGame([1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1], 2)

        dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, dir=test_dir)

        @test collect(keys(dict)) == [
            [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
            [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1],
            [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0],
            [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1]
        ]

        @test dict[[1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0]]["skipped"] == false
    end

    @testset "41-2-14 polytope filter out positivity vertices" begin
        scenario = LocalSignaling(4,4,2)
        vertices = LocalPolytope.vertices(scenario, rank_d_only = true)

        BG = BellGame([1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1], 2)

        dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, dir=test_dir)

        @test issetequal( collect(keys(dict)), [
            [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
            [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1],
            [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
            [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0],
            [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1],
            [1 1 1 1;0 0 0 0;0 0 0 0;0 0 0 0],
        ])
    end

    @testset "31-2-13 polytope skip by  max_vertices" begin
        scenario = LocalSignaling(3,3,2)
        vertices = LocalPolytope.vertices(scenario)

        BG  = BellGame([1 0 0;1 0 0;0 0 0], 1)

        dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, max_vertices=5, dir=test_dir)

        @test dict[[1 0 0;0 1 0;0 0 1]]["skipped"]

        @test collect(keys(dict)) == [[1 0 0;1 0 0;0 0 0],[1 0 0;0 1 0;0 0 1]]
    end

    @testset "json logging" begin
        filename = "JSON_logging_test.json"
        if isfile(test_dir*filename)
            rm(test_dir*filename, force=true)
        end

        scenario = LocalSignaling(3,3,2)
        vertices = LocalPolytope.vertices(scenario)

        BG  = BellGame([1 0 0;1 0 0;0 0 0], 1)

        dict = LocalPolytope.adjacency_decomposition(
            vertices, BG, scenario,
            dir=test_dir, log=true, log_filename = filename
        )

        log_json = JSON.parsefile(test_dir*filename)

        @test log_json == Dict{String,Any}(
            "[1 0 0; 0 1 0; 0 0 1]" => Dict{String,Any}(
                "generator_facet" => Any[1, 0, 0, 1, -1, -1, 1],
                "norm_facet" => Any[1, 0, -1, -1, 0, 1, 1],
                "considered" => true,
                "num_vertices" => 6,
                "skipped" => false
            ),
            "[1 0 0; 1 0 0; 0 0 0]" => Dict{String,Any}(
                "generator_facet" => Any[1, 1, 0, 0, 0, 0, 1],
                "norm_facet" => Any[1, 1, 0, 0, 0, 0, 1],
                "considered" => true,
                "num_vertices" => 14,
                "skipped" => false
            )
        )

        rm(test_dir*filename, force=true)
    end
end

end


text:
module BoundaryConditions

export
    BCType, Periodic, Flux, Gradient, Value, NoPenetration, Dirchlet, Neumann,
    BoundaryCondition, bctype, getbc, setbc!,
    CoordinateBoundaryConditions,
    FieldBoundaryConditions, HorizontallyPeriodicBCs, ChannelBCs,
    SolutionBoundaryConditions, HorizontallyPeriodicSolutionBCs, ChannelSolutionBCs,
    TendenciesBoundaryConditions, PressureBoundaryConditions,
    ModelBoundaryConditions, BoundaryFunction,
    apply_z_bcs!, apply_y_bcs!,
    fill_halo_regions!, zero_halo_regions!

using Oceananigans: @hascuda
using Oceananigans.Architectures
using Oceananigans.Fields

@hascuda using CUDAnative

include("boundary_condition_types.jl")
include("boundary_condition.jl")
include("coordinate_boundary_conditions.jl")
include("field_boundary_conditions.jl")
include("solution_boundary_conditions.jl")
include("tracer_boundary_conditions.jl")
include("tendency_boundary_conditions.jl")
include("pressure_boundary_conditions.jl")
include("model_boundary_conditions.jl")
include("boundary_function.jl")
include("show_boundary_conditions.jl")

include("fill_halo_regions.jl")
include("zero_halo_regions.jl")

include("apply_flux_bcs.jl")
include("apply_flux_periodic_no_penetration_bcs.jl")
include("apply_value_gradient_bcs.jl")

end


text:
"""
Provides logging functionality.
"""
module Loggers

using Millboard, Dates, MiniLogging
using Genie

import Base.log
export log

"""
Map Genie log levels to MiniLogging levels
"""
const LOG_LEVEL_MAPPING = Dict(
  :debug    => MiniLogging.DEBUG,
  :info     => MiniLogging.INFO,
  :warn     => MiniLogging.WARN,
  :error    => MiniLogging.ERROR,
  :err      => MiniLogging.ERROR,
  :critical => MiniLogging.CRITICAL
)

"""
    log(message::Union{String,Symbol,Number,Exception}, level::Union{String,Symbol} = :debug; showst = false) :: Nothing

Logs `message` to all registered logger, taking into account the log `level`.
"""
function log(message::Union{String,Symbol,Number,Exception}, level::Union{String,Symbol} = :debug; showst = false) :: Nothing
  message = string(" ", message)
  level = string(level)

  try
    if isempty(get_logger().handlers)
      if Genie.config.log_to_file
        basic_config(LOG_LEVEL_MAPPING[Genie.config.log_level], log_path()) # file logger
        push!(get_logger().handlers, MiniLogging.Handler(stderr, "%Y-%m-%d %H:%M:%S")) # console logger
      else
        basic_config(LOG_LEVEL_MAPPING[Genie.config.log_level], date_format = "%Y-%m-%d %H:%M:%S")
      end
    end
  catch ex
    basic_config(LOG_LEVEL_MAPPING[Genie.config.log_level], date_format = "%Y-%m-%d %H:%M:%S")
  end

  loggo = get_logger()

  if level == "debug"
    MiniLogging.@debug(loggo, message)
  elseif level == "info"
    MiniLogging.@info(loggo, message)
  elseif level == "warn"
    MiniLogging.@warn(loggo, message)
  elseif level == "error" || level == "err"
    MiniLogging.@error(loggo, message)
  elseif level == "critical"
    MiniLogging.@critical(loggo, message)
  else
    MiniLogging.@debug(loggo, message)
  end

  if (level == "error") && showst
    println()
    stacktrace()
  end

  nothing
end

function log_path()
  "$(joinpath(Genie.LOG_PATH, ENV["GENIE_ENV"])).log"
end
function log_path!()
  initlogfile()
  log_path()
end

function initlogfile()
	lp = log_path()
	Base.@info "Logging to file at $(abspath(lp)) \n"
  dirname(lp) |> mkpath
	touch(lp)
end

end


text:
import .Utils.Heap: PriorityQueue, dequeue!

"""
    prim_mst_edges(n, edges, weights, [root = 1])

[Prim's algorithm](https://en.wikipedia.org/wiki/Prim%27s_algorithm)
for minimum spanning tree.
Start from node `root` (will become the root node of the spanning
tree).  Return the parent vector and a Boolean vector indicating for
each edge whether the edge is part of the spanning tree.
"""
function prim_mst_edges(
    n::Integer,
    edges::Vector{E},
    weights::Array{Float64},
    root::Int = 1,
) where {E}
    mem = PrimMstMem(edges, n)
    prim_mst_edges(weights, root, mem)
    # swap root element to front
    mem.selected[1], mem.selected[root] = mem.selected[root], mem.selected[1]
    return mem.parent, view(mem.selected, 2:length(mem.selected))
end

struct PrimMstMem
    dist::Vector{Float64}
    parent::Vector{Int}
    neighbors::IncidenceIndex
    selected::Vector{Int}
    pq::PriorityQueue{Int, Float64}

    PrimMstMem(g::EdgeGraph) =
        PrimMstMem(IncidenceIndex(num_nodes(g), g.edges))

    PrimMstMem(g::Graph) =
        PrimMstMem(EdgeGraph(g))

    PrimMstMem(edges::Vector{E}, n::Integer) where {E} =
        PrimMstMem(IncidenceIndex(n, edges))

    PrimMstMem(head::Vector{I}, tail::Vector{I}, n::Integer) where {I} =
        PrimMstMem(IncidenceIndex(n, head, tail))

    function PrimMstMem(neighbors::IncidenceIndex)
        local n = num_nodes(neighbors)
        dist = Vector{Float64}(undef, n)
        parent = Vector{Int}(undef, n)
        selected = Vector{Int}(undef, n)
        pq = PriorityQueue{Int, Float64}(n)
        return new(dist, parent, neighbors, selected, pq)
    end
end

prim_mst_edges(
    weights::Vector{Float64}, root::Int, edges::Vector{E}, n::Integer,
) where {E} =
    prim_mst_edges(weights, root, PrimMstMem(edges, n))

prim_mst_edges(weights::Vector{Float64}, root::Int, mem::PrimMstMem) =
    prim_mst_edges(
        weights,
        mem.dist,
        mem.parent,
        mem.neighbors,
        mem.selected,
        mem.pq,
        root
    )

function prim_mst_edges(
    edge_weight::Vector{Float64},
    dist::Vector{Float64},
    parent::Vector{Int},
    neighbors::IncidenceIndex,
    selected::Vector{Int},
    pq::PriorityQueue{Int, Float64},
    root::Int = 1,
)::Vector{Int}
    @assert root >= 1 "root = $root"
    @assert length(selected) >= length(parent)
    @assert isempty(pq) "length(pq) = $(length(pq))"
    sizehint!(pq, length(parent))

    # see CLRS page 572 (2nd edition)
    dist .= typemax(Float64)
    parent .= 0
    dist[root] = typemin(Float64)
    parent[root] = -root
    selected[root] = -1
    pq[root] = dist[root]
    while !isempty(pq)
        u = dequeue!(pq)
        for (v, eidx) in neighbors[u]
            v == u && continue
            if parent[v] <= 0 && edge_weight[eidx] < dist[v]
                dist[v] = edge_weight[eidx]
                pq[v] = dist[v]    # decrease_key!(v)
                parent[v] = -u
                selected[v] = eidx
            end
        end
        parent[u] *= -1
    end
    return parent
end

function _prim_mst_dbg(
    edge_weight::Vector{Float64},
    dist::Vector{Float64},
    parent::Vector{Int},
    neighbors::IncidenceIndex,
    pq::PriorityQueue{Int, Float64},
    root::Int = 1,
)::Vector{Int}
    @assert root >= 1 "root = $root"
    @assert isempty(pq) "length(pq) = $(length(pq))"
    sizehint!(pq, length(parent))

    dist .= typemax(Float64)
    parent .= 0
    dist[root] = typemin(Float64)
    parent[root] = -root
    pq[root] = dist[root]
    while !isempty(pq)
        u = dequeue!(pq)
        println("fin: ", u - 1, " pi = ", -parent[u] - 1, " was ", dist[u])
        parent[u] *= -1
        for (v, eidx) in neighbors[u]
            v == u && continue
            if parent[v] <= 0 && edge_weight[eidx] < dist[v]
                if isinf(dist[v])
                    println("new: ", v - 1, " --> ", edge_weight[eidx])
                else
                    println("upd: ", v - 1, " --> ", edge_weight[eidx], " was ", dist[v])
                end
                dist[v] = edge_weight[eidx]
                pq[v] = dist[v]
                parent[v] = -u
            end
        end
    end
    return parent
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
############################################################################################
function check_correctness(val::AbstractVector{T}) where {T}
    ## Get relevant fields
    s = getfield.(val, 1)
    d = getfield.(val, 2)
    ## Initate container that holds time when state changes
    StateIter = Int64[]
    DurationIter = Int64[]
    ## Compute all state changes
    statechanges = [s[iter] - s[iter - 1] for iter in 2:length(s)]
    durationchanges = [d[iter] - d[iter - 1] for iter in 2:length(d)]
    ## Get all iterations where s changes
    for iter in eachindex(statechanges)
        if statechanges[iter] != 0
            push!(StateIter, iter)
        end
    end
    ## Get all iterations where d changes
    for iter in eachindex(durationchanges)
        if durationchanges[iter] != -1
            push!(DurationIter, iter)
        end
    end
    ## Get all changes that are correct
    changes = [StateIter[iter] == DurationIter[iter] for iter in eachindex(StateIter)]
    ## Return total changes - correct changes (should b 0)
    return length(StateIter) - sum(changes)
end
#Check if HSMM has impossible transitions
function check_correctness(kernel::SemiMarkov, val::AbstractMatrix{T}) where {T}
    Nparticles = size(val, 1)
    return sum([check_correctness(val[iter, :]) for iter in Base.OneTo(Nparticles)])
end

############################################################################################
# run PF and check if correct state trajectories
semimarkov_pf = ParticleFilter(_rng, semimarkov_objective)
_vals, _diag = propose(_rng, semimarkov_pf, semimarkov_objective)

@testset "PF - Ancestors - SemiMarkov" begin
    @test check_correctness(semimarkov_pf.particles.kernel, semimarkov_pf.particles.val) ==
          0
end


text:
@testset "test Number" begin
    one = SUC(ZERO())
    two = SUC(SUC(ZERO()))

    @test repr(ZERO()) == "0"
    @test repr(one) == "1"
    @test repr(two) == "2"
end

@testset "test Add" begin
    one = SUC(ZERO())
    two = SUC(SUC(ZERO()))

    @test Add()(ZERO(), ZERO()) == ZERO()
    @test Add()(ZERO(), one) == one
    @test Add()(one, ZERO()) == one
    @test Add()(one, one) == two
end

@testset "test Sub" begin
    one = SUC(ZERO())
    two = SUC(SUC(ZERO()))

    @test Sub()(ZERO(), ZERO()) == ZERO()
    @test Sub()(one, ZERO()) == one
    @test Sub()(one, one) == ZERO()
    @test Sub()(two, one) == one
end


text:
using Reexport

@reexport module UncertainMathematics
    include("uncertainvalues/elementary_operations_uncertainvalues.jl")
    include("uncertainvalues/elementary_operations_uncertainvalues_special_cases.jl")

    include("uncertainvalues/trig_functions_uncertainvalues.jl")

    include("uncertaindataset/elementary_operations_uncertaindatasets.jl")
    include("uncertainvaluedataset/elementary_operations_uncertainvaluedatasets.jl")
    include("uncertainindexdataset/elementary_operations_uncertainindexdatasets.jl")

end # module


text:
include("concordance.jl")

IOIGAs19=Array{Union{Nothing, String}}(nothing, length(IOIG));
for i in eachindex(IOIG);
    IOIGAs19[i] = IOIGTo19[IOIG[i]]
end

IOHalf = DataFrame(copy(transpose(IOSource[[4:117; 119; 121:126; 128], 3:116])),:auto);
insertcols!(IOHalf ,1, :Industry => IOIGAs19);
IOHalfSplitIndustry = groupby(IOHalf, :Industry);
IOHalf = combine(IOHalfSplitIndustry, valuecols(IOHalfSplitIndustry) .=> sum);

sort!(IOHalf)
IOHalf = select!(IOHalf, Not(:Industry));
IOHalf = transpose(Matrix(IOHalf[:,:]));

IO19 = DataFrame(copy([IOHalf[1:114, :] IOSource[4:117,117:126]]),:auto);
insertcols!(IO19 ,1, :Industry => IOIGAs19);
IO19SplitIndustry = groupby(IO19, :Industry);
IO19 = combine(IO19SplitIndustry, valuecols(IO19SplitIndustry) .=> sum);

FullIO19 = [IO19; [IOHalf[115:122,:] IOSource[[119; 121:126; 128], 117:126]]];

max_stars_repo_name: uq-aibe/Capital

max_issues_repo_name: uq-aibe/Capital

max_forks_repo_name: uq-aibe/Capital


text:
## rearrange.jl  methods the manipulated the data inside an NamedArray

## (c) 2013 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## this does ' as well '
import Base.adjoint
function adjoint(a::NamedArray)
    ndims(a) ≤ 2 || error("Number of dimension must be ≤ 2")
    if ndims(a) == 1
        NamedArray(a.array', (["1"], names(a, 1)), ("'", a.dimnames[1]))
    else
        NamedArray(a.array', reverse(a.dicts), reverse(a.dimnames))
    end
end

import Base.reverse
function reverse(a::NamedArray{T,N}; dims::Int) where {T,N}
    vdicts = Array{OrderedDict}(undef, N)
    n = size(a,dims)+1
    for i=1:N
        dict = copy(a.dicts[i])
        if i==dims
            newnames = reverse(names(dict))
            empty!(dict)
            for (ind,k) in enumerate(newnames)
                dict[k] = ind
            end
        end
        vdicts[i] = dict
    end
    NamedArray(reverse(a.array, dims=dims), tuple(vdicts...), a.dimnames)
end

## circshift automagically works...
## :' automagically works, how is this possible? it is ctranspose!

import Base.permutedims

function permutedims(v::NamedVector)
    NamedArray(reshape(v.array, (1, length(v.array))),
        (["1"], names(v, 1)),
        ("_", v.dimnames[1]))
end

function permutedims(a::NamedArray, perm::Vector{Int})
    dicts = a.dicts[perm]
    dimnames = a.dimnames[perm]
    NamedArray(permutedims(a.array, perm), dicts, dimnames)
end

import Base.transpose
transpose(a::NamedVector) = permutedims(a)
transpose(a::NamedArray) = permutedims(a, [2,1])

import Base.vec
vec(a::NamedArray) = vec(a.array)

import Base.rotl90, Base.rot180, Base.rotr90
rotr90(n::NamedArray) = transpose(reverse(n, dims=1))
rotl90(n::NamedArray) = transpose(reverse(n, dims=2))
rot180(n::NamedArray) = NamedArray(rot180(n.array), tuple([reverse(name) for name in names(n)]...), n.dimnames)

import Combinatorics.nthperm, Combinatorics.nthperm!
import Base.permute!, Base.invpermute!, Base.reverse, Base.reverse!
import Random.shuffle, Random.shuffle!

function nthperm(v::NamedVector, n::Int)
    newnames = nthperm(names(v, 1), n)
    NamedArray(nthperm(v.array,n), (newnames,), v.dimnames)
end
function nthperm!(v::NamedVector, n::Int)
    setnames!(v, nthperm(names(v, 1), n), 1)
    nthperm!(v.array, n)
    return v
end
function permute!(v::NamedVector, perm::AbstractVector)
    setnames!(v, names(v, 1)[perm], 1)
    permute!(v.array, perm)
    return v
end
invpermute!(v::NamedVector, perm::AbstractVector) = permute!(v, invperm(perm))
shuffle(v::NamedVector) = permute!(copy(v), randperm(length(v)))
shuffle!(v::NamedVector) = permute!(v, randperm(length(v)))
reverse(v::NamedVector, start=1, stop=length(v)) = NamedArray(reverse(v.array, start, stop),  (reverse(names(v, 1), start, stop),), v.dimnames)
function reverse!(v::NamedVector, start=1, stop=length(v))
    setnames!(v, reverse(names(v, 1), start, stop), 1)
    reverse!(v.array, start, stop)
    v
end

import Base: _sortslices, _negdims, DimSelector

function my_compute_itspace(A, ::Val{dims}) where {dims}
    negdims = _negdims(ndims(A), dims)
    axs = Iterators.product(ntuple(DimSelector{dims}(A), ndims(A))...)
    vec(permutedims(collect(axs), (dims..., negdims...))), negdims
end

function _sortslices(A::NamedArray, d::Val{dims}; kws...) where dims
    itspace, negdims = my_compute_itspace(A, d)
    vecs = map(its->view(A, its...), itspace)
    p = sortperm(vecs; kws...)
    if ndims(A) == 2 && isa(dims, Integer) && isa(A.array, Array)
        # At the moment, the performance of the generic version is subpar
        # (about 5x slower). Hardcode a fast-path until we're able to
        # optimize this.
        return dims == 1 ? A[p, :] : A[:, p]
    else
        B = similar(A)
        for (x, its) in zip(p, itspace)
            B[its...] = vecs[x]
        end
        if ndims(A) == 2
            nd = negdims[1]
            setnames!(B, names(B, nd)[p], nd)
        else
            @warn("Can't keep dimnames in sortslices() for ndims ≠ 2")
        end
        B
    end
end


text:
using Printf
using Statistics

const LOGDIR = joinpath(@__DIR__, "../../logs/mittelman")
const PLATO_INSTANCES = readlines(joinpath(@__DIR__, "../plato.txt"))
const SOLVERS = ["Clp", "CPLEX", "ECOS", "GLPK", "Gurobi", "Mosek", "Tulip"]
const SUFFIX = Dict{String, String}(
    "Clp" => ".clp",
    "CPLEX" => ".cpx",
    "ECOS" => ".ecos",
    "GLPK" => ".glpk",
    "Gurobi" => ".grb",
    "Mosek" => ".msk",
    "Tulip" => ".tlp"
)

gmean(u, δ) = exp(mean(log.(u .+ δ))) - δ

function parse_res(flog::String, res=Dict{String, Any}())
    res["segfault"] = false
    open(flog, "r") do f
        for ln in eachline(f)
            if occursin("Solver name:", ln)
                res["solver"] = split(ln)[3]
                continue
            end
            if occursin("Solve time:", ln)
                res["time"] = parse(Float64, split(ln)[3])
                continue
            end
            if occursin("Termination status:", ln)
                res["status"] = split(ln)[3]
                continue
            end
            if occursin("Result count:", ln)
                res["rescount"] = split(ln)[3]
                continue
            end
            if occursin("Primal status:", ln)
                res["primal_status"] = split(ln)[3]
                continue
            end
            if occursin("Dual status:", ln)
                res["dual_status"] = split(ln)[3]
                continue
            end
            if occursin("Segmentation fault", ln)
                res["segfault"] = true
            end
        end
    end
    return res
end

function parse_all(instances, logdir)

    # TODO: parse all results in a single dict
    res_all = Dict()

    for finst in instances
        res_all[finst] = Dict()
        for solver in SOLVERS
            res_all[finst][solver] = Dict()
            # Parse log file
            flog = joinpath(logdir, lowercase(solver), finst*SUFFIX[solver])
            r = parse_res(flog)

            # Check result status
            t = get(r, "time", Inf)
            st = get(r, "status", "ERROR")

            if solver == "Gurobi" && st == "OTHER_LIMIT"
                st = "ALMOST_OPTIMAL"
            end
            if solver == "Clp" && st == "OTHER_LIMIT" && t >= 10_000
                st = "TIME_LIMIT"
            end

            res_all[finst][solver]["time"] = (st == "OPTIMAL" || st == "ALMOST_OPTIMAL") ? min(t, 10_000) : 10_000
            res_all[finst][solver]["status"] = st
            res_all[finst][solver]["segfault"] = r["segfault"]
        end
    end

    return res_all
end

function print_table(res)
    # TODO: count number of success for each solver, and geometric mean of total time
    # Header line
    @printf "\\toprule\n"
    @printf "%-16s" "Problem"
    for solver in SOLVERS
        @printf " & %6s" solver
    end
    @printf "\\\\\n"
    @printf "\\midrule\n"

    # Number of instances solved
    @printf "%-16s" "Solved"
    for solver in SOLVERS
        nsolved = 0
        for (finst, r) in res
            st = r[solver]["status"]
            nsolved += (st == "OPTIMAL" || st == "ALMOST_OPTIMAL")
        end
        @printf " & %6d" nsolved
    end
    @printf "\\\\\n"

    # Average time
    @printf "%-16s" "Average"
    for solver in SOLVERS
        times = [res[finst][solver]["time"] for finst in PLATO_INSTANCES]
        μ = gmean(times, 10.0)
        @printf " & %6.1f" μ
    end
    @printf "\\\\\n"

    # Individual times
    @printf "\\midrule\n"
    for finst in sort(PLATO_INSTANCES)
        @printf "%-16s" replace(finst[1:end-4], "_" => "\\_")
        for solver in SOLVERS
            @printf " & "
            # Parse log file
            flog = joinpath(LOGDIR, lowercase(solver), finst*SUFFIX[solver])
            r = parse_res(flog)

            # Check result status
            t = res[finst][solver]["time"]
            st = res[finst][solver]["status"]
            segfault = res[finst][solver]["segfault"]

            if segfault
                @printf "%s" "\\texttt{seg}"
            elseif t < 10_000 && st == "OPTIMAL"
                # Problem was solved within the time limit
                @printf "%6.1f" t
            elseif t < 10_000 && st == "ALMOST_OPTIMAL"
                @printf "%6.1f\\texttt{r}" t
            elseif st == "TIME_LIMIT"
                # Time limit reached
                @printf "%6s" "\\texttt{t}"
            else
                # Other error
                @printf "%6s" "\\texttt{f}"
            end
        end

        @printf "\\\\\n"
    end
    @printf "\\bottomrule\n"
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    res = parse_all(PLATO_INSTANCES, LOGDIR)
    print_table(res)
    return nothing
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
module RMB

importall NLsolve
import Base.rand
import Distributions: wsample
import Calculus

export solvedp!, u!, pi!
export InventoryModel

include("stochastic-matrices.jl")
include("rust-models.jl")
include("inv_model.jl")

end #module


text:
using LinearAlgebra: Diagonal

# Loss function should return a pair of (e, ∂e/∂ŷ).
export
    mse,
    softmax, log_softmax,
    logit_cross_entropy,
    logit_binary_cross_entropy

"Mean squared error loss"
@inline mse(ŷ, y) = (0.5 * sum((y - ŷ) .^ 2), ŷ - y)

#=====

**Softmax** turns (multinomial) logit into probability:

    softmax(x)_i = exp(x_i) / Σ exp(x_i)

**Cross entropy** evaluates the expected coding length of data following distribution p when a wrong
distribution q is assumed:

    H(p, q) = -E_p[log(q)]
            = -Σ p log(q)
            = H(p) + D_KL(p || q)
            = -Σ p log(p) + Σ p log(p/q)

**Binary cross entropy** is the special case of cross entropy for 2-class classification. It is
applied on a single sigmoid output rather than a vector of size 2 as in the softmax case.

To mitigate numerical instabilities, computing by a composed function (prefixed with "logit_") is
usually preferred:

    logit_cross_entropy(ŷ, y) := cross_entropy(softmax(ŷ), y)
    logit_binary_cross_entropy(ŷ, y) := binary_cross_entropy(logistic(ŷ), y)

Softmax is invariant under translation by the same value in each coordinate. The standard logistic
function is a special case of softmax where one variable is fixed at 0.

=====#

function softmax(x::AbstractMatrix{Float64})::Matrix{Float64}
    max_x = maximum(x, dims=1)
    exp_x = exp.(x .- max_x)
    exp_x ./ sum(exp_x, dims=1)
end

function log_softmax(x::AbstractMatrix{Float64})::Matrix{Float64}
    x_ = x .- maximum(x, dims=1)
    x_ .- log.(sum(exp.(x_), dims=1))
end

#=

Derivation:

    p_i = softmax(y)_i  ⟹
    ∂p_i/∂y_j = [i = j] p_i - p_i p_j

    ∂e/∂y_j = Σ_i ∂e/∂p_i ⋅ ∂p_i/∂y_j

Or in the form of matrix calculus,

    ∂e/∂y = ∂e/∂p ⋅ ∂p/∂y

And

    e = -sum(y .* log(p)) / n  ⟹
    ∂e/∂p = -y / n * (1/p)

Also note that $∂p/∂y ∈ ℜ_{d × d}$, so if mini-batch is used it is impossible to be expressed in
matrix form. We directly write `∂e/∂p ⋅ ∂p/∂y` of one sample as

      (∂e/∂p) ⋅ (diag(p) - p ⋅ p')
    = ∇p .* p - sum(∇p .* p) * p    # a ℜ_d vector

The mini-batch case can be written as

    ... = ∇p .* p - sum(∇p .* p, dims=1) .* p    # a ℜ_{d × n} matrix

=#
function logit_cross_entropy(
    ŷ::AbstractMatrix{Float64},
    y::AbstractMatrix{Float64},
)::Tuple{Float64, Matrix{Float64}}
    n_batches = size(ŷ, 2)
    log_p = log_softmax(ŷ)
    e = -sum(y .* log_p) / n_batches

    p = softmax(ŷ)
    ∇p = @. -y / n_batches * inv(p)
    t = ∇p .* p
    ∇y = t - sum(t, dims=1) .* p

    (e, ∇y)
end

#=
    logit_binary_cross_entropy(ŷ, 0)
        = -log(1 - logistic(ŷ))
        = -log(1 / (1 + exp(ŷ)))
        = -(0 - log(1 + exp(ŷ)))
        = -(ŷ - log(1 + exp(ŷ))) + ŷ
        = -log_logistic(ŷ) + ŷ
=#
function logit_binary_cross_entropy(
    ŷ::AbstractMatrix{Float64},
    y::AbstractMatrix{Bool},
)::Tuple{Float64, Matrix{Float64}}
    @assert size(ŷ, 1) == 1
    t = -log_logistic.(ŷ)
    ny = .!y
    t[ny] .+= ŷ[ny]
    e = sum(t) / size(ŷ, 2)

    ∂t_over_∂y = -∇log_logistic.(ŷ, -t)
    ∂t_over_∂y[ny] .+= 1.0
    ∇y = ∂t_over_∂y

    (e, ∇y)
end

#=

Addendum
--------

The derivation of the gradient of softmax + cross-entropy can be found in:

- https://stats.stackexchange.com/questions/235528/backpropagation-with-softmax-cross-entropy
- https://charlee.li/how-to-compute-the-derivative-of-softmax-and-cross-entropy/
- http://denizyuret.github.io/Knet.jl/stable/softmax/#Softmax-1

Copy the LaTeX here:

$$
\begin{aligned}
\frac{\partial}{\partial z_k}\text{CE} &= \frac{\partial}{\partial z_k}\sum_{j=1}^n
\big(- y_j \log \sigma(z_j) \big) \\  &= -\sum_{j=1}^n y_j \frac{\partial}{\partial z_k} \log
\sigma(z_j) && \cdots \text{addition rule, } -y_j \text{ is constant}\\  &= -\sum_{j=1}^n y_j
\frac{1}{\sigma(z_j)} \cdot \frac{\partial}{\partial z_k} \sigma(z_j) && \cdots \text{chain rule}\\
&= -y_k \cdot \frac{\sigma(z_k)(1-\sigma(k))}{\sigma(z_k)} –  \sum_{j \neq k} y_j \cdot
\frac{-\sigma(z_j)\sigma(z_k)}{\sigma(z_j)} && \cdots \text{consier both }j = k \text{ and } j \neq
k \\  &= -y_k \cdot (1-\sigma(z_k)) + \sum_{j \neq k} y_j \sigma(z_k) \\  &= -y_k + y_k \sigma(z_k)
+  \sum_{j \neq k} y_j \sigma(z_k) \\  &= -y_k + \sigma(z_k) \sum_j y_j.
\end{aligned}
$$

=#

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
__precompile__()

module AlmostBandedMatrices
using BandedMatrices, ToeplitzMatrices, Base, Compat, FillArrays
if VERSION ≥ v"0.7-"
    using LinearAlgebra, SparseArrays, Random
    using LinearAlgebra.LAPACK
    import LinearAlgebra: BlasInt,
                        BlasReal,
                        BlasFloat,
                        BlasComplex,
                        axpy!,
                        A_mul_B!,
                        Ac_mul_B,
                        Ac_mul_B!,
                        A_mul_Bc,
                        A_mul_Bc!,
                        Ac_mul_Bc,
                        Ac_mul_Bc!,
                        At_mul_B,
                        At_mul_B!,
                        A_mul_Bt,
                        A_mul_Bt!,
                        At_mul_Bt,
                        At_mul_Bt!,
                        copy_oftype,
                        checksquare,
                        adjoint,
                        transpose
   import LinearAlgebra.BLAS: libblas
   import LinearAlgebra.LAPACK: liblapack
   import LinearAlgebra: cholfact, cholfact!, norm, diag, eigvals!, eigvals,
                qr, qrfact, axpy!, ldiv!
   import SparseArrays: sparse

   const lufact = LinearAlgebra.lu # TODO: Remove once 0.6 is dropped
else
    import Base.LinAlg: BlasInt,
                        BlasReal,
                        BlasFloat,
                        BlasComplex,
                        axpy!,
                        A_mul_B!,
                        Ac_mul_B,
                        Ac_mul_B!,
                        A_mul_Bc,
                        A_mul_Bc!,
                        Ac_mul_Bc,
                        Ac_mul_Bc!,
                        At_mul_B,
                        At_mul_B!,
                        A_mul_Bt,
                        A_mul_Bt!,
                        At_mul_Bt,
                        At_mul_Bt!,
                        A_ldiv_B!,
                        At_ldiv_B!,
                        Ac_ldiv_B!,
                        copy_oftype,
                        checksquare
   using Base.LAPACK
   import Base.BLAS: libblas
   import Base.LAPACK: liblapack
   import Base: lufact, cholfact, cholfact!, norm, diag, eigvals!, eigvals,
                At_mul_B, Ac_mul_B, A_mul_B!, qr, qrfact, axpy!
   import Base: sparse, rank

   rmul!(A::AbstractArray, b::Number) = scale!(A, b)
   lmul!(a::Number, B::AbstractArray) = scale!(a, B)
   parentindices(A) = parentindexes(A)
end

import Base: getindex, setindex!, *, +, -, ==, <, <=, >, >=, /, ^, \,
             transpose, showerror, reindex, checkbounds, @propagate_inbounds

import Base: convert, view, indices, unsafe_indices, indices1, first,
             last, size, length, unsafe_length, start, next, done, step,
             to_indices, to_index, indices, show, fill!, copy!, promote_op

import FillArrays: AbstractFill

export AlmostBandedMatrix, SemiSeparableMatrix

# include("blas.jl")
# include("lapack.jl")

include("LowRankMatrix.jl")
include("SemiSeparableMatrix.jl")
include("AlmostBandedMatrix.jl")

end #module


text:
module Data
include("utils.jl")
using .Zip

const _base_kaggle_error =
"""
You MUST have kaggle installed to download the dataset or manually download and extract it yourself"

To manually downlead it yourself:
Navigate to the place you downleaded this project and create a directory called data.
Download the BirdClef2022 dataset from kaggle.com and extract the birdclef-2022 directory into the
directory you just created.
This project should now have the file structure:

data/
    birdclef-2022/
        ...
src/
    ...
...

Using Kaggle:
Download and install kaggle using your python's pip.
Once you have done that, go to https://www.kaggle.com/docs/api to get your credentials.
Put the credentials at
Linux/Mac:  $(get(ENV,"HOME","\$HOME"))/.kaggle/kaggle.json
Windows  : C:\\Users<Windows-username>.kaggle\\kaggle.json
"""

kaggle_error(pythoncmd) = isnothing(pythoncmd) ? _base_kaggle_error : "$_base_kaggle_error\nTo download kaggle, run this in your terminal:\n\t\$ $pythoncmd -m pip install kaggle.\n"

function foldersize(folder)
    total_size = 0
    function appendwork!(work,folder)
        for (root,dirs,files) in walkdir(folder)

            append!(work,map(dirs) do dir
                joinpath(root,dir)
            end)
            files = map(files) do file
                joinpath(root,file)
            end
            for file in files
                total_size+=filesize(file)
            end
        end
    end
    work =[folder]
    while length(work)>0
        current_job = pop!(work)
        appendwork!(work,current_job)
    end
    total_size
end

function validate_directory(directories...)
    is_valid,h = false,false
    for directory in directories
        if !isdir(directory)
            if !isfile(directory)
                mkpath(directory)
                is_valid=false
                h=true
            else
                throw(ErrorException("$directory must not be a file"))
            end
        else
            if !h
                is_valid=true
            end
        end
    end
    is_valid
end

function has_dataset(path::String,dataset::String)
    data_path = joinpath(path,dataset)
    validate_directory(data_path) && foldersize(data_path)>=6600000000
end

function find_python_kaggle(python)
    try
    run(`$python -m kaggle --v`)
    return "$python -m kaggle"
    catch e
        throw(ErrorException(kaggle_error(python)))
    end
end
function find_kaggle()
    kaggle = Sys.which("kaggle")
    if !isnothing(kaggle)
        return kaggle
    end
    python3 = Sys.which("python3")
    if !isnothing(python3)
        return find_python_kaggle(python3)
    end
    python = Sys.which("python")
    if !isnothing(python)
        return find_python_kaggle(python)
    end
    throw(ErrorException(kaggle_error(nothing)))
end

function download_dataset(path::String,dataset::String)
    @info "Downloading dataset..." path dataset
    @info "Searching for kaggle..."
    kaggle = find_kaggle()
    cmd=`$kaggle competitions download -c $dataset -p $path`
    @info "Starting data installation using" cmd
    run(cmd)
    zipped,unzipped = joinpath(path,"$dataset.zip"),joinpath(path,dataset)
    @info "extracting files...." zipped unzipped
    extract(zipped,unzipped)
    @info "cleaning up..."
    rm(zipped,force=true)
    @info "Done downloading dataset" dataset
end

export download_dataset,has_dataset,validate_directory,foldersize

end

max_stars_repo_name: BirdClef/birdclef

max_issues_repo_name: BirdClef/birdclef

max_forks_repo_name: BirdClef/birdclef


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

@with_kw mutable struct ChabocheDriverState <: AbstractMaterialState
    time :: Float64 = zero(Float64)
    strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
end

@with_kw struct ChabocheParameterState <: AbstractMaterialState
    E :: Float64 = 0.0
    nu :: Float64 = 0.0
    R0 :: Float64 = 0.0
    Kn :: Float64 = 0.0
    nn :: Float64 = 0.0
    C1 :: Float64 = 0.0
    D1 :: Float64 = 0.0
    C2 :: Float64 = 0.0
    D2 :: Float64 = 0.0
    Q :: Float64 = 0.0
    b :: Float64 = 0.0
end

@with_kw struct ChabocheVariableState <: AbstractMaterialState
    stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    X1 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    X2 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
    cumeq :: Float64 = zero(Float64)
    R :: Float64 = zero(Float64)
    jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64})
end

@with_kw mutable struct Chaboche <: AbstractMaterial
    drivers :: ChabocheDriverState = ChabocheDriverState()
    ddrivers :: ChabocheDriverState = ChabocheDriverState()
    variables :: ChabocheVariableState = ChabocheVariableState()
    variables_new :: ChabocheVariableState = ChabocheVariableState()
    parameters :: ChabocheParameterState = ChabocheParameterState()
    dparameters :: ChabocheParameterState = ChabocheParameterState()
end

function integrate_material!(material::Chaboche)
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R, jacobian = v

    jacobian = isotropic_elasticity_tensor(lambda, mu)

    stress += dcontract(jacobian, dstrain)
    seff = stress - X1 - X2
    seff_dev = dev(seff)
    f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
    if f > 0.0
        g! = create_nonlinear_system_of_equations(material)
        x0 = [tovoigt(stress); R; tovoigt(X1); tovoigt(X2)]
        F = similar(x0)
        res = nlsolve(g!, x0; autodiff = :forward)
        x = res.zero
        res.f_converged || error("Nonlinear system of equations did not converge!")

        stress = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[1:6])
        R = x[7]
        X1 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[8:13])
        X2 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[14:19])
        seff = stress - X1 - X2
        seff_dev = dev(seff)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        plastic_strain += dp*n
        cumeq += dp
        # Compute Jacobian
        function residuals(x)
            F = similar(x)
            g!(F, x)
            return F
        end
        drdx = ForwardDiff.jacobian(residuals, x)
        drde = zeros((length(x),6))
        drde[1:6, 1:6] = -tovoigt(jacobian)
        jacobian = fromvoigt(SymmetricTensor{4,3}, (drdx\drde)[1:6, 1:6])
    end
    variables_new = ChabocheVariableState(stress = stress,
                                          X1 = X1,
                                          X2 = X2,
                                          R = R,
                                          plastic_strain = plastic_strain,
                                          cumeq = cumeq,
                                          jacobian = jacobian)
    material.variables_new = variables_new
end

function create_nonlinear_system_of_equations(material::Chaboche)
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    mu = E/(2.0*(1.0+nu))
    lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R = v

    function g!(F, x::Vector{T}) where {T} # System of non-linear equations
        jacobian = isotropic_elasticity_tensor(lambda, mu)
        stress_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[1:6])
        R_ = x[7]
        X1_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[8:13])
        X2_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[14:19])

        seff = stress_ - X1_ - X2_
        seff_dev = dev(seff)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R_)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)
        dstrain_plastic = dp*n
        tovoigt!(view(F, 1:6), stress - stress_ + dcontract(jacobian, dstrain - dstrain_plastic))
        F[7] = R - R_ + b*(Q-R_)*dp
        if isapprox(C1, 0.0)
            tovoigt!(view(F,8:13),X1 - X1_)
        else
            tovoigt!(view(F,8:13), X1 - X1_ + 2.0/3.0*C1*dp*(n - 1.5*D1/C1*X1_))
        end
        if isapprox(C2, 0.0)
            tovoigt!(view(F,14:19), X2 - X2_)
        else
            tovoigt!(view(F, 14:19), X2 - X2_ + 2.0/3.0*C2*dp*(n - 1.5*D2/C2*X2_))
        end
    end
    return g!
end


text:
""" Cell type specificity plots (ODE model), Figure 2. """

using Gadfly;
using Statistics;
using gcSolver;
using DataFrames;
using ForwardDiff;
using StatsFuns;
import LinearAlgebra: diag;

gdf = Gadfly;

# Plot of dose response curves
function doseResPlot2(ligandName, cellType, date, unkVec, alphaCov = true)
    responseDF = gcSolver.importData(true)
    DateFrame = gcSolver.importConvFrame()
    tps = [0.5, 1, 2, 4] .* 60
    doseVec = unique(responseDF, "Dose")
    doseVec = doseVec[!, :Dose]
    sigma = gcSolver.getSigma(cellType)
    if alphaCov
        predictDF = DataFrame(Dose = Float64[], tps = Float64[], sensitivity = Float64[])
    else
        predictDF = DataFrame(Dose = Float64[], tps = Float64[], Variance = Float64[])
    end

    filtFrame = filter(row -> row["Ligand"] .== ligandName, responseDF)
    filter!(row -> row["Cell"] .== cellType, filtFrame)
    filter!(row -> string(row["Date"]) .== date, filtFrame)

    if alphaCov
        realDataDF = filtFrame[!, [:Dose, :Time, :Mean, :alphStatCov]]
        realDataDF = groupby(realDataDF, [:Time, :Dose])
        realDataDF = combine(realDataDF, :Mean => mean, :alphStatCov => mean)
    else
        realDataDF = filtFrame[!, [:Dose, :Time, :Mean, :Variance]]
        realDataDF = groupby(realDataDF, [:Time, :Dose])
        realDataDF = combine(realDataDF, :Mean => mean, :Variance => mean)
    end

    for (i, dose) in enumerate(doseVec)
        #check if ligand name is IL2

        if ligandName == "IL15"
            #put into second slot
            doseLevel = [0, dose, 0]
        else
            #put ILdose into first slot
            doseLevel = [dose, 0, 0]
        end

        idxx = findfirst(responseDF.Cell .== cellType)
        recAbunds = 10.0 .^ Vector{Float64}(responseDF[idxx, [:IL2Ra, :IL2Rb, :gc, :IL15Ra, :IL7Ra]])
        iterParams = gcSolver.fitParams(doseLevel, unkVec, recAbunds, cellType)

        iterParams = gcSolver.mutAffAdjust(iterParams, responseDF[findfirst(responseDF.Ligand .== ligandName), [:IL2RaKD, :IL2RBGKD]])

        if alphaCov
            JacResults = runAlphJac(tps, iterParams, sigma, recAbunds)
            for indx = 1:length(tps)
                #use dataframe and push row into it - enter data into data frame
                push!(predictDF, (dose, tps[indx] / 60, JacResults[indx]))
            end
        else
            VarResults = runCkineVarProp(tps, iterParams, sigma, recAbunds)
            for indx = 1:length(tps)
                #use dataframe and push row into it - enter data into data frame
                push!(predictDF, (dose, tps[indx] / 60, VarResults[indx]))
            end
        end
    end

    if alphaCov

        predictDF.sensitivity .*= filter(row -> row.Date ∈ [date], DateFrame).Conv .^ 2
        realDataDF.alphStatCov_mean = abs.(realDataDF.alphStatCov_mean)

        pl1 = gdf.plot(
            layer(predictDF, x = :Dose, y = :sensitivity, color = :tps, Geom.line),
            layer(realDataDF, x = :Dose, y = :alphStatCov_mean, color = :Time, Geom.point),
            Scale.x_log10,
            Scale.y_log10,
            Guide.title(string(cellType, " Response to ", ligandName)),
            Guide.xlabel("Dose"),
            Guide.ylabel("Sensitivity"),
            Scale.color_discrete(),
            Guide.colorkey(title = "Time (hr)", labels = ["4", "2", "1", "0.5"]),
            Coord.cartesian(ymin = 1, ymax = 10),
        )
    else
        predictDF.Variance .*= filter(row -> row.Date ∈ [date], DateFrame).Conv .^ 2
        pl1 = gdf.plot(
            layer(predictDF, x = :Dose, y = :Variance, color = :tps, Geom.line),
            layer(realDataDF, x = :Dose, y = :Variance_mean, color = :Time, Geom.point),
            Scale.x_log10,
            Guide.title(string(cellType, " Response to ", ligandName)),
            Guide.xlabel("Dose"),
            Guide.ylabel("Variance"),
            Scale.color_discrete(),
            Guide.colorkey(title = "Time (hr)", labels = ["4", "2", "1", "0.5"]),
            Scale.y_log10,
            Coord.cartesian(ymin = 1, ymax = 10),
        )
    end
    return pl1
end

function runAlphJac(tps::Vector, params::Vector, sigma, recAbund)
    gcSolver.checkInputs(tps, params)

    # Sigma is the covariance matrix of the input parameters
    function jacF(x)
        pp = vcat(params[1:24], x[1], params[26:end])
        return runCkine(tps, pp, pSTAT5 = true)
    end

    jac = zeros(1, length(tps))
    ForwardDiff.jacobian!(jac, jacF, [params[25]])
    jac ./= (recAbund[1] / params[25]) #maybe
    return diag(transpose(jac) * sigma[1, 1] * jac)
end

"""Use this if you want to change the parameters here and not input any in the command line"""
function figureJ2()
    fitVec = gcSolver.importFit()
    fitVec = convert(Vector{Float64}, fitVec[!, :Fit])
    fitVec = softplus.(fitVec)

    p1 = doseResPlot2("IL2", "Treg", "4/18/2019", fitVec, true)
    p2 = doseResPlot2("IL2", "Thelper", "4/18/2019", fitVec, true)
    p3 = doseResPlot2("IL2", "NK", "3/15/2019", fitVec, false)
    p4 = doseResPlot2("IL2", "CD8", "3/15/2019", fitVec, false)
    p5 = doseResPlot2("N88D C-term", "Treg", "12/5/2019", fitVec, true)
    p6 = doseResPlot2("N88D C-term", "Thelper", "12/5/2019", fitVec, true)
    p7 = doseResPlot2("N88D C-term", "NK", "12/5/2019", fitVec, false)
    p8 = doseResPlot2("N88D C-term", "CD8", "12/5/2019", fitVec, false)
    p9 = doseResPlot2("WT N-term", "Treg", "11/8/2019", fitVec, true)
    p10 = doseResPlot2("WT N-term", "Thelper", "11/8/2019", fitVec, true)
    p11 = doseResPlot2("WT N-term", "NK", "11/8/2019", fitVec, false)
    p12 = doseResPlot2("WT N-term", "CD8", "11/8/2019", fitVec, false)
    p13 = doseResPlot2("V91K C-term", "Treg", "11/27/2019", fitVec, false)
    p14 = doseResPlot2("V91K C-term", "Thelper", "11/27/2019", fitVec, false)
    p15 = doseResPlot2("V91K C-term", "NK", "11/27/2019", fitVec, false)
    p16 = doseResPlot2("V91K C-term", "CD8", "11/27/2019", fitVec, false)
    p17 = doseResPlot2("R38Q N-term", "Treg", "12/5/2019", fitVec, false)
    p18 = doseResPlot2("R38Q N-term", "Thelper", "12/5/2019", fitVec, false)
    p19 = doseResPlot2("R38Q N-term", "NK", "12/5/2019", fitVec, false)
    p20 = doseResPlot2("R38Q N-term", "CD8", "12/5/2019", fitVec, false)
    p21 = doseResPlot2("F42Q N-Term", "Treg", "12/5/2019", fitVec, false)
    p22 = doseResPlot2("F42Q N-Term", "Thelper", "12/5/2019", fitVec, false)
    p23 = doseResPlot2("F42Q N-Term", "NK", "12/5/2019", fitVec, false)
    p24 = doseResPlot2("F42Q N-Term", "CD8", "12/5/2019", fitVec, false)

    draw(
        SVG("figureJ2.svg", 4000px, 2400px),
        gridstack([p1 p2 p3 p4; p5 p6 p7 p8; p9 p10 p11 p12; p13 p14 p15 p16; p17 p18 p19 p20; p21 p22 p23 p24]),
    )
end


text:
# StackedLSTM.jl
#
# Layers for stacked LSTM

mutable struct StackedLSTMCell{A}
	chain::A
end

"""
    HiddenRecur(cell)

`HiddenRecur` takes a recurrent cell and makes it stateful, managing the hidden state
in the background. As opposed to `Recur`, it returns the hidden state from the cell when
called, and not the cell output.

See also: [`Recur`](@ref)
"""
mutable struct HiddenRecur{T}
  cell::T
  init
  state
end

HiddenRecur(m, h = Flux.hidden(m)) = HiddenRecur(m, h, h)

function (m::HiddenRecur)(xs...)
  h, y = m.cell(m.state, xs...)
  m.state = h
  return h[1]
end

Flux.@functor HiddenRecur cell, init

Base.show(io::IO, m::HiddenRecur) = print(io, "HiddenRecur(", m.cell, ")")

Flux.trainable(m::StackedLSTMCell) = (m.chain)
Flux.reset!(m::HiddenRecur) = (m.state = m.init)

"""
	StackedLSTM(in, out, hiddensize, layers)

Stacked LSTM network. Feeds the data through a chain of LSTM layers, where the hidden state
of the previous layer gets fed to the next one. The first layer corresponds to
`LSTM(in, hiddensize)`, the hidden layers to `LSTM(hiddensize, hiddensize)`, and the final
layer to `LSTM(hiddensize, out)`. Takes the keyword argument `init` for the initialization
of the layers.
"""
function StackedLSTM(in::Integer, out::Integer, hiddensize::Integer, layers::Integer;
			init = Flux.glorot_uniform)
	if layers == 1
		chain = Chain(LSTM(in, out; init = init))
	elseif layers == 2
		chain = Chain(HiddenRecur(Flux.LSTMCell(in, hiddensize; init = init)),
					LSTM(hiddensize, out; init = init))
	else
		chain_vec=[HiddenRecur(Flux.LSTMCell(in, hiddensize; init = init))]
		for i=1:layers-2
			push!(chain_vec, HiddenRecur(Flux.LSTMCell(hiddensize, hiddensize; init = init)))
		end
		chain = Chain(chain_vec..., LSTM(hiddensize, out; init = init))
	end
	return StackedLSTMCell(chain)
end

function (m::StackedLSTMCell)(x)
	return m.chain(x)
end

function Base.show(io::IO, l::StackedLSTMCell)
	if length(l.chain)==1
		print(io, l.chain[1])
	else
		print(io, "StackedLSTMCell(", size(l.chain[1].cell.Wi, 2), ", ", size(l.chain[end].cell.Wi, 1)÷4)
		print(io, ", ", size(l.chain[end].cell.Wi, 2))
		print(io, ", ", length(l.chain), ")")
	end
end

max_stars_repo_name: AndyTGardner/FluxArchitectures

max_issues_repo_name: AndyTGardner/FluxArchitectures

max_forks_repo_name: AndyTGardner/FluxArchitectures