text:
include("filtering.jl");
gr(size=(600,600), tick_orientation = :out, grid = false,
    linecolor = :black,
    markerstrokecolor = :black,
    thickness_scaling = 2,
    markersize = 6)
################################ Adjust Streaks table ##################################
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
filter!(r->r.Treatment in list &&
    r.Trial < 51 &&
    r.MouseID != "pc7",
    streaks)

streaks[streaks.Treatment .== "PreVehicle",:Treatment] .= "Control"
streaks[streaks.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in streaks[streaks.Treatment .== "Saline",:Stim]]
streaks[streaks.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
streaks[streaks.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in streaks[streaks.Treatment .== "SB242084_opt",:Stim]]
streaks[!,:ROILeavingTime] = streaks.Stop_trial .- streaks.Stop_poking
streaks[!,:Treatment] = categorical(streaks.Treatment, ordered = false)
levels!(streaks.Treatment,["Control",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Citalopram",
    "Optogenetic",
    "Methysergide",
    "SB242084_opt"
    ])
levels(streaks.Treatment)
################################ Adjust Pokes table ##################################
filter!(r->r.Treatment in list &&
    r.Trial < 51 &&
    r.MouseID != "pc7",
    pokes)
pokes[pokes.Treatment .== "PreVehicle",:Treatment] .= "Control"
pokes[pokes.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in pokes[pokes.Treatment .== "Saline",:Stim]]
pokes[pokes.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
pokes[pokes.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in pokes[pokes.Treatment .== "SB242084_opt",:Stim]]

filt_1 = filter(r->r.Treatment == "Control",pokes)
gd1 = groupby(filt_1,[:MouseID,:TimeFromLeaving,:Protocol])
df2 = combine(gd1, :CumRewTrial => mean => :CumRewTrial)
gd2 = groupby(df2,[:Protocol,:TimeFromLeaving])
df3 = combine(gd2, :CumRewTrial => mean, :CumRewTrial => sem)
Protocol_colors!(df3)

@df df3 scatter(:TimeFromLeaving,:CumRewTrial_mean, group = :Protocol, xlims = (0,30), markersize = 4, color = :color)

################################ Fig2 scatter  controls ##################################
Pstate= Prew(1:20)
open_html_table(Pstate)
Pstate[!,:Color] = [get(protocol_colors,x,:grey) for x in Pstate.Protocol]
@df Pstate plot(:Poke,:Prew, group = :Protocol, linecolor = :Color, legend = false)
@df Pstate scatter!(:Poke,:Prew, group = :Protocol, color = :Color)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","Distributions.pdf"))

gd1 = groupby(streaks,[:Protocol,:MouseID,:Phase,:Treatment])
df1 = combine(gd1, :Num_Rewards => mean => :Num_Rewards,
    :Num_pokes => mean => :Num_pokes,
    :AfterLast => mean => :AfterLast,
    :Leaving_NextPrew => mean => :Leaving_NextPrew,
    :AverageRewRate => mean => :AverageRewRate)

gd2 = groupby(df1,:Treatment)
df2 = combine(gd2, :AverageRewRate => mean, :AverageRewRate => sem,
    :Leaving_NextPrew => mean, :Leaving_NextPrew => sem)
filt_2 = filter(r->r.Treatment == "Control",df2)
res = DataFrame(Condition = ["Average", "At Leaving"],
    Mean = [filt_2[1,:AverageRewRate_mean], filt_2[1,:Leaving_NextPrew_mean]],
    Sem = [filt_2[1,:AverageRewRate_sem], filt_2[1,:Leaving_NextPrew_sem]])
@df res scatter(:Condition,:Mean, yerror = :Sem, color = :grey, xlims = (0.25,1.75), label = false,ylims = (0.1,0.3), yticks = 0.0:0.1:0.3)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","RateAtLeaving.pdf"))

filt_2 = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control"],df2)
Drug_colors!(filt_2)
res = combine(groupby(filt_2, :Treatment)) do dd
    DataFrame(Condition = ["Average", "At Leaving"],
        Mean = [dd[1,:AverageRewRate_mean], dd[1,:Leaving_NextPrew_mean]],
        Sem = [dd[1,:AverageRewRate_sem], dd[1,:Leaving_NextPrew_sem]])
    end
Drug_colors!(res)
@df res scatter(:Condition,:Mean, yerror = :Sem, color = :color, xlims = (0.25,1.75), fillalpha=0.5, label = false)

gd3 = groupby(df1,[:Protocol,:Treatment])
df3 = combine(gd3, :Num_Rewards => mean, :Num_Rewards => sem,
    :Num_pokes => mean, :Num_pokes => sem,
    :AfterLast => mean, :AfterLast => sem,
    :Leaving_NextPrew => mean, :Leaving_NextPrew => sem)

filt_3 = filter(r->r.Treatment == "Control",df3)
sort!(filt_3, :Protocol)
filt_3[!,:Color] = [get(protocol_colors,x,:grey) for x in filt_3.Protocol]
@df filt_3 scatter(:Protocol, :Num_Rewards_mean, yerror = :Num_Rewards_sem, label = false, color = :Color, xlims = (0.25,2.75), ylims = (1,5), yticks = 0:0.5:10)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","NumRewards.pdf"))
@df filt_3 scatter(:Protocol, :Num_pokes_mean, yerror = :Num_pokes_sem, label = false, color = :Color, xlims = (0.25,2.75), ylims = (11,15))
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","NumPokes.pdf"))
@df filt_3 scatter(:Protocol, :AfterLast_mean, yerror = :AfterLast_sem, label = false, color = :Color, xlims = (0.25,2.75), ylims = (5,7))
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","AfterLast.pdf"))
@df filt_3 scatter(:Protocol, :Leaving_NextPrew_mean, yerror = :Leaving_NextPrew_sem,label = false, color = :Color, xlims = (0.25,2.75), ylims = (0.1,0.3), yticks = 0.0:0.1:0.3)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig2","Prew.pdf"))

testdf = filter(r->r.Treatment == "Control",streaks)
testdf[!,:NumericalProt] = parse.(Float64, testdf.Protocol)
NR_m1 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + (1|MouseID)),testdf))
NR_m2 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + NumericalProt + (1|MouseID)),testdf))
Likelyhood_Ratio_test(NR_m1,NR_m2)
NP_m1 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + (1|MouseID)),testdf))
NP_m2 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt + (1|MouseID)),testdf))
Likelyhood_Ratio_test(NP_m1,NP_m2)
AL_m1 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + (1|MouseID)),testdf))
AL_m2 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt + (1|MouseID)),testdf))
Likelyhood_Ratio_test(AL_m1,AL_m2)
PR_m1 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + (1|MouseID)),testdf))
PR_m2 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt + (1|MouseID)),testdf))
Likelyhood_Ratio_test(PR_m1,PR_m2)

# fm4 = fit(MixedModel,@formula(AfterLast ~ 1 + (1|MouseID)),age_df,Poisson())
################################ Fig3 selective protocol effect ##################################
df1 = combine(groupby(streaks,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:ROI_Leaving_Time)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end

df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(filter!(r->r.Phase in ["SB242084","Altanserin"],df2),showmice = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_Wilcoxon_ROILeaving.pdf"))

df2 = combine(groupby(df1,[:Protocol,:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(filter!(r->r.Phase in ["SB242084","Altanserin"],df2),showmice = false)
################################ Fig3 scatter  selective ##################################
gd1 = groupby(streaks,[:Protocol,:MouseID,:Phase,:Treatment])
df1 = combine(gd1, :Num_Rewards => mean => :Num_Rewards,
    :Num_pokes => mean => :Num_pokes,
    :AfterLast => mean => :AfterLast,
    :Leaving_NextPrew => mean => :Leaving_NextPrew,
    :AverageRewRate => mean => :AverageRewRate)

gd2 = groupby(df1,[:Protocol,:Treatment])

df2 = combine(gd2, [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> mean,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> sem,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> median,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> CIq
    )

dp = filter!(r->r.Treatment in ["SB242084","Altanserin", "Control"],df2)
Drug_colors!(dp)
@df dp scatter(:Protocol, :Num_Rewards_mean, yerror = :Num_Rewards_sem, group = :Treatment, color = :color, xflip = true, yticks = 1:0.5:6, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Num_Rewards_median, yerror = :Num_Rewards_CIq, group = :Treatment, color = :color, xflip = true, yticks = 1:0.5:6, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_Prot_NumRewards.pdf"))

@df dp scatter(:Protocol, :Num_pokes_mean, yerror = :Num_pokes_sem, group = :Treatment, color = :color, xflip = true, yticks = 8:1:17, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Num_pokes_median, yerror = :Num_pokes_CIq, group = :Treatment, color = :color, xflip = true, yticks = 8:1:17, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_Prot_NumPokes.pdf"))

@df dp scatter(:Protocol, :AfterLast_mean, yerror = :AfterLast_sem, group = :Treatment, color = :color, xflip = true, yticks = 4:0.5:10, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :AfterLast_median, yerror = :AfterLast_CIq, group = :Treatment, color = :color, xflip = true, yticks = 4:0.5:10, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_Prot_AfterLast.pdf"))

@df dp scatter(:Protocol, :Leaving_NextPrew_mean, yerror = :Leaving_NextPrew_sem, group = :Treatment, color = :color, xflip = true, yticks = 0.1:0.1:0.5, ylims = (0.1,0.5), xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Leaving_NextPrew_median, yerror = :Leaving_NextPrew_CIq, group = :Treatment, color = :color, xflip = true, yticks = 0.1:0.1:0.5, ylims = (0.1,0.5), xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_Prot_PrewLeaving.pdf"))
################################ Fig4 global protocol effect ##################################
df1 = combine(groupby(streaks,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_Rewards)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end

df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(filter!(r->r.Phase in ["Optogenetic","Citalopram"],df2),showmice = false)
plot_wilcoxon(filter!(r->r.Phase in ["Optogenetic","Citalopram"],df2),showmice = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Wilcoxon_ROILeaving.pdf"))

streaks.Poking_Travel_to
streaks.Trial_Travel_to
streaks.ROI_Leaving_Time

df2 = combine(groupby(df1,[:Protocol,:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(filter!(r->r.Phase in ["Optogenetic","Citalopram"],df2),showmice = false)
################################ Fig4 scatter  global ##################################
gd1 = groupby(streaks,[:Protocol,:MouseID,:Phase,:Treatment])
df1 = combine(gd1, :Num_Rewards => mean => :Num_Rewards,
    :Num_pokes => mean => :Num_pokes,
    :AfterLast => mean => :AfterLast,
    :Leaving_NextPrew => mean => :Leaving_NextPrew,
    :AverageRewRate => mean => :AverageRewRate)

gd2 = groupby(df1,[:Protocol,:Treatment])

df2 = combine(gd2, [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> mean,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> sem,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> median,
    [:Num_Rewards,:Num_pokes,:AfterLast,:Leaving_NextPrew,:AverageRewRate] .=> CIq
    )

dp = filter!(r->r.Treatment in ["Optogenetic","Citalopram","Control"],df2)
Drug_colors!(dp)
@df dp scatter(:Protocol, :Num_Rewards_mean, yerror = :Num_Rewards_sem, group = :Treatment, color = :color, xflip = true, yticks = 1:0.5:6, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Num_Rewards_median, yerror = :Num_Rewards_CIq, group = :Treatment, color = :color, xflip = true, yticks = 1:0.5:6, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Prot_NumRewards.pdf"))

@df dp scatter(:Protocol, :Num_pokes_mean, yerror = :Num_pokes_sem, group = :Treatment, color = :color, xflip = true, yticks = 8:1:17, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Num_pokes_median, yerror = :Num_pokes_CIq, group = :Treatment, color = :color, xflip = true, yticks = 8:1:17, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Prot_NumPokes.pdf"))

@df dp scatter(:Protocol, :AfterLast_mean, yerror = :AfterLast_sem, group = :Treatment, color = :color, xflip = true, yticks = 4:0.5:10, xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :AfterLast_median, yerror = :AfterLast_CIq, group = :Treatment, color = :color, xflip = true, yticks = 4:0.5:10, xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Prot_AfterLast.pdf"))

@df dp scatter(:Protocol, :Leaving_NextPrew_mean, yerror = :Leaving_NextPrew_sem, group = :Treatment, color = :color, xflip = true, yticks = 0.1:0.1:0.6, ylims = (0.1,0.6), xlims = (0.25,2.75), label = false)
@df dp scatter(:Protocol, :Leaving_NextPrew_median, yerror = :Leaving_NextPrew_CIq, group = :Treatment, color = :color, xflip = true, yticks = 0.1:0.1:0.5, ylims = (0.1,0.5), xlims = (0.25,2.75), label = false)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Prot_PrewLeaving.pdf"))
################################### Selective Stats #####################################
test_SEL = filter(r->r.Phase in ["SB242084", "Altanserin"],streaks)
test_SEL[!,:NumericalProt] = parse.(Float64, test_Sel.Protocol)

NR_SEL_m1 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + (1|MouseID)),test_SEL))
NR_SEL_m2 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + NumericalProt + (1|MouseID)),test_SEL))
NR_SEL_m3 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + Treatment + NumericalProt + (1|MouseID)),test_SEL))
Likelyhood_Ratio_test(NR_SEL_m1,NR_SEL_m2)
Likelyhood_Ratio_test(NR_SEL_m2,NR_SEL_m3)

NP_SEL_m1 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + (1|MouseID)),test_SEL))
NP_SEL_m2 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt + (1|MouseID)),test_SEL))
NP_SEL_m3 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_SEL))
Likelyhood_Ratio_test(NP_SEL_m1,NP_SEL_m2)
Likelyhood_Ratio_test(NP_SEL_m2,NP_SEL_m3)

AL_SEL_m1 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + (1|MouseID)),test_SEL))
AL_SEL_m2 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt + (1|MouseID)),test_SEL))
AL_SEL_m3 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_SEL))
AL_SEL_m4 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt * Treatment + (1|MouseID)),test_SEL))
Likelyhood_Ratio_test(AL_SEL_m1,AL_SEL_m2)
Likelyhood_Ratio_test(AL_SEL_m2,AL_SEL_m3)
Likelyhood_Ratio_test(AL_SEL_m3,AL_SEL_m4)

PR_SEL_m1 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + (1|MouseID)),test_SEL))
PR_SEL_m2 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt + (1|MouseID)),test_SEL))
PR_SEL_m3 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_SEL))
PR_SEL_m4 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt * Treatment + (1|MouseID)),test_SEL))
Likelyhood_Ratio_test(PR_SEL_m1,PR_SEL_m2)
Likelyhood_Ratio_test(PR_SEL_m2,PR_SEL_m3)
Likelyhood_Ratio_test(PR_SEL_m3,PR_SEL_m4)
############################### Global stats #################################
test_GLO = filter(r->r.Phase in ["Optogenetic", "Citalopram"],streaks)
test_GLO[!,:NumericalProt] = parse.(Float64, test_GLO.Protocol)

NR_GLO_m1 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + (1|MouseID)),test_GLO))
NR_GLO_m2 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + NumericalProt + (1|MouseID)),test_GLO))
NR_GLO_m3 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + Treatment + NumericalProt + (1|MouseID)),test_GLO))
NR_GLO_m4 = fit!(LinearMixedModel(@formula(Num_Rewards ~ 1 + Treatment * NumericalProt + (1|MouseID)),test_GLO))
Likelyhood_Ratio_test(NR_GLO_m1,NR_GLO_m2)
Likelyhood_Ratio_test(NR_GLO_m2,NR_GLO_m3)
Likelyhood_Ratio_test(NR_GLO_m3,NR_GLO_m4)

NP_GLO_m1 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + (1|MouseID)),test_GLO))
NP_GLO_m2 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt + (1|MouseID)),test_GLO))
NP_GLO_m3 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_GLO))
NP_GLO_m4 = fit!(LinearMixedModel(@formula(Num_pokes ~ 1 + NumericalProt * Treatment + (1|MouseID)),test_GLO))
Likelyhood_Ratio_test(NP_GLO_m1,NP_GLO_m2)
Likelyhood_Ratio_test(NP_GLO_m2,NP_GLO_m3)
Likelyhood_Ratio_test(NP_GLO_m3,NP_GLO_m4)

AL_GLO_m1 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + (1|MouseID)),test_GLO))
AL_GLO_m2 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt + (1|MouseID)),test_GLO))
AL_GLO_m3 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_GLO))
AL_GLO_m4 = fit!(LinearMixedModel(@formula(AfterLast ~ 1 + NumericalProt * Treatment + (1|MouseID)),test_GLO))
Likelyhood_Ratio_test(AL_GLO_m1,AL_GLO_m2)
Likelyhood_Ratio_test(AL_GLO_m2,AL_GLO_m3)
Likelyhood_Ratio_test(AL_GLO_m3,AL_GLO_m4)

PR_GLO_m1 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + (1|MouseID)),test_GLO))
PR_GLO_m2 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt + (1|MouseID)),test_GLO))
PR_GLO_m3 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt + Treatment + (1|MouseID)),test_GLO))
PR_GLO_m4 = fit!(LinearMixedModel(@formula(Leaving_NextPrew ~ 1 + NumericalProt * Treatment + (1|MouseID)),test_GLO))
Likelyhood_Ratio_test(PR_GLO_m1,PR_GLO_m2)
Likelyhood_Ratio_test(PR_GLO_m2,PR_GLO_m3)
Likelyhood_Ratio_test(PR_GLO_m3,PR_GLO_m4)

################################ Fig3 test selective ##################################
s = filter(r -> r.Phase in ["Altanserin", "SB242084","Way_100135"], streaks)
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_test_Numpokes.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_test_AfterLast.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_test_PokingTravel.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig3","Sel_test_ROITravel.pdf"))

################################ Fig4 test Global ##################################
s = filter(r -> r.Phase in ["Optogenetic","Citalopram"], streaks)
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Test_Numpokes.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Test_AfterLast.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Test_PokingTravel.pdf"))

df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2; sorting = Plotting_position)
savefig(joinpath(figs_loc,"LabMeetingJan2021","Fig4","Glob_Test_ROITravel.pdf"))


text:
"""
$(TYPEDEF)

Defines a total Liouvillian to feed to the solver using the `DiffEqOperator` interface. It contains both closed-system and open-system Liouvillians.

# Fields

$(FIELDS)
"""
struct DiffEqLiouvillian{diagonalization,adiabatic_frame}
    "Hamiltonian"
    H::AbstractHamiltonian
    "Open system in eigenbasis"
    opensys_eig::Vector{AbstractLiouvillian}
    "Open system in normal basis"
    opensys::Vector{AbstractLiouvillian}
    "Levels to truncate"
    lvl::Int
    "Internal cache"
    u_cache::Any
end

"""
$(SIGNATURES)

The constructor of the `DiffEqLiouvillian` type. `opensys_eig` is a list of open-system Liouvillians that which require diagonalization of the Hamiltonian. `opensys` is a list of open-system Liouvillians which does not require diagonalization of the Hamiltonian. `lvl` is the truncation levels of the energy eigenbasis if the method supports the truncation.
"""
function DiffEqLiouvillian(H, opensys_eig, opensys, lvl)
    # for DenseHamiltonian smaller than 10×10, do not truncate
    if !(typeof(H) <: AbstractSparseHamiltonian) && (size(H, 1) <= 10)
        lvl = size(H, 1)
        u_cache = similar(get_cache(H))
    else
        lvl = size(H, 1) < lvl ? size(H, 1) : lvl
        # for SparseHamiltonian, we will create dense matrix cache
        # for the truncated subspace
        u_cache = Matrix{eltype(H)}(undef, lvl, lvl)
    end
    diagonalization = isempty(opensys_eig) ? false : true
    adiabatic_frame = typeof(H) <: AdiabaticFrameHamiltonian
    DiffEqLiouvillian{diagonalization,adiabatic_frame}(H, opensys_eig, opensys, lvl, u_cache)
end

function (Op::DiffEqLiouvillian{false,false})(du, u, p, t)
    s = p(t)
    Op.H(du, u, p, s)
    for lv in Op.opensys
        lv(du, u, p, t)
    end
end

function update_cache!(cache, Op::DiffEqLiouvillian{false,false}, p, t)
    update_cache!(cache, Op.H, p, p(t))
    for lv in Op.opensys
        update_cache!(cache, lv, p, t)
    end
end

function update_vectorized_cache!(cache, Op::DiffEqLiouvillian{false,false}, p, t)
    update_vectorized_cache!(cache, Op.H, p, p(t))
    for lv in Op.opensys
        update_vectorized_cache!(cache, lv, p, t)
    end
end

function (Op::DiffEqLiouvillian{true,false})(du, u, p, t)
    s = p(t)
    w, v = Op.H.EIGS(Op.H, s, Op.lvl)
    ω_ba = transpose(w) .- w
    ρ = v' * u * v
    H = Diagonal(w)
    Op.u_cache .= -1.0im * (H * ρ - ρ * H)
    for lv in Op.opensys_eig
        lv(Op.u_cache, ρ, ω_ba, v, s)
    end
    mul!(du, v, Op.u_cache * v')
    for lv in Op.opensys
        lv(du, u, p, t)
    end
end

function update_cache!(cache, Op::DiffEqLiouvillian{true,false}, p, t::Real)
    s = p(t)
    w, v = Op.H.EIGS(Op.H, s, Op.lvl)
    ω_ba = transpose(w) .- w
    # initialze the cache as Hamiltonian in eigenbasis
    fill!(Op.u_cache, 0.0)
    for i = 1:length(w)
        @inbounds Op.u_cache[i, i] = -1.0im * w[i]
    end
    for lv in Op.opensys_eig
        update_cache!(Op.u_cache, lv, ω_ba, v, s)
    end
    mul!(cache, v, Op.u_cache * v')
    for lv in Op.opensys
        update_cache!(cache, lv, p, t)
    end
end

function (Op::DiffEqLiouvillian{true,true})(du, u, p, t)
    s = p(t)
    H = Op.H(p.tf, s)
    w = diag(H)
    ω_ba = transpose(w) .- w
    du .= -1.0im * (H * u - u * H)
    for lv in Op.opensys_eig
        lv(du, u, ω_ba, s)
    end
end

# """
# $(TYPEDEF)

# Defines an open system DiffEqOperator in the multiple eigenbasis.

# # Fields

# $(FIELDS)
# """
# struct OpenSysOpHybrid{adiabatic_frame}
#     """Hamiltonian"""
#     H::AbstractHamiltonian
#     """Open system in energy eigenbasis"""
#     opensys_eig::Vector{AbstractLiouvillian}
#     """Open system in normal basis"""
#     opensys::Vector{AbstractLiouvillian}
#     """Levels to truncate"""
#     lvl::Int
#     """Internal cache"""
#     u_cache::Any
# end

# function OpenSysOpHybrid(H, open_eig, op, lvl)
#     # for DenseHamiltonian smaller than 10×10, do not truncate
#     if !(typeof(H) <: AbstractSparseHamiltonian) && (size(H, 1) <= 10)
#         lvl = size(H, 1)
#         u_cache = similar(get_cache(H))
#     else
#         # for SparseHamiltonian, we will create dense matrix cache
#         # for the truncated subspace
#         u_cache = Matrix{eltype(H)}(undef, lvl, lvl)
#     end
#     adiabatic_frame = typeof(H) <: AdiabaticFrameHamiltonian
#     OpenSysOpHybrid{adiabatic_frame}(H, open_eig, op, lvl, u_cache)
# end

# function update_cache!(cache, Op::OpenSysOpHybrid{false}, p, t::Real)
#     s = p(t)
#     w, v = Op.H.EIGS(Op.H, s, Op.lvl)
#     ω_ba = transpose(w) .- w
#     # initialze the cache as Hamiltonian in eigenbasis
#     fill!(Op.u_cache, 0.0)
#     for i = 1:length(w)
#         @inbounds Op.u_cache[i, i] = -1.0im * w[i]
#     end
#     for lv in Op.opensys_eig
#         update_cache!(Op.u_cache, lv, ω_ba, v, s)
#     end
#     mul!(cache, v, Op.u_cache * v')
#     for lv in Op.opensys
#         update_cache!(cache, lv, p, t)
#     end
# end

ext: jl


text:
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))), "test", "files", "simple_plant.mtg")

@testset "delete_node!" begin
    # Delete a node:
    mtg = read_mtg(file)
    length_start = length(mtg)
    delete_node!(get_node(mtg, 3))
    @test length(mtg) == length_start - 1

    # Delete a node with two children:
    @test_logs (:warn, "Scale of the child node branched but its deleted parent was decomposing. Keep branching, please check if the decomposition is still correct.") delete_node!(get_node(mtg, 4))
    # There is a warning because we don't know how to properly link the leaf (node 5) to its new parent.

    # Furthermore, the link of the second child (node 6) was modified from following to decomposing:
    @test get_node(mtg, 6).MTG.link == "/"
end

@testset "delete_node! with a user function" begin
    mtg = read_mtg(file)
    delete_node!(get_node(mtg, 4), child_link_fun = node -> node.MTG.link)
    # The link of the children didn't change:
    @test get_node(mtg, 6).MTG.link == "<"

    delete_node!(get_node(mtg, 3), child_link_fun = node -> "+")
    # Both children links changes to branching:
    @test get_node(mtg, 5).MTG.link == "+"
    @test get_node(mtg, 6).MTG.link == "+"
end

@testset "delete_nodes!" begin
    mtg = read_mtg(file)
    length_start = length(mtg)
    delete_nodes!(mtg, scale = 2) # Will remove all nodes of scale 2
    @test get_node(mtg, 3) === nothing
    @test length(mtg) === length_start - 1

    # Delete the leaves:
    mtg = read_mtg(file)
    delete_nodes!(mtg, symbol = "Leaf")
    @test length(mtg) === length_start - 2

    # Delete with a function, here we delete all nodes that have a parent with Length < 0.2:
    mtg = read_mtg(file)
    delete_nodes!(
        mtg,
        filter_fun = function (node)
            if !isroot(node)
                node.parent[:Length] !== nothing ? node.parent[:Length] < 0.2 : false
            else
                false
            end
        end
    )
    @test length(mtg) === length_start - 3
end

@testset "prune! a node" begin
    file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))), "test", "files", "simple_plant.mtg")

    # Prune a node in the middle:
    mtg = read_mtg(file)
    prune!(get_node(mtg, 6))
    @test length(mtg) == 5

    # Prune a leaf node:
    mtg = read_mtg(file)
    prune!(get_node(mtg, 5))
    @test length(mtg) == 6

    # Prune the root node (we expect an error):
    @test_throws ErrorException prune!(read_mtg(file))
end


text:
using SpecialFunctions

abstract type Integrator end
abstract type GaussianQuadrature end
abstract type GaussianOrthogonalPolynomial end
abstract type IntegrationFunction end

struct GaussLaguerrePolynomial <: GaussianOrthogonalPolynomial
  s::Float64
end

get_alpha(poly::GaussLaguerrePolynomial, i::Int) = 2 * i + 1 + poly.s
get_beta(poly::GaussLaguerrePolynomial, i::Int) = i * (i + poly.s)
get_mu_0(poly::GaussLaguerrePolynomial) = exp(lgamma(poly.s + 1.0))
get_w(poly::GaussLaguerrePolynomial, x::Float64) = ^(x, poly.s) * exp(-x)

mutable struct SegmentIntegral <: Integrator
  absoluteAccuracy::Float64
  absoluteError::Float64
  maxEvals::Int
  evals::Int
  intervals::Int
end

SegmentIntegral(intervals::Int) = SegmentIntegral(1.0, 0.0, 1, 0, intervals)

function operator(integrator::Integrator, f::Function, a::Float64, b::Float64)
  integrator.evals = 0
  if a == b
    return 0.0
  end

  if (b > a)
    return integrate(integrator, f, a, b)
  else
    return -integrate(integrator, f, b, a)
  end
end

function (integrator::SegmentIntegral)(f::Function, a::Float64, b::Float64)
  integrator.evals = 0
  if a == b
    return 0.0
  end

  if (b > a)
    return integrate(integrator, f, a, b)
  else
    return -integrate(integrator, f, b, a)
  end
end

function Math.integrate(integrator::SegmentIntegral, f::Union{Function, IntegrationFunction}, a::Float64, b::Float64)
  dx = (b - a) / (integrator.intervals * 1.0)
  sum_ = 0.5 * (f(a) + f(b))
  end_ = b - 0.5 * dx
  x = a + dx
  while x < end_
    sum_ += f(x)
    x += dx
  end

  return sum_ * dx
end

struct GaussLaguerreIntegration <: GaussianQuadrature
  x::Vector{Float64}
  w::Vector{Float64}
end

function (gli::GaussLaguerreIntegration)(f::IntegrationFunction)
  sum_ = 0.0
  @inbounds @simd for i = length(gli.x):-1:1
    sum_ += gli.w[i] * f(gli.x[i])
  end

  return sum_
end

function GaussLaguerreIntegration(n::Int, s::Float64 = 0.0)
  x, w = build_gaussian_quadrature(n, GaussLaguerrePolynomial(s))

  return GaussLaguerreIntegration(x, w)
end

get_order(integration::GaussianQuadrature) = length(integration.x)

function build_gaussian_quadrature(n::Int, poly::GaussianOrthogonalPolynomial)
  x = zeros(n)
  w = zeros(n)

  e = zeros(n - 1)

  @simd for i = 2:n
    @inbounds x[i] = get_alpha(poly, i - 1)
    @inbounds e[i - 1] = sqrt(get_beta(poly, i - 1))
  end

  x[1] = get_alpha(poly, 0)

  # sym tridiagonal
  # tridiag = SymTridiagonal(x, e)
  # eigtri = eigfact(tridiag)

  tqr = TqrEigenDecomposition(x, e, OnlyFirstRowEigenVector(), Overrelaxation())
  x_new = tqr.d
  ev = tqr.ev

  mu_0 = get_mu_0(poly)

  @simd for i = 1:n
    @inbounds w[i] = mu_0 * ev[i, 1] * ev[i, 1] / get_w(poly, x_new[i])
  end

  return x_new, w
end


text:
module Bodies

import Base: parent

using AstroTime: SECONDS_PER_DAY, SECONDS_PER_CENTURY, value, j2000, seconds
using ItemGraphs: ItemGraph, SimpleGraph, add_edge!, add_vertex!, items

export CelestialBody,
    NAIFId,
    SolarSystemBarycenter,
    Sun,
    declination,
    declination_rate,
    euler_angles,
    euler_rates,
    from_naifid,
    grav_param,
    naifid,
    parent,
    right_ascension,
    right_ascension_rate,
    rotation_angle,
    rotation_rate,
    ssb,
    sun

const NAIFId = Int
const BODIES = ItemGraph{NAIFId, NAIFId}(SimpleGraph())

abstract type CelestialBody end

function grav_param end
function mean_radius end
function polar_radius end
function equatorial_radius end

Base.show(io::IO, body::CelestialBody) = print(io, string(nameof(typeof(body))))

from_naifid(id::NAIFId) = from_naifid(Val(id))

path_ids(from::CelestialBody, to::CelestialBody) = items(BODIES, naifid(from), naifid(to))

abstract type Barycenter <: CelestialBody end

struct SolarSystemBarycenter <: Barycenter end
const ssb = SolarSystemBarycenter()
Base.show(io::IO, ::SolarSystemBarycenter) = print(io, "Solar System Barycenter")
naifid(::SolarSystemBarycenter) = 0
from_naifid(::Val{0}) = ssb

struct Sun <: CelestialBody end
const sun = Sun()
parent(::Sun) = ssb
naifid(::Sun) = 10
from_naifid(::Val{10}) = sun

alpha0(::CelestialBody) = 0.0
alpha1(::CelestialBody) = 0.0
alpha2(::CelestialBody) = 0.0
delta0(::CelestialBody) = 0.0
delta1(::CelestialBody) = 0.0
delta2(::CelestialBody) = 0.0
omega0(::CelestialBody) = 0.0
omega1(::CelestialBody) = 0.0
omega2(::CelestialBody) = 0.0
alpha(::CelestialBody) = zeros(1)
delta(::CelestialBody) = zeros(1)
omega(::CelestialBody) = zeros(1)
theta0(::CelestialBody) = zeros(1)
theta1(::CelestialBody) = zeros(1)

const Δtc = SECONDS_PER_CENTURY
const Δtd = SECONDS_PER_DAY

theta(b::CelestialBody, t) = theta0(b) .+ theta1(b) .* t / Δtc

function right_ascension(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    α = alpha(b)
    α₀ = alpha0(b)
    α₁ = alpha1(b)
    α₂ = alpha2(b)
    θ = theta(b, t)
    α₀ + α₁ * t / Δtc + α₂ * t^2 / Δtc^2 + sum(α .* sin.(θ))
end

function declination(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    δ = delta(b)
    δ₀ = delta0(b)
    δ₁ = delta1(b)
    δ₂ = delta2(b)
    θ = theta(b, t)
    δ₀ + δ₁ * t / Δtc + δ₂ * t^2 / Δtc^2 + sum(δ .* cos.(θ))
end

function rotation_angle(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    ω = omega(b)
    ω₀ = omega0(b)
    ω₁ = omega1(b)
    ω₂ = omega2(b)
    θ = theta(b, t)
    ω₀ + ω₁ * t / Δtd + ω₂ * t^2 / Δtd^2 + sum(ω .* sin.(θ))
end

function right_ascension_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    α = alpha(b)
    α₁ = alpha1(b)
    α₂ = alpha2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    α₁ / Δtc + 2 * α₂ * t / Δtc^2 + sum(α .* θ₁ ./ Δtc .* cos.(θ))
end

function declination_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    δ = delta(b)
    δ₁ = delta1(b)
    δ₂ = delta2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    δ₁ / Δtc + 2 * δ₂ * t / Δtc^2 - sum(δ .* θ₁ ./ Δtc .* sin.(θ))
end

function rotation_rate(b::CelestialBody, ep)
    t = value(seconds(j2000(ep)))
    ω = omega(b)
    ω₁ = omega1(b)
    ω₂ = omega2(b)
    θ = theta(b, t)
    θ₁ = theta1(b)
    ω₁ / Δtd + 2 * ω₂ * t / Δtd^2 + sum(ω .* θ₁ ./ Δtc .* cos.(θ))
end

function euler_angles(b::CelestialBody, ep)
    right_ascension(b, ep) + π/2, π/2 - declination(b, ep), mod2pi(rotation_angle(b, ep))
end

function euler_rates(b::CelestialBody, ep)
    right_ascension_rate(b, ep), -declination_rate(b, ep), rotation_rate(b, ep)
end

include("planets.jl")
include("minor.jl")
include("satellites.jl")

const ALL_NAMES = Tuple([["Sun", "SolarSystemBarycenter"]
                         collect(PLANET_NAMES .* "Barycenter");
                         collect(PLANET_NAMES);
                         ["Luna", "Phobos", "Deimos"];
                         collect(JUPITER_SATELLITE_NAMES);
                         collect(SATURN_SATELLITE_NAMES);
                         collect(URANUS_SATELLITE_NAMES);
                         collect(NEPTUNE_SATELLITE_NAMES);
                         collect(MINOR_BODY_NAMES);
                         ["Pluto", "PlutoBarycenter"];
                         collect(PLUTO_SATELLITE_NAMES)])

gm = joinpath(@__DIR__, "..", "..", "gen", "gm.jl")
isfile(gm) || error("`gm.jl` has not been generated, yet. Run `gen/gen_gm.jl`")
include(gm)

pck = joinpath(@__DIR__, "..", "..", "gen", "pck.jl")
isfile(pck) || error("`pck.jl` has not been generated, yet. Run `gen/gen_pck.jl`")
include(pck)

function __init__()
    # Sun and SSB
    add_vertex!(BODIES, 0)
    add_edge!(BODIES, 0, 10)

    # Planets and barycenters
    for i in 1:length(PLANET_NAMES)
        id = 100i + 99
        add_edge!(BODIES, 0, i)
        add_edge!(BODIES, i, id)
    end

    # Earth satellite
    add_edge!(BODIES, 3, 301)

    # Mars satellites
    add_edge!(BODIES, 4, 401)
    add_edge!(BODIES, 4, 402)

    # Jupiter satellites
    for i in 1:length(JUPITER_SATELLITE_NAMES[1:end-1])
        id = 500 + i
        add_edge!(BODIES, 5, id)
    end
    add_edge!(BODIES, 5, 553)

    # Saturn satellites
    for i in 1:length(SATURN_SATELLITE_NAMES)
        id = 600 + i
        add_edge!(BODIES, 6, id)
    end

    # Uranus satellites
    for i in 1:length(URANUS_SATELLITE_NAMES)
        id = 700 + i
        add_edge!(BODIES, 7, id)
    end

    # Neptune satellites
    for i in 1:length(NEPTUNE_SATELLITE_NAMES)
        id = 800 + i
        add_edge!(BODIES, 8, id)
    end

    # Pluto satellites
    for i in 1:length(PLUTO_SATELLITE_NAMES)
        id = 900 + i
        add_edge!(BODIES, 9, id)
    end

    # Minor bodies
    add_edge!(BODIES, 0, 9)
    add_edge!(BODIES, 9, 999)
    add_edge!(BODIES, 0, 2000001)
    add_edge!(BODIES, 0, 2000002)
    add_edge!(BODIES, 0, 2000004)
    add_edge!(BODIES, 0, 2000016)
    add_edge!(BODIES, 0, 2000021)
    add_edge!(BODIES, 0, 2431010)
    add_edge!(BODIES, 0, 2000433)
    add_edge!(BODIES, 0, 2000511)
    add_edge!(BODIES, 0, 9511010)
    add_edge!(BODIES, 0, 2002867)
    add_edge!(BODIES, 0, 2025143)
    add_edge!(BODIES, 0, 1000093)
    add_edge!(BODIES, 0, 1000005)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ dce56652-31d4-11eb-365b-6100a847b13d
using LightGraphs

# ╔═╡ 9c4b53ae-31d4-11eb-052a-abd726285a2d
using GraphRecipes

# ╔═╡ a9737288-31d5-11eb-1ca8-596a083b9dbb
using PGFPlotsX

# ╔═╡ c510f6ae-31d4-11eb-1901-5563116763d3
using Plots

# ╔═╡ 237d50d0-31c5-11eb-0055-6f6e97b5f479
md" ## _Lesson 1: Variable and Types_
In this lesson, you will learn about variables and types and how to utilize them in Julia. You will also learn about basic mathematical operations
"

# ╔═╡ 836d87b2-31c5-11eb-203e-0b3425e33f4c
md" ## Variables
A program is a series of instructions to be performed on various kind of data. The data can range from names, to heights, to actual scientific data. To better organize the data, you can use variables, which are labels that the data is assigned to. To better understand this, lets look at some examples.
"

# ╔═╡ 4404e28e-31c6-11eb-2dfe-fd3c843495f3
md" **Example 1**
Lets say that you have a collection of data about all your favorites like favorite color, integer, and favorite scientific constant. To collect in a organized fasion you can assign the data to specific variables with labels, or names that tell people what kind of data the variable contains."

# ╔═╡ 9be8b67a-31c7-11eb-3117-d900b81590f4
favorite_color = "blue"

# ╔═╡ 8aa726e8-31c8-11eb-2866-a12e7f701e49
md" This variable contains a string because it uses quotation marks."

# ╔═╡ c6a234a4-31c7-11eb-114e-77493d1e92b4
favorite_integer = 11

# ╔═╡ 8d8ae334-31c8-11eb-0c90-2339c7117f4d
md" This variable contains an integer."

# ╔═╡ ce4f411a-31c7-11eb-2748-3d045a819acb
favorite_science_constant = 0.0821

# ╔═╡ 935e4d16-31c8-11eb-3b0b-495838d64c69
md" This variable contains a floating point number. It is a floating point number because it is a decimal."

# ╔═╡ 587630e2-31d2-11eb-0e2a-11d70ebdb31b
md" Then, to actually print that value, just type `print(variable name)`. In this notebook however, you do not need the print, just type the variable name"

# ╔═╡ 62f3be04-31d2-11eb-37f0-89d5d01c3e7c
favorite_color, favorite_integer, favorite_science_constant

# ╔═╡ cc46e5f2-31d2-11eb-3b29-e1555c25860f
md" **Example 2**
Using variables, you can perform basic math problems and assign the answers to those problems to new variables. Look below to see how variables can be manipulated to perform arithmatic. First we assign the integers 2 and 3 to the variables a and b.
"

# ╔═╡ f5213ffe-31d2-11eb-3b28-abe9bfdb25bf
a=2

# ╔═╡ 164132ca-31d3-11eb-1865-cb1f15921ebe
b=3

# ╔═╡ 2d261cd0-31d3-11eb-287b-2d9d5546f682
md" Once that is finished, you can now add, subtract, multiply and divide using those variables and you will get values.
"

# ╔═╡ 437f37f0-31d3-11eb-3eac-0fbffcd63384
sum = a + b

# ╔═╡ 51c0d60c-31d3-11eb-13de-4935e1795602
difference = a - b

# ╔═╡ 584dc1a6-31d3-11eb-0c72-a5797f0ff9a3
product = a * b

# ╔═╡ 604a1e18-31d3-11eb-17b1-f1da04ca048e
quotient = a / b

# ╔═╡ 74d7ad46-31d3-11eb-2a26-9788acd07194
power = a^2

# ╔═╡ 7a44d362-31d3-11eb-2006-2f79e5383f95
modulus = a%b

# ╔═╡ 84cf7bca-31d3-11eb-3594-a722da40b1bf
md" Additionally, you can use `a//b` if you are looking for a value exact, instead of in decimal form. "

# ╔═╡ 90c953fa-31c5-11eb-0916-b7b3828e12e1
md" ## Types
Evey element, or variable must have a type. The type depends on the hierarchial system.

"

# ╔═╡ 98c12592-31c5-11eb-0be0-bdf59ba1274a
begin
g = DiGraph(13)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
add_edge!(g, 3, 4)
add_edge!(g, 3, 7)
# Instead of adding an edge 5 -> 3, have an edge that goes from 5 to
# the ghost node, then to 3.
add_edge!(g, 4, 5)
add_edge!(g, 4, 6)
add_edge!(g, 7, 8)
add_edge!(g, 7, 9)
add_edge!(g, 8, 10)
add_edge!(g, 8, 11)
add_edge!(g, 9, 12)
add_edge!(g, 9, 13)

# Add an edge by itself that will force GraphRecipes to have a decent
# window size.
add_edge!(g, 7, 7)
add_edge!(g, 8, 8)
add_edge!(g, 9, 9)
add_edge!(g, 10, 10)
add_edge!(g, 11, 11)
add_edge!(g, 12, 12)
add_edge!(g, 13, 13)

names = [:Any, :Number, :Real, :abstractfloat, :float64, :float32, :Integer, :signed, :unsigned, :Int64, :Int32, :UInt64, :UInt32]
x = [0,0,0,0,-1,1,5,4,6,3.25,4.75,6,8]
y = [5,4,3,2,1,1,3,2,2,1,1,1,1]
graphplot(g, x=x, y=y, curvature_scalar=0.0, nodesize=1.0,
          node_weights=[1.0,1.0,1.0,1.0,1.0,1.0,1.0 ,1.0,1.0,1.0,1.0,1.0,1.0],
          names=names, nodecolor=:lightyellow, color=:black,
          nodeshape=:rect, self_edge_size=0.0)
end

# ╔═╡ 9ff95ace-31e2-11eb-3e4a-d1fa507ff913
md" If a variable is declared as type number, anything below number on the hiearchy chart can be data in the type number since number is an overarching branch. To determine the type of variable, just use `typeof(number)`."

# ╔═╡ f76a077c-31e2-11eb-2bd0-399350c2e820
typeof(0.0921)

# ╔═╡ ff24c0a6-31e2-11eb-25e0-c5f530e27ed0
typeof(123)

# ╔═╡ 0a316b98-31e3-11eb-1efd-1554bc7b78b2
md"You do need to be careful about initializing the variable to the correct type. Julia takes care of most of it but there are several instances where the results you want will not happen. For example, if you assign a = 2, but later reassign it to 1.42, the decimal part will be cut off because since you originally assigned a to an int, that means that a is an int, so storing a floating point number will be truncated. It is easier just to plan ahead but incase you forget, you can utilize `convert` to change the type to match the number."

# ╔═╡ cb4b9cd6-31e3-11eb-1ce1-453f4af82cee
convert(Float64,2)

# ╔═╡ dc16eef8-31e3-11eb-1eba-737e85238d07
c=4

# ╔═╡ e05ed4e4-31e3-11eb-1f62-d987a3d416fe
d = convert(Float64,c)

# ╔═╡ 86cd5ffc-31c4-11eb-1394-c7ddfa5dc534
md" ## Acknowledgements
1. Learning Julia notebooks based on content from [Zero to Julia] (https://techytok.com/from-zero-to-julia/)
2. https://discourse.julialang.org/t/flow-charts-lightgraphs/31139/9"

# ╔═╡ a36c5818-31d4-11eb-0de0-35c9674d733f
import Pkg

# ╔═╡ c9c55d7a-31d4-11eb-0a65-abd881438331
Pkg.add("LightGraphs")

# ╔═╡ a7a054ca-31d4-11eb-317b-cf624d47791a
Pkg.add("GraphRecipes")

# ╔═╡ b4575984-31d4-11eb-0045-fb0824a08c20
Pkg.add("Plots")

# ╔═╡ 9fae0234-31d5-11eb-21ba-abf6e13932c1
Pkg.add("PGFPlotsX")

# ╔═╡ Cell order:
# ╟─237d50d0-31c5-11eb-0055-6f6e97b5f479
# ╟─836d87b2-31c5-11eb-203e-0b3425e33f4c
# ╟─4404e28e-31c6-11eb-2dfe-fd3c843495f3
# ╠═9be8b67a-31c7-11eb-3117-d900b81590f4
# ╟─8aa726e8-31c8-11eb-2866-a12e7f701e49
# ╠═c6a234a4-31c7-11eb-114e-77493d1e92b4
# ╟─8d8ae334-31c8-11eb-0c90-2339c7117f4d
# ╠═ce4f411a-31c7-11eb-2748-3d045a819acb
# ╟─935e4d16-31c8-11eb-3b0b-495838d64c69
# ╟─587630e2-31d2-11eb-0e2a-11d70ebdb31b
# ╠═62f3be04-31d2-11eb-37f0-89d5d01c3e7c
# ╟─cc46e5f2-31d2-11eb-3b29-e1555c25860f
# ╠═f5213ffe-31d2-11eb-3b28-abe9bfdb25bf
# ╠═164132ca-31d3-11eb-1865-cb1f15921ebe
# ╟─2d261cd0-31d3-11eb-287b-2d9d5546f682
# ╠═437f37f0-31d3-11eb-3eac-0fbffcd63384
# ╠═51c0d60c-31d3-11eb-13de-4935e1795602
# ╠═584dc1a6-31d3-11eb-0c72-a5797f0ff9a3
# ╠═604a1e18-31d3-11eb-17b1-f1da04ca048e
# ╠═74d7ad46-31d3-11eb-2a26-9788acd07194
# ╠═7a44d362-31d3-11eb-2006-2f79e5383f95
# ╟─84cf7bca-31d3-11eb-3594-a722da40b1bf
# ╟─90c953fa-31c5-11eb-0916-b7b3828e12e1
# ╠═98c12592-31c5-11eb-0be0-bdf59ba1274a
# ╟─9ff95ace-31e2-11eb-3e4a-d1fa507ff913
# ╠═f76a077c-31e2-11eb-2bd0-399350c2e820
# ╠═ff24c0a6-31e2-11eb-25e0-c5f530e27ed0
# ╟─0a316b98-31e3-11eb-1efd-1554bc7b78b2
# ╠═cb4b9cd6-31e3-11eb-1ce1-453f4af82cee
# ╠═dc16eef8-31e3-11eb-1eba-737e85238d07
# ╠═e05ed4e4-31e3-11eb-1f62-d987a3d416fe
# ╟─86cd5ffc-31c4-11eb-1394-c7ddfa5dc534
# ╠═c9c55d7a-31d4-11eb-0a65-abd881438331
# ╠═dce56652-31d4-11eb-365b-6100a847b13d
# ╠═a7a054ca-31d4-11eb-317b-cf624d47791a
# ╠═b4575984-31d4-11eb-0045-fb0824a08c20
# ╠═9c4b53ae-31d4-11eb-052a-abd726285a2d
# ╠═9fae0234-31d5-11eb-21ba-abf6e13932c1
# ╠═a9737288-31d5-11eb-1ca8-596a083b9dbb
# ╠═c510f6ae-31d4-11eb-1901-5563116763d3
# ╠═a36c5818-31d4-11eb-0de0-35c9674d733f


text:
module GpSelection
    export fixMatrix, calcSigma, generateZ, generateData

    using Distributions
    using MatrixDepot
    using Combinatorics

    include("utils.jl")
    include("sampledata.jl")
    include("varBayes.jl")
end


text:
function check_method(local_search_method::String)
    if !(local_search_method in ["mutation", "powell", "de", "cmaes"])
        error(
            "\"$local_search_method\": Invalid local_search_method. "
            * "Should be one of [\"mutaion\", \"Powell\", \"DE\", \"CMAES\"]"
        )
    elseif local_search_method in ["powell", "de"] && !isinstalled("scipy.optimize")
        error(
            "Cannnot import scipy.optimize functions. Use \"mutation\" for local search."
        )
    end
end

function optimize(
        model::Model,
        nth_param_set::Int64;
        popsize::Int64=5,
        max_generation::Int64=10000,
        initial_threshold::Float64=1e12,
        allowable_error::Float64=0.0,
        n_children::Int64=50,
        maxiter::Int64=100,
        local_search_method::String="mutation")
    check_method(
        lowercase(local_search_method)
    )

    for dir in ["fitparam", "logs"]
        if !isdir(
            joinpath(
                model.path,
                dir
            )
        )
            mkdir(
                joinpath(
                    model.path,
                    dir
                )
            )
        end
    end

    try
        files = readdir(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set"
            )
        )
        for file in files
            if occursin(".dat", file)
                rm(
                    joinpath(
                        model.path,
                        "fitparam",
                        "$nth_param_set",
                        "$file"
                    )
                )
            end
        end
    catch
        mkdir(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set"
            )
        )
    end

    search_rgn::Matrix{Float64} = model.search_region()

    n_population::Int64 = popsize * size(search_rgn, 2)
    n_gene::Int64 = size(search_rgn, 2)

    ga_v2(
        model=model,
        nth_param_set=nth_param_set,
        max_generation=max_generation,
        n_population=n_population,
        n_gene=n_gene,
        initial_threshold=initial_threshold,
        allowable_error=allowable_error,
        local_search_method=lowercase(local_search_method),
        n_children=n_children,
        maxiter=maxiter
    )
end

function optimize_continue(
        model::Model,
        nth_param_set::Int64;
        popsize::Int64=5,
        max_generation::Int64=10000,
        initial_threshold::Float64=1e12,
        allowable_error::Float64=0.0,
        n_children::Int64=50,
        maxiter::Int64=100,
        local_search_method::String="mutation",
        p0_bounds::Vector{Float64}=[0.1, 10.0])
    check_method(
        lowercase(local_search_method)
    )

    search_rgn::Matrix{Float64} = model.search_region()

    n_population::Int64 = popsize * size(search_rgn, 2)
    n_gene::Int64 = size(search_rgn, 2)

    if !isdir(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set"
        )
    )
        mkdir(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set"
            )
        )

        ga_v2(
            model=model,
            nth_param_set=nth_param_set,
            max_generation=max_generation,
            n_population=n_population,
            n_gene=n_gene,
            initial_threshold=initial_threshold,
            allowable_error=allowable_error,
            local_search_method=lowercase(local_search_method),
            n_children=n_children,
            maxiter=maxiter
        )
    else
        ga_v2_continue(
            model=model,
            nth_param_set=nth_param_set,
            max_generation=max_generation,
            n_population=n_population,
            n_gene=n_gene,
            initial_threshold=initial_threshold,
            allowable_error=allowable_error,
            local_search_method=lowercase(local_search_method),
            n_children=n_children,
            maxiter=maxiter,
            p0_bounds=p0_bounds
        )
    end
end


text:
struct Once{T,F}
    factory::F
    promise::Promise{T}
end

function Once{T}() where {T}
    T isa Type || _once_invalid_type_parameter(T)
    return Once{T,Type{T}}(T, Promise{T}())
end

function Once{T}(f) where {T}
    T isa Type || _once_invalid_type_parameter(T)
    return Once{T,_typeof(f)}(f, Promise{T}())
end

@noinline _once_invalid_type_parameter(@nospecialize T) =
    error("`Once{T}`` expcet a type for `T`; got T = $T")

function Once(f)
    value = f()
    T = typeof(value)
    promise = Promise{T}()
    # Not using the `value` calculated above in `Promise{T}` to avoid serializing into the
    # precompiled module just in case it is used as in `const O = Once(f)` in a package.
    return Once{T,_typeof(f)}(f, promise)
end

Base.getindex(once::Once) = race_put_with!(once.factory, once.promise)

function Base.delete!(once::Once)
    @atomic :monotonic once.promise.value = NOTSET
    return once
end


text:
Base.compilecache("CUDex")
try success(`nvcc --version`)
    cd("../src/cuda") do
        run(`make libcudex.so`)
    end
catch
    warn("CUDA not installed, GPU support will not be available.")
end


text:
function likelihood(
    examinees::Vector{<:AbstractExaminee}
    )
    return mapreduce( e -> _log_c(e.latent.likelihood), +, examinees)::Float64
end

function _likelihood(
    latents::Vector{<:AbstractLatent}
    )
    return mapreduce( l -> _log_c(l.likelihood), +, latents)::Float64
end


text:
struct MixedModelBootstrap{T<:AbstractFloat}
    m::LinearMixedModel{T}
    bstr::Tables.ColumnTable
    cnamedict::Dict{String,Int}
end

function MixedModelBootstrap(m::LinearMixedModel, bstr::NamedTuple) where {T}
    dict = Dict{String,Int}()
    k = findfirst(isequal(:β₁), keys(bstr))
    isnothing(k) && throw(ArgumentError("No :β₁ in keys(bstr)"))
    for nm in coefnames(m)
        dict[nm] = k
        k += 1
    end
    MixedModelBootstrap(m, bstr, dict)
end

function Base.getproperty(bsamp::MixedModelBootstrap, s::Symbol)
    if s == :model
        getfield(bsamp, :m)
    elseif s ∈ [:objective, :σ, :θ]
        getproperty(getfield(bsamp, :bstr), s)
    elseif s == :σs
        σs(bsamp)
    elseif s == :σρs
        σρs(bsamp)
    elseif haskey(getfield(bsamp, :cnamedict), string(s))
        getfield(bsamp, :bstr)[getfield(bsamp, :cnamedict)[string(s)]]
    else
        getfield(bsamp, s)
    end
end

"""
    parametricbootstrap(rng::AbstractRNG, nsamp::Integer, m::LinearMixedModel,
        props=(:objective, :σ, :β, :θ); β = m.β, σ = m.σ, θ = m.θ)
    parametricbootstrap(nsamp::Integer, m::LinearMixedModel,
        props=(:objective, :σ, :β, :θ); β = m.β, σ = m.σ, θ = m.θ)

Perform `nsamp` parametric bootstrap replication fits of `m`, returning a
`Tables.ColumnTable` of parameter estimates of the refit model.

The default random number generator is `Random.GLOBAL_RNG`.

# Named Arguments

`β`, `σ`, and `θ` are the values of `m`'s parameters for simulating the responses.
"""
function parametricbootstrap(rng::AbstractRNG, nsamp::Integer, m::LinearMixedModel{T};
    β = m.β, σ = m.σ, θ = m.θ) where {T}
    y₀ = copy(response(m))  # to restore original state of m
    θscr = copy(θ)
    βscr = copy(β)
    k = length(θ)
    bnms = Symbol.(subscriptednames("β", length(β)))
    vnms = (:objective, :σ, bnms..., :θ)
    value = (;(nm => nm == :θ ? SVector{k,T}[] : Vector{T}(undef, nsamp)
        for nm in vnms)...)
    try
        @showprogress 1 for i in 1:nsamp
            refit!(simulate!(rng, m, β = β, σ = σ, θ = θ))
            value.objective[i] = objective(m)
            value.σ[i] = sdest(m)
            fixef!(βscr, m)
            for (j, bnm) in enumerate(bnms)
                getproperty(value, bnm)[i] = βscr[j]
            end
            push!(value.θ, SVector{k, T}(getθ!(θscr, m)))
        end
    finally
        refit!(m, y₀)
    end
    MixedModelBootstrap(deepcopy(m), value)
end

function mktable(nsamp, p, k, T)
    nms = (:objective, :σ, Symbol.(subscriptednames("β", p))..., :θ)
    (; (nm => nm == :θ ? Vector{SVector{k,T}}(undef, nsamp) : Vector{T}(undef, nsamp) for nm in nms)...)
end

function parametricbootstrap(nsamp::Integer, m::LinearMixedModel, β = m.β, σ = m.σ, θ = m.θ)
    parametricbootstrap(Random.GLOBAL_RNG, nsamp, m, β = β, σ = σ, θ = θ)
end

function Base.propertynames(bsamp::MixedModelBootstrap)
    append!([:model, :objective, :σ, :θ, :σs, :σρs], Symbol.(keys(bsamp.cnamedict)))
end

function byreterm(bsamp::MixedModelBootstrap{T}, f::Function) where {T}
    m = bsamp.m
    oldθ = getθ(m)     # keep a copy to restore later
    retrms = m.reterms
    value = [typeof(v)[] for v in f.(retrms, m.σ)]
    for (σ,θ) in zip(bsamp.bstr.σ, bsamp.bstr.θ)
        setθ!(m, θ)
        for (i, v) in enumerate(f.(retrms, σ))
            push!(value[i], v)
        end
    end
    refit!(setθ!(m, oldθ))
    NamedTuple{(Symbol.(fnames(m))...,)}(value)
end

σs(bsamp::MixedModelBootstrap) = byreterm(bsamp, σs)

σρs(bsamp::MixedModelBootstrap) = byreterm(bsamp, σρs)

"""
    shortestCovInt(v, level = 0.95)

Return the shortest interval containing `level` proportion of the values of `v`
"""
function shortestCovInt(v, level = 0.95)
    n = length(v)
    0 < level < 1 || throw(ArgumentError("level = $level should be in (0,1)"))
    vv = issorted(v) ? v : sort(v)
    ilen = Int(ceil(n * level))   # the length of the interval in indices
    len, i = findmin([vv[i + ilen - 1] - vv[i] for i in 1:(n + 1 - ilen)])
    vv[[i, i + ilen - 1]]
end

"""
    simulate!(rng::AbstractRNG, m::LinearMixedModel{T}; β=m.β, σ=m.σ, θ=T[])
    simulate!(m::LinearMixedModel; β=m.β, σ=m.σ, θ=m.θ)

Overwrite the response (i.e. `m.trms[end]`) with a simulated response vector from model `m`.
"""
function simulate!(rng::AbstractRNG, m::LinearMixedModel{T}; β = coef(m), σ=m.σ, θ=T[]) where {T}
    isempty(θ) || setθ!(m, θ)
    y = randn!(rng, response(m))      # initialize y to standard normal
    for trm in m.reterms              # add the unscaled random effects
        unscaledre!(rng, y, trm)
    end
                    # scale by σ and add fixed-effects contribution
    BLAS.gemv!('N', one(T), m.X, β, σ, y)
    m
end

function simulate!(m::LinearMixedModel{T}; β=m.β, σ=m.σ, θ=T[]) where {T}
    simulate!(Random.GLOBAL_RNG, m, β=β, σ=σ, θ=θ)
end

"""
    unscaledre!(y::AbstractVector{T}, M::ReMat{T}, b) where {T}
    unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, M::ReMat{T}) where {T}

Add unscaled random effects defined by `M` and `b` to `y`.  When `rng` is present the `b`
vector is generated as `randn(rng, size(M, 2))`
"""
function unscaledre! end

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, b::AbstractVector{T}) where {T}
    m, n = size(A)
    length(y) == m && length(b) == n || throw(DimensionMismatch(""))
    z = A.z
    @inbounds for (i, r) in enumerate(A.refs)
        y[i] += b[r] * z[i]
    end
    y
end

unscaledre!(y::AbstractVector{T}, A::ReMat{T,1}, B::AbstractMatrix{T}) where {T} =
    unscaledre!(y, A, vec(B))

function unscaledre!(y::AbstractVector{T}, A::ReMat{T,S}, b::AbstractMatrix{T}) where {T,S}
    Z = A.z
    k, n = size(Z)
    l = nlevs(A)
    length(y) == n && size(b) == (k, l) || throw(DimensionMismatch(""))
    @inbounds for (i, ii) in enumerate(A.refs)
        for j in 1:k
            y[i] += Z[j, i] * b[j, ii]
        end
    end
    y
end

unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,S}) where {T,S} =
    unscaledre!(y, A, lmul!(A.λ, randn(rng, S, nlevs(A))))

unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,1}) where {T} =
    unscaledre!(y, A, lmul!(first(A.λ), randn(rng, 1, nlevs(A))))

unscaledre!(y::AbstractVector, A::ReMat) = unscaledre!(Random.GLOBAL_RNG, y, A)


text:
using BSON
using RCall
@rimport graphics as plt
@rimport grDevices as dev
@rimport cytof3
@rimport rcommon
@rimport base
include("sample/sample_lam.jl")

function plotpng(fname, s=10, w=480, h=480, ps=12; kw...)
  dev.png(fname, w=w*s, h=h*s, pointsize=ps*s, kw...)
end
addcut(clus, s=1) = plt.abline(h=base.cumsum(base.table(clus)) .+ .5, lwd=3*s, col="yellow")

# For BSON
using Cytof5, Flux, Distributions

function post_process(output_path; thresh=.99)
  out = BSON.load(output_path)
  if :simdat in keys(out)
    has_simdat = true
    simdat = decompress_simdat!(out[:simdat])
  else
    has_simdat = false
    y = Matrix{Float64}.(out[:y])
  end
  RESULTS_DIR = join(split(output_path, "/")[1:end-1], "/")
  IMG_DIR = "$(RESULTS_DIR)/img/"
  mkpath(IMG_DIR)

  c = out[:c]
  elbo = out[:metrics][:elbo] / sum(c.N)
  state = out[:state]
  state_hist = out[:state_hist]
  metrics = out[:metrics]

  # ELBO
  dev.pdf("$(IMG_DIR)/elbo.pdf")
  plt.plot(elbo[200:end], xlab="iteration", ylab="ELBO / sum(N)", typ="l")
  dev.dev_off()

  # Metrics
  dev.pdf("$(IMG_DIR)/metrics.pdf")
  plt.par(mfrow=[length(metrics), 1], oma=rcommon.oma_ts(), mar=rcommon.mar_ts())
  for (k, m) in metrics
    plt.plot(metrics[k][1:end]/sum(c.N), xlab="iter", ylab=string(k), typ="l",
             xaxt="n", las=2)
  end
  plt.axis(1)
  plt.par(mfrow=[1, 1], oma=rcommon.oma_default(), mar=rcommon.mar_ts())
  dev.dev_off()

  NSAMPS = 200
  samples = [Cytof5.VB.rsample(state)[2] for n in 1:NSAMPS]
  trace = [Cytof5.VB.rsample(s)[2] for s in state_hist]

  # y_samps
  if has_simdat
    m = [isnan.(yi) for yi in simdat[:y]]
  else
    m = [isnan.(yi) for yi in y]
  end
  if has_simdat
    y_samps = [Tracker.data.(Cytof5.VB.rsample(state, simdat[:y], c)[3]) for n in 1:10]
  else
    y_samps = [Tracker.data.(Cytof5.VB.rsample(state, y, c)[3]) for n in 1:10]
  end
  dev.pdf("$(IMG_DIR)/y_hist.pdf")
  for i in 1:c.I
    plt.hist(vec(y_samps[end][i][m[i]]), xlab="", ylab="", main="", xlim=[-8, 1]);
  end
  dev.dev_off()

  # Prob miss
  ygrid = collect(range(-8, stop=0, length=500))
  pm = hcat([Cytof5.VB.prob_miss(ygrid, c.beta[i]...) for i in 1:c.I]...)
  col=2:4
  dev.pdf("$(IMG_DIR)/prob_miss.pdf")
  plt.matplot(ygrid, pm, xlab="y", ylab="prob. of missing", lw=3, typ="l", lty=1, col=col)
  plt.legend("topright", legend=1:3, col=col, lwd=3)
  dev.dev_off()

  # Z
  if c.use_stickbreak
    Z = [Int.(cumprod(reshape(s.v, 1, c.K)) .> s.H) for s in samples]
  else
    Z = [Int.(reshape(s.v, 1, c.K) .> s.H) for s in samples]
  end
  Z_mean = mean(Z).data
  dev.pdf("$(IMG_DIR)/Z.pdf")
  fZ(z, s_png=1) = plt.abline(h=collect(1:size(z, 1)) .+ .5,
                              v=collect(1:size(z, 2)) .+ .5, col="grey", lwd=s_png)
  cytof3.my_image(Z_mean, xlab="features", ylab="markers",
                  col=cytof3.greys(10), addL=true, f=fZ)
  dev.dev_off()

  # True Z
  if has_simdat
    dev.pdf("$(IMG_DIR)/Z_true.pdf")
    cytof3.my_image(simdat[:Z], xlab="features", ylab="markers", col=cytof3.greys(10), addL=true, f=fZ)
    dev.dev_off()
  end

  # TODO: Add truth

  # W
  W = [s.W.data for s in samples]
  dev.pdf("$(IMG_DIR)/W.pdf")
  for i in 1:c.I
    Wi = hcat([w[i, :] for w in W]...)
    plt.boxplot(Wi');
    if has_simdat
      plt.abline(h=simdat[:W][i, :], lty=2, col="grey")
    end
  end
  dev.dev_off()

  # mu
  mu = hcat([[-cumsum(s.delta0.data); cumsum(s.delta1.data)] for s in samples]...)
  dev.pdf("$(IMG_DIR)/mu.pdf")
  plt.boxplot(mu');
  plt.abline(h=0, v=c.L[0]+.5, col="grey");
  if has_simdat
    plt.abline(h=simdat[:mus][0], lty=2, col="grey");
    plt.abline(h=simdat[:mus][1], lty=2, col="grey");
  end
  dev.dev_off()

  # sig2
  sig2 = hcat([s.sig2.data for s in samples]...)
  dev.pdf("$(IMG_DIR)/sig2.pdf")
  plt.boxplot(sig2');
  if has_simdat
    plt.abline(h=simdat[:sig2], lty=2, col="grey");
  end
  dev.dev_off()

  # eps
  eps = hcat([s.eps.data for s in samples]...)
  dev.pdf("$(IMG_DIR)/eps.pdf")
  plt.boxplot(eps');
  dev.dev_off()

  # v
  v = hcat([s.v.data for s in samples]...)
  dev.pdf("$(IMG_DIR)/v.pdf")
  plt.boxplot(v');
  dev.dev_off()

  # v cumprod
  v_cumprod = hcat([cumprod(s.v.data) for s in samples]...)
  dev.pdf("$(IMG_DIR)/v_cumprod.pdf")
  plt.boxplot(v_cumprod');
  dev.dev_off()

  # alpha
  alpha = vcat([s.alpha.data for s in samples]...);
  dev.pdf("$(IMG_DIR)/alpha.pdf")
  rcommon.plotPost(alpha, main="alpha");
  dev.dev_off()

  ### trace plots ###
  mkpath("$(IMG_DIR)/trace/")

  # Z trace
  if c.use_stickbreak
    Z_trace = [Int.(reshape(cumprod(t.v), 1, c.K) .> t.H).data for t in trace]
  else
    Z_trace = [Int.(reshape(t.v, 1, c.K) .> t.H).data for t in trace]
  end
  dev.pdf("$(IMG_DIR)/trace/Z.pdf")
  for z in Z_trace
    cytof3.my_image(z, xlab="features", ylab="markers")
    plt.abline(h=collect(1:c.J) .+ .5, v=collect(1:c.K) .+ .5, col="grey")
  end
  dev.dev_off()

  # mu_trace
  mu_trace = hcat([[-cumsum(s.delta0.data); cumsum(s.delta1.data)] for s in trace]...)
  dev.pdf("$(IMG_DIR)/trace/mu.pdf")
  plt.matplot(mu_trace', xlab="iter", ylab="mu", typ="l", lwd=2)
  dev.dev_off()

  # sig2_trace
  sig2_trace = hcat([s.sig2.data for s in trace]...)
  dev.pdf("$(IMG_DIR)/trace/sig2.pdf")
  plt.matplot(sig2_trace', xlab="iter", ylab="sig2", typ="l", lwd=2)
  dev.dev_off()

  # eps_trace
  eps_trace = hcat([s.eps.data for s in trace]...)
  dev.pdf("$(IMG_DIR)/trace/eps.pdf")
  plt.matplot(eps_trace', xlab="iter", ylab="eps", typ="l", lwd=2)
  dev.dev_off()

  # W trace
  W_trace = cat([s.W.data for s in trace]..., dims=3)
  for i in 1:c.I
    dev.pdf("$(IMG_DIR)/trace/W$(i).pdf")
    plt.matplot(W_trace[i, :, :]', xlab="iter", ylab="W$(i)", typ="l", lwd=2)
    dev.dev_off()
  end

  # v trace
  v_trace = hcat([s.v.data for s in trace]...)
  dev.pdf("$(IMG_DIR)/trace/v.pdf")
  plt.matplot(v_trace', xlab="iter", ylab="v", typ="l", lwd=2)
  dev.dev_off()

  # alpha trace
  alpha_trace = vcat([s.alpha.data for s in trace]...)
  dev.pdf("$(IMG_DIR)/trace/alpha.pdf")
  plt.plot(alpha_trace, xlab="iter", ylab="alpha", typ="l", lwd=2)
  dev.dev_off()

  ### yZ ###
  if has_simdat
    lam = [sample_lam(state, simdat[:y], c) for b in 1:10]
  else
    lam = [sample_lam(state, y, c) for b in 1:10]
  end
  lam_mode = lam_f(lam, mode)
  W_mean= dropdims(mean(cat(W..., dims=3), dims=3), dims=3)

  mkpath("$(IMG_DIR)/yz")
  s_png = 10
  for i in 1:c.I
    # Yi
    plotpng("$(IMG_DIR)/yz/y$(i)_post.png", s_png)
    lami_est, k_ord = relabel_lam(lam_mode[i], W_mean[i, :])
    if has_simdat
      yi = simdat[:y][i]
    else
      yi = y[i]
    end
    cytof3.my_image(yi[sortperm(lami_est), :], na="black",
                    zlim=[-4,4], col=cytof3.blueToRed(9),
                    f=x->addcut(lami_est, s_png),
                    addL=true, xlab="markers", ylab="cells");
    dev.dev_off()

    # Zi
    dev.pdf("$(IMG_DIR)/yz/Z$(i)_post.pdf")
    k_top = argmax(cumsum(W_mean[i, k_ord]) .> thresh)
    cytof3.my_image(Z_mean[:, k_ord[1:k_top]]', f=z->fZ(z, 1), addL=true,
                    col=cytof3.greys(10), xlab="markers", ylab="cell types")
    dev.dev_off()
  end
end


text:
# Autogenerated wrapper script for qr_mumps_jll for x86_64-apple-darwin-libgfortran5
export libcqrm, libdqrm, libqrm_common, libsqrm, libzqrm

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
using SuiteSparse_jll
using METIS_jll
using SCOTCH_jll
JLLWrappers.@generate_wrapper_header("qr_mumps")
JLLWrappers.@declare_library_product(libcqrm, "@rpath/libcqrm.dylib")
JLLWrappers.@declare_library_product(libdqrm, "@rpath/libdqrm.dylib")
JLLWrappers.@declare_library_product(libqrm_common, "@rpath/libqrm_common.dylib")
JLLWrappers.@declare_library_product(libsqrm, "@rpath/libsqrm.dylib")
JLLWrappers.@declare_library_product(libzqrm, "@rpath/libzqrm.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll, SuiteSparse_jll, METIS_jll, SCOTCH_jll)
    JLLWrappers.@init_library_product(
        libcqrm,
        "lib/libcqrm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdqrm,
        "lib/libdqrm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqrm_common,
        "lib/libqrm_common.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsqrm,
        "lib/libsqrm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzqrm,
        "lib/libzqrm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
apply_filter(u, filter, domain)

Apply `filter` to a spatial function `u` defined on `domain`.
"""
function apply_filter end

apply_filter(u, ::IdentityFilter, ::Domain) = u

function apply_filter(u, filter::TopHatFilter, domain::ClosedIntervalDomain)
    h = filter.width
    a, b = domain.left, domain.right

    ufun = Fun(u, domain.left..domain.right)
    u_int = integrate(ufun)

    # Exact filtered solution
    function ū(x)
        α = max(a, x - h(x))
        β = min(b, x + h(x))
        1 / (β - α) * (u_int(β) - u_int(α))
    end

    ū
end

function apply_filter(u, filter::TopHatFilter, domain::PeriodicIntervalDomain)
    h = filter.width
    a, b = domain.left, domain.right
    L = b - a

    ufun = Fun(u, domain.left..domain.right)
    u_int = integrate(ufun)
    i(α, β) = u_int(β) - u_int(α)

    # Exact filtered solution
    function ū(x)
        α = x - h(x)
        β = x + h(x)

        xₗ = max(a, α)
        xᵣ = min(β, b)

        I = i(xₗ, xᵣ)
        if α < a
            I += i(α + L, b)
        end
        if b < β
            I += i(a, β - L)
        end

        I / (β - α)
    end

    ū
end

function apply_filter(u, filter::ConvolutionalFilter, domain::ClosedIntervalDomain)
    h = filter.width
    G = filter.kernel
    a, b = domain.left, domain.right

    # Exact filtered solution
    function ū(x)
        α = max(a, x - h(x))
        β = min(b, x + h(x))

        ufun = Fun(u, α..β)
        Gfun = Fun(ξ -> G(ξ - x), α..β)

        sum(Gfun * ufun) / sum(Gfun)
    end

    ū
end

function apply_filter(u, filter::ConvolutionalFilter, domain::PeriodicIntervalDomain)
    h = filter.width
    G = filter.kernel
    a, b = domain.left, domain.right
    L = b - a

    # Exact filtered solution
    function ū(x)
        α = x - h(x)
        β = x + h(x)

        u_ext(x) = x < a ? u(x + L) : (x < b ? u(x) : u(x - L))

        ufun = Fun(u, α..β)
        Gfun = Fun(ξ -> G(ξ - x), α..β)

        sum(Gfun * ufun) / sum(Gfun)
    end

    ū
end


text:
function safe_eval(ex::Expr, ed::Acorn.Editor)
    wait_limit = 15.0

    success = false
    expr = :(try catch err return err end)
    expr.args[1] = ex
    inject_yield!(expr)

    try
        t = @async eval(expr)

        tend = time() + wait_limit
        while !istaskdone(t) && (time() <= tend)
            sleep(0.2)
            yield()
        end

        if !istaskdone(t)
            killt = @task try
                Base.throwto(t, InterruptException())
            catch
            end
            try
                schedule(killt, InterruptException(), error = false)
                Acorn.setStatusMessage(ed, "Task killed!")
            catch
            end
        else
            if isa(t.result, Exception)
                throw(t.result)
            else
                success = true
            end
        end

    catch err
        bt = catch_backtrace()
        msg = sprint(showerror, err, bt)
        Acorn.setStatusMessage(ed, msg)
    end

    return success
end

function inject_yield!(ex::Expr)
    if ex.head == :while || ex.head == :for
        # the loop body is a block, so make sure the block ends with a yield
        # i think there *should* ever only be one block, but just to be sure...
        for k in eachindex(ex.args)
            if isa(ex.args[k], Expr) && ex.args[k].head == :block
                push!(ex.args[k].args, :(yield()))
            end
        end
    else
        for k in eachindex(ex.args)
            if isa(ex.args[k], Expr)
                inject_yield!(ex.args[k])
            end
        end
    end
    return ex
end

function evalCommand(str::AbstractString, ed::Acorn.Editor)
    src = "let\n" * str * "\nend"
    exprs = Vector{Expr}()
    k = 1
    while k <= length(src)
        ex, k = Meta.parse(src, k, greedy=true, raise=false)
        ex != nothing && push!(exprs, ex)
    end

    bcontinue = true
    tend = time() + 20
    while !game_over() && bcontinue
        for k in eachindex(exprs)
            bcontinue = safe_eval(exprs[k], ed)
            if !bcontinue
                break
            end
            sleep(0.002)
        end
        bcontinue = bcontinue ? time() <= tend : false
        yield()
    end
end

function commandRun(ed::Acorn.Editor, args::String)
    refresh()
    evalCommand(Acorn.rowsToString(ed.rows), ed)
end

function commandReload(ed::Acorn.Editor, args::String)
    # safe_eval(:(Test.refresh()), ed)
    refresh()
end

function commandLoadLevel(ed::Acorn.Editor, args::String)
    level = ""
    loadLevelCB(ed::Acorn.Editor, buf::String, key::Char) = begin
        level = strip(buf)
    end

    Acorn.editorPrompt(ed, "Level to load: ", callback=loadLevelCB, buf="",
        showcursor=true)

    try
        if occursin(r"\d+", level)
            k = parse(Int, level)
            load_level(k)
            ifile = joinpath(@__DIR__, "..", "programs", @sprintf("level_%02d.jl", k))
        else
            # load an existing level
            if !endswith(level, ".pzl")
                level *= ".pzl"
            end
            pzlfile = joinpath(@__DIR__, "..", "levels", level)
            if isfile(pzlfile)
                ifile = joinfile(@__DIR__, "..", "programs",
                    splitext(level)[1] * ".jl")

                load_level(pzlfile)
            else
                Acorn.setStatusMessage(ed, "[ERROR]: invalid level $(level)")
                return nothing
            end

        end

        Acorn.editorOpen(ed, ifile)
    catch
        Acorn.setStatusMessage(ed, "[ERROR]: invalid level $(level)")
    end
    return nothing
end

function commandLoadPuzzle(ed::Acorn.Editor, args::String)
    level = ""
    loadLevelCB(ed::Acorn.Editor, buf::String, key::Char) = begin
        level = strip(buf)
    end

    Acorn.editorPrompt(ed, "Name of puzzle to load: ", callback=loadLevelCB, buf="",
        showcursor=true)

    if occursin(r"\d+", level)
        k = parse(Int, level)
        ifile = joinpath(@__DIR__, "..", "levels", @sprintf("%02d.pzl", k))
    elseif occursin(r"^[\w\-]+$", level)
        level *= ".pzl"
        ifile = joinpath(@__DIR__, "..", "levels", levl)
    else
        Acorn.setStatusMessage(ed, "[ERROR]: invalid level $(level)")
        return
    end

    Acorn.editorOpen(ed, ifile)
end

function acorn_init()
    Acorn.addCommand(:run, commandRun,
               help="run: run the current script (crtl-x)")

    Acorn.setKeyBinding('x', "run")

    Acorn.addCommand(:reload, commandReload,
               help="reload: reload the current board (crtl-r)")

    Acorn.setKeyBinding('r', "reload")

    Acorn.addCommand(:level, commandLoadLevel,
               help="level: load a level (ctrl-l)")

    Acorn.setKeyBinding('l', "level")

    Acorn.addCommand(:load, commandLoadPuzzle,
               help="load: load a puzzle (ctrl-o)")

    Acorn.rmKeyBinding('o')
    Acorn.setKeyBinding('o', "load")
end


text:
using CalibrationTests
using CalibrationErrors
using Distributions
using StatsBase

using Random
using Statistics
using Test

Random.seed!(1234)

@testset "CalibrationTests" begin
    @testset "Binary trend" begin
        include("binary_trend.jl")
    end

    @testset "Consistency test" begin
        include("consistency.jl")
    end

    @testset "SKCE" begin
        @testset "Asymptotic" begin
            include("skce/asymptotic.jl")
        end
        @testset "Asymptotic block" begin
            include("skce/asymptotic_block.jl")
        end
        @testset "Distribution-free" begin
            include("skce/distribution_free.jl")
        end
    end

    @testset "Asymptotic CME" begin
        include("cme.jl")
    end

    @testset "Deprecations" begin
        include("deprecated.jl")
    end
end


text:
using PortAudio, LibSndFile

stream = PortAudioStream("Microphone (USB Microphone)", 1, 0) # 44100 samples/sec
buf = read(stream, 441000)
save("recorded10sec.wav", buf)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import DirectDepWithoutCompatEntry
import ScientificTypes
import Test

Test.@testset "DirectDepWithoutCompatEntry.jl" begin
    Test.@test DirectDepWithoutCompatEntry.f(1) == 2
end


text:
@testset "Map, reduce, mapreduce, broadcast" begin
    @testset "map and map!" begin
        v1 = @SVector [2,4,6,8]
        v2 = @SVector [4,3,2,1]
        mv = MVector{4, Int}()

        normal_v1 = [2,4,6,8]
        normal_v2 = [4,3,2,1]

        @test map(-, v1) === @SVector [-2, -4, -6, -8]
        @test map(+, v1, v2) === @SVector [6, 7, 8, 9]
        @test map(+, normal_v1, v2) === @SVector [6, 7, 8, 9]
        @test map(+, v1, normal_v2) === @SVector [6, 7, 8, 9]

        map!(+, mv, v1, v2)
        @test mv == @MVector [6, 7, 8, 9]
    end

    @testset "reduce" begin
        v1 = @SVector [2,4,6,8]
        @test reduce(+, v1) === 20
        @test reduce(+, 0, v1) === 20
    end

    @testset "mapreduce" begin
        v1 = @SVector [2,4,6,8]
        v2 = @SVector [4,3,2,1]
        @test mapreduce(-, +, v1) === -20
        @test mapreduce(-, +, 0, v1) === -20
        @test mapreduce(*, +, v1, v2) === 40
        @test mapreduce(*, +, 0, v1, v2) === 40
    end

    @testset "broadcast and broadcast!" begin
        c = 2
        v1 = @SVector [2,4,6,8]
        v2 = @SVector [4,3,2,1]
        mv = MVector{4, Int}()
        mm = MMatrix{4, 2, Int}()
        M = @SMatrix [1 2; 3 4; 5 6; 7 8]

        @test broadcast(-, v1) === map(-, v1)

        @test broadcast(+, v1, c) === @SVector [4, 6, 8, 10]
        @test broadcast(+, v1, v2) === map(+, v1, v2)
        @test broadcast(+, v1, M) === @SMatrix [3 4; 7 8; 11 12; 15 16]

        broadcast!(-, mv, v1)
        @test mv == @MVector [-2, -4, -6, -8]

        broadcast!(+, mv, v1, v2)
        @test mv == @MVector [6, 7, 8, 9]

        broadcast!(+, mm, v1, M)
        @test mm == @MMatrix [3 4; 7 8; 11 12; 15 16]
    end
end


text:
function DiffEqBase.__solve(
  jump_prob::DiffEqBase.AbstractJumpProblem{P},
  alg::DiffEqBase.DEAlgorithm,timeseries=[],ts=[],ks=[],recompile::Type{Val{recompile_flag}}=Val{true};
  kwargs...) where {P,recompile_flag}

  integrator = init(jump_prob,alg,timeseries,ts,ks,recompile;kwargs...)
  solve!(integrator)
  integrator.sol
end

function DiffEqBase.__init(
  _jump_prob::DiffEqBase.AbstractJumpProblem{P},
  alg::DiffEqBase.DEAlgorithm,timeseries=[],ts=[],ks=[],recompile::Type{Val{recompile_flag}}=Val{true};
  callback=nothing, seed = nothing,alias_jump = Threads.threadid() == 1,
  kwargs...) where {P,recompile_flag}

  if alias_jump
    jump_prob = _jump_prob
    reset_jump_problem!(jump_prob,seed)
  else
    jump_prob = resetted_jump_problem(_jump_prob,seed)
  end

  # DDEProblems do not have a recompile_flag argument
  if jump_prob.prob isa DiffEqBase.AbstractDDEProblem
    integrator = init(jump_prob.prob,alg,timeseries,ts,ks;
                      callback=CallbackSet(callback,jump_prob.jump_callback),
                      kwargs...)
  else
    integrator = init(jump_prob.prob,alg,timeseries,ts,ks,recompile;
                      callback=CallbackSet(callback,jump_prob.jump_callback),
                      kwargs...)
  end
end

function resetted_jump_problem(_jump_prob,seed)
  jump_prob = deepcopy(_jump_prob)
  if !isempty(jump_prob.jump_callback.discrete_callbacks)
    if seed === nothing
      Random.seed!(jump_prob.jump_callback.discrete_callbacks[1].condition.rng,seed_multiplier()*rand(UInt64))
    else
      Random.seed!(jump_prob.jump_callback.discrete_callbacks[1].condition.rng,seed)
    end
  end

  if !isempty(jump_prob.variable_jumps)
    @assert jump_prob.prob.u0 isa ExtendedJumpArray
    jump_prob.prob.u0.jump_u .= -randexp.(_jump_prob.rng,eltype(_jump_prob.prob.tspan))
  end
  jump_prob
end

function reset_jump_problem!(jump_prob,seed)
  if seed !== nothing && !isempty(jump_prob.jump_callback.discrete_callbacks)
    Random.seed!(jump_prob.jump_callback.discrete_callbacks[1].condition.rng,seed)
  end

  if !isempty(jump_prob.variable_jumps)
    @assert jump_prob.prob.u0 isa ExtendedJumpArray
    jump_prob.prob.u0.jump_u .= -randexp.(jump_prob.rng,eltype(jump_prob.prob.tspan))
  end
end


text:
@testset "BFFPSV18 algorithm: type constructors" begin

    # interval (1D) decomposition, one variable
    alg = BFFPSV18(δ=1e-3, setrep=Interval, vars=[2], dim=4);
    @test alg.vars == [2]
    @test alg.block_indices == [2]
    @test alg.row_blocks == [[2]]
    @test alg.column_blocks == [[1], [2], [3], [4]]
    @test setrep(alg) == Interval{Float64, IA.Interval{Float64}}

    # interval (1D) decomposition, one variable, without explicitly passing Interval
    alg = BFFPSV18(δ=1e-3, vars=[2], dim=4);
    @test alg.vars == [2]
    @test alg.block_indices == [2]
    @test alg.row_blocks == [[2]]
    @test alg.column_blocks == [[1], [2], [3], [4]]
    @test setrep(alg) == Interval{Float64, IA.Interval{Float64}}

    # interval (1D) decomposition, two variables
    alg = BFFPSV18(δ=1e-3, vars=[1, 4], dim=4);
    @test alg.vars == [1, 4]
    @test alg.block_indices == [1, 4]
    @test alg.row_blocks == [[1], [4]]
    @test alg.column_blocks == [[1], [2], [3], [4]]
    @test setrep(alg) == Interval{Float64, IA.Interval{Float64}}

    # hyperrectangle decomposition
    alg = BFFPSV18(δ=1e-3, setrep=Hyperrectangle, vars=[3], partition=[[1], [2, 3, 4]]);
    @test alg.vars == [3]
    @test alg.block_indices == [2] # variable 3 is in the second block
    @test alg.row_blocks == [[2, 3, 4]]
    @test alg.column_blocks == [[1], [2, 3, 4]]
    @test setrep(alg) == Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}

    # hyperrectangle decomposition, without explicitly passing the set type
    alg = BFFPSV18(δ=1e-3, vars=[3], partition=[[1], [2, 3, 4]]);
    @test alg.vars == [3]
    @test alg.block_indices == [2] # variable 3 is in the second block
    @test alg.row_blocks == [[2, 3, 4]]
    @test alg.column_blocks == [[1], [2, 3, 4]]
    @test setrep(alg) == Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}

    # hyperrectangular decomposition with different block sizes and a single target variable
    alg = BFFPSV18(δ=1e-3, setrep=Hyperrectangle, vars=[25], partition = [1:24, [25], 26:48])
    @test alg.vars == [25]
    @test alg.block_indices == [2]
    @test alg.row_blocks == [[25]]
    @test alg.column_blocks == [1:24, [25], 26:48]

    # hyperrectangular decomposition with different block sizes, with varables belonging to different blocks
    alg = BFFPSV18(δ=1e-3, setrep=Hyperrectangle, vars=[25, 37, 40], partition = [1:24, [25], 26:48])
    @test alg.vars == [25, 37, 40]
    @test alg.block_indices == [2, 3]
    @test alg.row_blocks == [[25], 26:48]
    @test alg.column_blocks == [1:24, [25], 26:48]

    # hyperrectangle decomposition with different variables in the same partition
    alg = BFFPSV18(δ=1e-3, vars=[1, 2], partition=[1:2, 3:5]);
    @test alg.vars == [1, 2]
    @test alg.block_indices == [1]
    @test alg.row_blocks == [[1, 2]]
    @test alg.column_blocks == [[1, 2], [3, 4, 5]]
    @test setrep(alg) == Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}

    # hyperrectangle decomposition with different variables in different partitions
    alg = BFFPSV18(δ=1e-3, vars=[1, 3], partition=[1:2, 3:5]);
    @test alg.vars == [1, 3]
    @test alg.block_indices == [1, 2]
    @test alg.row_blocks == [[1, 2], [3, 4, 5]]
    @test alg.column_blocks == [[1, 2], [3, 4, 5]]
    @test setrep(alg) == Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}
end


text:
using GeometryTypes, StaticArrays, ModernGL
using Transpiler, Visualize
using Visualize: GLRasterization, TextUniforms, get_texture_atlas, Sprite2
using Visualize.GLRasterization: show!, destroy!, swapbuffers!

function vertex_main(vertex)
    return vertex
end

function geometry_main(emit!, vertex_out)
    pos = vertex_out[1]
    quad = Vec4f0(pos[1], pos[2], pos[1] + 0.4f0, pos[2] + 0.4f0)
    v1 = quad[Vec(1, 2)]
    v2 = quad[Vec(1, 4)]
    v3 = quad[Vec(3, 2)]
    v4 = quad[Vec(3, 4)]
    emit!(Vec4f0(v1[1], v1[2], 0f0, 1f0), Vec4f0(1, 0, 0, 1))
    emit!(Vec4f0(v2[1], v2[2], 0f0, 1f0), Vec4f0(1, 0, 0, 1))
    emit!(Vec4f0(v3[1], v3[2], 0f0, 1f0), Vec4f0(1, 0, 0, 1))
    emit!(Vec4f0(v4[1], v4[2], 0f0, 1f0), Vec4f0(1, 0, 0, 1))
    return
end
fragment_main(geom_out) = (geom_out,)

resolution = (1024, 1024)
window = GLFWWindow(Area => resolution, Visualize.Debugging => true)
for event in Visualize.NativeWindowEvents
    add!(window, event)
end
show!(window)
window[Visualize.Open] = true
N = 20

vertices = VertexArray(
    [Point2f0(sin(2pi * (i / N)) , cos(2pi * (i / N))) for i = 1:N],
    face_type = Face{1, OffsetInteger{1, GLint}}
);
draw = GLRasterizer(
    vertices, (),
    vertex_main, fragment_main,
    geometryshader = geometry_main
)

glClearColor(1, 1, 1, 1)

while isopen(window)
    GLFW.PollEvents()
    glViewport(0, 0, widths(window[Area])...)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glUseProgram(draw.program)
    glBindVertexArray(vertices.id)
    GLRasterization.draw_vbo(vertices)
    swapbuffers!(window)
    sleep(0.001)
    yield()
end
destroy!(window)


text:
module DocumentationModule
export Citation, JournalCitation, BookCitation, GenericCitation, summarize, @add_doc
import Base.display
import Base.convert
import Base.promote_rule

"""Abstract class supertyping all objects that are entries of a bibliography (for documenting functions). Each subtype has its corresponding fields and defines its own display function."""
abstract type Citation end

"""
Documentation is a generic type to which all types which are to be used with `@doc` can be automatically converted.
At the moment this means: Any type of Citation and Markdown
"""
type Documentation
  doc::Union{Citation,Base.Markdown.MD}
end

"""Citation class for articles in scientific jounals."""
type JournalCitation <: Citation
  authors::Array{AbstractString,1}
  journalabbr::AbstractString
  vol::UInt
  page::UInt
  year::UInt
end

function display(cite::JournalCitation)
  display(Base.Markdown.parse(""" - $(join(cite.authors,", ")), $(cite.journalabbr) **$(cite.vol)**, *$(cite.page)* ($(cite.year))"""))
end

"""Citation class for complete Books. Refer to BookSectionCitation for more specific parts of a book."""
type BookCitation <: Citation
  authors::Array{AbstractString,1}
  title::AbstractString
  ISBN::AbstractString
end

function display(cite::BookCitation)
  display(Base.Markdown.parse(""" - $(join(cite.authors,", ")): *$(cite.title)* (ISBN $(cite.ISBN))"""))
end

"""Citation class for anything that can not be described by any of the other subclasses of Citation. Make sure none of the others fits the purpose before using this one."""
type GenericCitation <: Citation
  string::AbstractString
end

function display(cite::GenericCitation)
  display(Base.Markdown.parse(""" - $(cite.string)"""))
end

function display(docuArray::Array{Documentation,1})
  local cites = Array{Citation,1}()
  display(Base.Markdown.parse("""
  # Documentation
  """))
  for document in docuArray
    docu = document.doc
    if typeof(docu)==Base.Markdown.MD
      display(docu)
    elseif isa(docu,Citation)
      push!(cites,docu)
    end
  end

  print('\n')

  display(Base.Markdown.parse("""
  # Citations:
  """))
  for citekey in cites
    display(citekey)
  end
end

function display(docuArray::Array{Citation,1})
  display(Base.Markdown.parse("""
  # Citations:
  """))
  for citekey in docuArray
    display(citekey)
  end
end

function display{T<:Citation}(docuArray::Array{T,1})
  display(Base.Markdown.parse("""
  # Citations:
  """))
  for citekey in docuArray
    display(citekey)
  end
end

"""
Taken from base/docs/Docs.jl, this function shows a summary of the fields, type and supertypes of a DataType.
In Docs.jl the function `typesummary` is only meant to be called when no other documentation can be found for it.
If `typesummary` doesn't print "No documentation found" anymore in the future, then this function becomes obsolete.
"""
function summarize(f::DataType)
    parts = [
    """
    **Summary:**
    ```julia
    $(f.abstract ? "abstract" : f.mutable ? "type" : "immutable") $f <: $(supertype(f))
    ```
    """
    ]
    if !isempty(fieldnames(f))
        pad    = maximum([length(string(f)) for f in fieldnames(f)])
        fields = ["$(rpad(f, pad)) :: $(t)" for (f, t) in zip(fieldnames(f), f.types)]
        push!(parts,
        """
        **Fields:**
        ```julia
        $(join(fields, "\n"))
        ```
        """)
    end
    if !isempty(subtypes(f))
        push!(parts,
        """
        **Subtypes:**
        ```julia
        $(join(subtypes(f), "\n"))
        ```
        """)
    end
    Markdown.parse(join(parts, "\n"))
end

# promotion and conversion
Base.promote_rule(::Type{Base.Markdown.MD},::Type{Documentation}) 	= Documentation
Base.promote_rule(::Type{Base.Markdown.MD},::Type{Citation}) 		= Documentation
Base.promote_rule(::Type{Base.Markdown.MD},::Type{GenericCitation}) 	= Documentation
Base.promote_rule(::Type{Base.Markdown.MD},::Type{BookCitation}) 	= Documentation
Base.promote_rule(::Type{Base.Markdown.MD},::Type{JournalCitation}) 	= Documentation
Base.promote_rule(::Type{Documentation},::Type{GenericCitation}) 	= Documentation
Base.promote_rule(::Type{Documentation},::Type{BookCitation}) 		= Documentation
Base.promote_rule(::Type{Documentation},::Type{JournalCitation}) 	= Documentation
Base.convert(::Type{Documentation},md::Base.Markdown.MD) = Documentation(md)
Base.convert(::Type{Documentation},cite::Citation) = Documentation(cite)

macro ignore_warnings(expr)
  quote
    STDERRorig = STDERR
    drop = redirect_stderr()
    result = $expr
    redirect_stderr(STDERRorig)
    return result
  end
end

macro add_doc(newentry,obj)
  _curmod = current_module()
  _obj = :($_curmod.$obj)
  _olddoc = eval(:(@doc $_obj))
  return quote
    if isa($_olddoc, Vector)
      @ignore_warnings @doc [$newentry,$_olddoc...] $_obj  # making sure we promote to the correct type
    else
      @ignore_warnings @doc [$newentry,$_olddoc] $_obj     # making sure we promote to the correct type
    end
  end
end

end # module


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

module TestTypes

using RelationalDatasets
using Test

function test_initialize_relational_dataset()
    data = RelationalDatasets.RelationalDataset((
        ["cancer(alice)."],
        ["cancer(dan)."],
        ["friends(alice,bob)."]
    ))
    return data.pos[1], data.neg[1], data.facts[1]
end

@testset "RelationalDataset Type" begin
    @test test_initialize_relational_dataset() == (
        "cancer(alice).", "cancer(dan).", "friends(alice,bob)."
    )
end

end # module