text:
using Pkg
#  If NetABM is not available then it should be installed via
#  Uncomment if not installed
#  ] add https://github.com/loopulos/NetABM.jl

function useit(list::Array{Symbol})
    installed = [key for key in keys(Pkg.installed())]
    strpackages = @. string(list)
    uninstalled = setdiff(strpackages,installed)

    map(Pkg.add,uninstalled)
    for package ∈ list
        @eval using $package
    end
end

thep = [:LightGraphs, :StatsBase, :Plots, :GraphPlot, :Statistics, :LaTeXStrings, :Bootstrap, :NetABM]
useit(thep)

num_agents = 1000
#  g = erdos_renyi(num_agents,0.4)
g = barabasi_albert(num_agents,9)
agents = [Agent(i) for i = 1:num_agents]
init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
set_coop_agents!(agents;p_cop=0.1)
map(x -> assign_contacts!(g,x), agents);
get_coop!(agents)

S = Array{Float64}(undef,0)
I = Array{Float64}(undef,0)
R = Array{Float64}(undef,0)

for i in 1:100
    next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.7, R=true)
    update_state!(agents)
    push!(S,([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
    push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
    push!(R,([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)
end

plot(1:length(S),S,label="S")
plot!(1:length(S),I,label="I")
plot!(1:length(S),R,label="R")

#  What about levels of cooperation?
#

n_boot = num_agents
some_data = randn(100)
bs1 = bootstrap(std, some_data, BasicSampling(n_boot))
cil = 0.95
bci1 = confint(bs1, BasicConfInt(cil))
II = [zeros(100) for i in 1:100]

g = barabasi_albert(5000,10)
global p = plot()
begin
    Threads.@threads for j in [0.0,0.3,0.6,0.9]
        agents = [Agent(i) for i = 1:5000]
        init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
        set_coop_agents!(agents;p_cop=j)
        map(x -> assign_contacts!(g,x), agents);
        get_coop!(agents)
        I = Array{Float64}(undef,0)
        Threads.@threads for i in 1:100
            next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.7, R=true)
            update_state!(agents)
            push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/5000)
        end
        p = plot!(1:length(I),I,label=string(j))
    end
end
p

#  Changing the average degree of the network

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [2, 6, 10, 20]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,j)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.9)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.5, R=true)
                update_state!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

#  step = 1
#  for step in eachindex(Iall[1])
#      bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
#      bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
#      bci1 = confint(bs1, BasicConfInt(cil))
#      bci2 = confint(bs2, BasicConfInt(cil))
#      push!(Imean,bci1[1][1])
#      push!(Ilow,bci2[1][2])
#      push!(Ihigh,bci2[1][3])
#  end
#
#  Imean
#  Ilow
#  Ihigh
#  plot(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.5)
#  plot(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.5)

#  Changing the cooperation according to neighboors

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true)
                update_state!(agents)
                update_coop!(agents,0.3)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

#  Changing the cooperation according to previous infections

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j = 0.1
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_infections!(agents,i;lrt=false)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sis_adapting_infections.png")

#  Changing the cooperation according to previous infections including lrt

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_infections!(agents,i;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sir_adapting_infections_incl_lrt.png")

#  Changing the cooperation according to neighboors given distance

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [1,2,3]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=0.4)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true)
                update_state!(agents)
                update_coop_distance!(agents,g,j,0.3;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"../figs/sir_adapting_distance.png")

#  Changing the cooperation according to neighboors infections given distance

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0,1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=0.4)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_distance_inf!(agents,g,j,i;lrt=true)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"figs/sis_adapting_distance_infections_neighborhood.png")

#  Changing the cooperation according to first neighboors infection

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
global j
        #  Threads.@threads
global p = plot()
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents]
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1])
            set_coop_agents!(agents;p_cop=0.4)
            set_adapt_agents!(agents;p_cop=j)
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents)
            I = Array{Float64}(undef,0)
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false)
                update_state!(agents)
                update_coop_distance_inf!(agents,g,1,i;lrt=false)
                get_coop!(agents)
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
            end
            push!(Iall,I)
            #  p = plot!(1:length(I),I,label=string(j))
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

p=plot!(xlims=(0,100),ylims=(0,1))
p=plot!(xlabel = L"t_i", ylabel = L"I/N")
savefig(p,"figs/sis_adapting_first_neighboors.png")

#  Changing the cooperation using infection counter
#  first infection.

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
global adoptprob = (1,[0.5,0.3])
        #  Threads.@threads
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.4);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=false);
                update_state!(agents);
                update_coop_given_distance!(agents,g,adoptprob[1],1,adoptprob[2];lrt=false);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sis_adapting_d1.png")
p

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 10
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
global adoptprob = (1,[0.5,0.3])
        #  Threads.@threads
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7, 1]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.4);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;inf_prob=0.1, rec_prob=0.1, coop_red=0.8, R=true);
                update_state!(agents);
                update_coop_given_distance!(agents,g,adoptprob[1],1,adoptprob[2];lrt=false);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Imean .- Ilow, Ihigh .- Imean), fillalpha=0.35,label=string(j))
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sir_adapting_d1.png")
p

###################################################
###################################################
###################################################
###################################################
# Adaptive responses with risk attitudes
# first let's just increase the effectiveness
# given the times they got infected themselves
# this would be an equivalent to increase biosecurity

num_agents = 1000
global n_boot = 100
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
begin
    for j in [0.1, 0.4, 0.7]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            #  the_selected = rand(1:num_agents)
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            #  init_demographics!(agents[setdiff(1:num_agents,the_selected)];states=["S","I"],initial=[1.,0.]);
            #  init_demographics!([agents[the_selected]];states=["S","I"],initial=[0.,1.]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.1, R=true);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
            end
            push!(Iall,I)
        end
        #  end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sir_attitudes_1I.png")
p

###################################################
###################################################
###################################################
###################################################
# This is the same as the previous but looking at the
# mean ρ of each group. Didn't want to mess with the
# already working code.

num_agents = 1000
global Imean = Array{Float64}(undef,0)
global Ilow = Array{Float64}(undef,0)
global Ihigh = Array{Float64}(undef,0)
global Iall = Array{Array}(undef,0)
global n_boot = 100
global cil = 0.95
#  global d = 0
#  global probs = [1]
global j
#  global adoptprob = (1,[0.5,0.3])
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
        global rtall = Array{Float64}(undef,0);
        global lrtall = Array{Float64}(undef,0);
        global raall = Array{Float64}(undef,0);
        global lraall = Array{Float64}(undef,0);
begin
    for j in [0.4]
        global Imean = Array{Float64}(undef,0)
        global Ilow = Array{Float64}(undef,0)
        global Ihigh = Array{Float64}(undef,0)
        global Iall = Array{Array}(undef,0)
        Threads.@threads for a in 1:100
            g = barabasi_albert(num_agents,10)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.7,0.3]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            rt = Array{Float64}(undef,0);
            ra = Array{Float64}(undef,0);
            lrt = Array{Float64}(undef,0);
            lra = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.1, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                #  agents[1].adapter = true
                #  agents[3].counter = 1
                #  agents[1].coop_effect
                #  update_single_effect_distance!(agents,g,3,0,1,adoptstep[2])
                #  agents[1].coop_effect
                push!(rt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && !ag.adapter]));
                push!(lrt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && ag.adapter]));
                push!(ra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && !ag.adapter]));
                push!(lra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && ag.adapter]));
            end
            push!(rtall,I)
            push!(lrtall,I)
            push!(raall,I)
            push!(lraall,I)
        end
        end
        #  for step in eachindex(Iall[1])
        #      bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
        #      bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
        #      bci1 = confint(bs1, BasicConfInt(cil))
        #      bci2 = confint(bs2, BasicConfInt(cil))
        #      push!(Imean,bci1[1][1])
        #      push!(Ilow,bci2[1][2])
        #      push!(Ihigh,bci2[1][3])
        #  end
        #  p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
        #  p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
end
p

savefig(p,"figs/sis_attitudes_1I.png")
p

####################################
# Memory test
#
using NetABM
using LightGraphs
using StatsBase
using Plots
using Statistics
using Bootstrap
using LaTeXStrings

global num_agents = 1000
global n_boot = 100
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")

@time let
    #  for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    for j in [0.1,0.3,0.7]
        Imean = Array{Float64}(undef,0)
        Ilow = Array{Float64}(undef,0)
        Ihigh = Array{Float64}(undef,0)
        Iall = Array{Array}(undef,0)
        for a in 1:100
            g = barabasi_albert(num_agents,3)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            I = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                #  GC.gc()
            end
            push!(Iall,I)
            GC.gc()
        end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
            #  GC.gc()
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
    GC.gc()
end
p
savefig(p,"figs/sis_attitudes_01.png")

#  Track number of infections

global num_agents = 1000
global n_boot = 1000
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global maxI = Array{Float64}(undef,0)
global sumI = []
#  global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")

@time let
    #  for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    for copera in 0.0:0.1:1
        for adap in 0.0:0.1:1
            global meantotI = Array{Float64}(undef,0)
            global Imean = Array{Float64}(undef,0)
            global Ilow = Array{Float64}(undef,0)
            global Ihigh = Array{Float64}(undef,0)
            global Itotall = Array{Float64}(undef,0)
            global Iall = Array{Array}(undef,0)
            global Sall = Array{Array}(undef,0)
            global Rall = Array{Array}(undef,0)
            for a in 1:50#100
                g = barabasi_albert(num_agents,3)
                agents = [Agent(i) for i = 1:num_agents];
                init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
                set_coop_agents!(agents;p_cop=copera);
                set_adapt_agents!(agents;p_cop=adap);
                map(x -> assign_contacts!(g,x), agents);
                #  I = Array{Float64}(undef,0);
                tau = []
                num_S = []
                num_I = []
                num_R = []
                #  global oldS = 1.
                global oldS = ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents
                #  global oldI = 0
                global oldI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                #  global totI = 0
                global totI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                for i in 1:100
                    #  push!(tau, params.now_t)
                    push!(num_S, ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
                    push!(num_I, ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
                    push!(num_R, ([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)

                    next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                    update_state!(agents);
                    update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);

                    newI = oldS - num_S[end]
                    global totI = totI + newI
                    oldS = num_S[end]

                    #  push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                end
                #  push!(Iall,I)
                push!(Iall,num_I)
                push!(Itotall,totI)
                GC.gc()
            end
            for step in eachindex(Iall[1])
                #  bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                #  bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                #  bci1 = confint(bs1, BasicConfInt(cil))
                #  bci2 = confint(bs2, BasicConfInt(cil))
                #  push!(Imean,bci1[1][1])
                #  push!(Ilow,bci2[1][2])
                #  push!(Ihigh,bci2[1][3])
                bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                bci1 = confint(bs1, BasicConfInt(cil))
                bci2 = confint(bs2, BasicConfInt(cil))
                push!(Imean,bci1[1][1])
                push!(Ilow,bci2[1][2])
                push!(Ihigh,bci2[1][3])
            end
            bs1 = bootstrap(mean, Itotall,BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            push!(sumI,bci1[1][1])
            push!(maxI,maximum(Imean))
            #  p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
        end
        GC.gc()
    end
end
p
savefig(p,"figs/sis_attitudes_01.png")

using Colors
] add ColorSchemes
] add PlotThemes
using ColorSchemes
using PlotThemes
theme(:default)
Plots.showtheme(:sand)
pyplot()
gr()

plt = heatmap(
        #  xlabel = L"p_c",
        xlabel = "Cooperative Agents",
        #  ylabel = L"\tau_l",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(maxI, (11, 11))'
#  plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :YlOrRd_9)
plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :algae,dpi=300)
#  plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :thermal)
#  plt = heatmap!(xs, ys, z, aspect_ratio = 1)
plt
png(plt,"figs/heatmap_scaled_maximum.png")
savefig(plt,"figs/heatmap_scaled_maximum.pdf")

plt = heatmap(
        #  xlabel = L"p_c",
        xlabel = "Cooperative Agents",
        #  ylabel = L"\tau_l",
        ylabel = "Adaptive Agents",
        title = "Total Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(sumI, (11, 11))'
plt = heatmap!(xs, ys, z, aspect_ratio = 1,c = :algae,dpi=300)
#  plt = heatmap!(xs, ys, z, aspect_ratio = 1)
plt
png(plt,"figs/heatmap_scaled_total.png")
savefig(plt,"figs/heatmap_scaled_total.pdf")

global rtall = Array{Array}(undef,0);
global lrtall = Array{Array}(undef,0);
global raall = Array{Array}(undef,0);
global lraall = Array{Array}(undef,0);
global adoptstep = (1,[0.1,0.05])

global num_agents = 1000
global n_boot = 100
global cil = 0.95
global j
global adoptstep = (0,[0.1])
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
global p = plot(xlims=(0,100),ylims=(0,1),xlabel = L"t_i", ylabel = L"I/N")
@time let
    #  for j in [0.0, 0.1, 0.4, 0.7, 1.0]
    for j in [0.5]
        Imean = Array{Float64}(undef,0)
        Ilow = Array{Float64}(undef,0)
        Ihigh = Array{Float64}(undef,0)
        Iall = Array{Array}(undef,0)
        for a in 1:100
            g = barabasi_albert(num_agents,3)
            agents = [Agent(i) for i = 1:num_agents];
            init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
            set_coop_agents!(agents;p_cop=0.5);
            set_adapt_agents!(agents;p_cop=j);
            map(x -> assign_contacts!(g,x), agents);
            #  get_coop!(agents);
            I = Array{Float64}(undef,0);
            rt = Array{Float64}(undef,0);
            ra = Array{Float64}(undef,0);
            lrt = Array{Float64}(undef,0);
            lra = Array{Float64}(undef,0);
            for i in 1:100
                next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=false);
                update_state!(agents);
                update_effect_given_distance!(agents,g,adoptstep[1],1,adoptstep[2]);
                #  get_coop!(agents);
                push!(I,([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents);
                push!(rt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && !ag.adapter]));
                push!(lrt,mean([ag.coop_effect for ag in agents if ag.attitude == "rt" && ag.adapter]));
                push!(ra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && !ag.adapter]));
                push!(lra,mean([ag.coop_effect for ag in agents if ag.attitude == "ra" && ag.adapter]));
            end
            push!(Iall,I)
            push!(rtall,rt)
            push!(lrtall,lrt)
            push!(raall,ra)
            push!(lraall,lra)
            GC.gc()
        end
        for step in eachindex(Iall[1])
            bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
            bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            bci2 = confint(bs2, BasicConfInt(cil))
            push!(Imean,bci1[1][1])
            push!(Ilow,bci2[1][2])
            push!(Ihigh,bci2[1][3])
            #  GC.gc()
        end
        p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label=string(j))
    end
    GC.gc()
end

Imean = Array{Float64}(undef,0)
Ilow = Array{Float64}(undef,0)
Ihigh = Array{Float64}(undef,0)
Iall = Array{Array}(undef,0)
for step in eachindex(lraall[1])
    bs1 = bootstrap(mean,getindex.(lraall,step),BasicSampling(n_boot))
    bs2 = bootstrap(std,getindex.(lraall,step),BasicSampling(n_boot))
    bci1 = confint(bs1, BasicConfInt(cil))
    bci2 = confint(bs2, BasicConfInt(cil))
    push!(Imean,bci1[1][1])
    push!(Ilow,bci2[1][2])
    push!(Ihigh,bci2[1][3])
end
pyplot()
p=plot!(1:length(Imean), Imean, ribbon = (Ilow, Ihigh), fillalpha=0.35,label="lrt")
p

p
#  savefig(p,"figs/sis_attitudes_01.png")

ext: jl


text:
function init_types(model::Model,geq_expressions::Array=[],leq_expressions::Array=[],eq_expressions::Array=[],constraint_per_iteration::Int = 1)
    heristic = data(model,geq_expressions,leq_expressions,eq_expressions)

    return heristic, parameters
end

function step_solve!(heristic::data, parameters::options)
    optimize!(heristic.model)
    nothing
    #return heristic, parameters
end

function evaluate_violations!(heristic::data, parameters::options)
    geq_error = - value.(heristic.geq_expressions)
    leq_error = value.(heristic.leq_expressions)
    eq_error = abs.(value.(heristic.eq_expressions))

    values, indices = findextrema(geq_error,parameters.constraint_per_iteration)

    for i in 1:indices
        con = @build_constraint(heristic.geq_expressions[i] >= 0)
        add_constraint(heristic.model, con)

    end
    nothing
end

function findextrema(v::Array,n::Int)
    idx = sortperm(v; rev=true)
    return v[idx][1:n], idx[1:n]
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using InteractiveDynamics
using DynamicalSystems
using GLMakie

using OrdinaryDiffEq: Tsit5
ds = Systems.lorenz()

u0s =  [[10,20,40.0] .+ rand(3) for _ in 1:6]

diffeq = (alg = Tsit5(), dtmax = 0.01)

fig, obs = interactive_evolution(
    ds, u0s; tail = 1000, diffeq, colors = to_color.(COLORS)
)

ax3D = content(fig[1,1])
ax3D.elevation = 0.24269908169872415
ax3D.azimuth = 5.645530633326986

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 10)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Test, Distributions, ForwardDiff, ReverseDiff, KernelGoodnessOfFit

# setup the problem
rbf_ = GaussianRBF(1.0)

μ = 0.0
σ² = 1.0
p_ = MultivariateNormal([μ, 0], [σ² 0; 0 σ²])  # true
q_ = MultivariateNormal([μ, 5.0], [σ² 0; 0 σ²])  # false

vs_ = [0.0 ; 0.0]
vs_ = reshape(vs_, 2, 1)

n_ = 400
xs_ = rand(p_, n_)

# new vs old: verify that new implementation is correct
@testset "fssd vs fssd_old: eq" begin
    @test KernelGoodnessOfFit.fssd(rbf_, q_, xs_, vs_) == KernelGoodnessOfFit.fssd_old(rbf_, q_, xs_, vs_)
    @test KernelGoodnessOfFit.fssd(rbf_, p_, xs_, vs_) == KernelGoodnessOfFit.fssd_old(rbf_, p_, xs_, vs_)
end

# ensure new method is faster :)
@testset "fssd vs fssd_old: perf" begin
    t1 = @elapsed KernelGoodnessOfFit.fssd(rbf_, q_, xs_, vs_)
    t2 = @elapsed KernelGoodnessOfFit.fssd_old(rbf_, q_, xs_, vs_)

    @test t1 < t2
end

# sizes: testing of sizes
@testset "shapes" begin
    # setup for tests
    Ξ_ = KernelGoodnessOfFit.compute_Ξ(rbf_, q_, xs_, vs_)

    J, = size(Ξ_)
    d, n = size(Ξ_[1])

    @test (d, J) == size(vs_)
    @test (d, n) == size(xs_)

    τ = KernelGoodnessOfFit.τ_from_Ξ(Ξ_)
    @test size(τ) == (d * J, n)

    μ, Σ = KernelGoodnessOfFit.Σₚ(τ)
    @test size(μ) == (d * J, 1)
    @test size(Σ) == (d * J, d * J)

    σ = KernelGoodnessOfFit.σ²_H₁(μ, Σ)
    @test size(σ) == ()
end

# gradients: verify gradients
@testset "gradient shapes" begin
    Δv = ForwardDiff.gradient(v -> KernelGoodnessOfFit.fssd_H₁_opt_factor(rbf_, q_, xs_, v), vs_)
    @test size(Δv) == size(vs_)

    Δv = ReverseDiff.gradient(v -> KernelGoodnessOfFit.fssd_H₁_opt_factor(rbf_, q_, xs_, v), vs_)
    @test size(Δv) == size(vs_)

    # fails for some reason; but does not fail when I actually run the test...interesting, eh?
    # Δσₖ = ForwardDiff.derivative(σₖ -> fssd_H₁(GaussianRBF(σₖ), q_, xs_, vs_), 1.0)
end

# values: testing specific values
@testset "specific values" begin
    V = reshape([0; 0], (2, 1))

    X = [1.0 2.0 3.0; 2.0 3.0 4.0]

    Ξ_1_true = [
        -0.116086 -0.00425237 -1.58109e-5;
        0.0580429 -0.00106309 -7.90543e-6
    ]
    Ξ_ = KernelGoodnessOfFit.compute_Ξ(rbf_, q_, X, V)

    @test Ξ_[1] ≈ Ξ_1_true atol=0.01

    fssd_true = 0.0004333865124975432

    @test KernelGoodnessOfFit.fssd(Ξ_) * size(X)[2] ≈ fssd_true

    σ₁² = 2.8732542480282987e-5

    τ_ = KernelGoodnessOfFit.τ_from_Ξ(Ξ_)
    μ_, Σ_ = KernelGoodnessOfFit.Σₚ(τ_)

    @test KernelGoodnessOfFit.σ²_H₁(μ_, Σ_) ≈ σ₁²
end

# optimizing the power
@testset "optimize_power" begin
    # because of compilation and stuff, running them in sequence is heavily biased towards
    # the ones being run last
    KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, q_; diff = :difference)
    KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, q_; diff = :backward)
    KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, q_; diff = :forward)

    μs = 1:10

    t_difference = mean([@elapsed KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, MultivariateNormal([μ, 5.0], [σ² 0; 0 σ²]); diff = :difference) for μ ∈ μs])
    t_forward = mean([@elapsed KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, MultivariateNormal([μ, 5.0], [σ² 0; 0 σ²]); diff = :forward) for μ ∈ μs])
    t_backward = mean([@elapsed KernelGoodnessOfFit.optimize_power(rbf_, vs_, xs_, MultivariateNormal([μ, 5.0], [σ² 0; 0 σ²]); diff = :backward) for μ ∈ μs])

    @test t_backward ≥ 0.0
    @test t_forward ≥ 0.0
    @test t_forward < t_difference
end


text:
ENV["PHAGEDATAPATH"]="/home/matteo/poliTo/Tesi/Data"

using UAF2, PhageFields, PhageData, NLopt
using Test

boyer = import_dataset(:boyer)

boyer = subdata(boyer, not_disappearing(boyer))
boyer = add_pseudocounts(boyer, 0.5)

epistasis = EpistasisMu{boyer.A, boyer.L}()
model = Model(epistasis, boyer.T)

ntrain = round(Int, boyer.S*0.8)
train, test = randtrain(boyer, ntrain)

starting_point!(model, train)
prior = set_prior(100.0, model)
opt = set_optimizer(:LD_MMA, model, train; ftol_rel = 1e-7, Δλ = 1e-10)

get_lambda(m ,p ,d) = m.λ
get_mu(m, p, d) = [m.x[end]];

model, ret, logP, properties = learn!(model, prior, opt, train, monitor = [get_lambda, get_mu])


text:
# ---
# title: 1681. Minimum Incompatibility
# id: problem1681
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Backtracking, Greedy
# link: <https://leetcode.com/problems/minimum-incompatibility/description/>
# hidden: true
# ---
#
# You are given an integer array `nums` and an integer `k`. You are asked to
# distribute this array into `k` subsets of **equal size** such that there are
# no two equal elements in the same subset.
#
# A subset's **incompatibility** is the difference between the maximum and
# minimum elements in that array.
#
# Return _the **minimum possible sum of incompatibilities** of the _`k` _subsets
# after distributing the array optimally, or return_`-1` _if it is not
# possible._
#
# A subset is a group integers that appear in the array with no particular
# order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,1,4], k = 2
#     Output: 4
#     Explanation: The optimal distribution of subsets is [1,2] and [1,4].
#     The incompatibility is (2-1) + (4-1) = 4.
#     Note that [1,1] and [2,4] would result in a smaller sum, but the first subset contains 2 equal elements.
#
# **Example 2:**
#
#
#
#     Input: nums = [6,3,8,1,3,1,2,2], k = 4
#     Output: 6
#     Explanation: The optimal distribution of subsets is [1,2], [2,3], [6,8], and [1,3].
#     The incompatibility is (2-1) + (3-2) + (8-6) + (3-1) = 6.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [5,3,3,6,3,3], k = 3
#     Output: -1
#     Explanation: It is impossible to distribute nums into 3 subsets where no two elements are equal in the same subset.
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= nums.length <= 16`
#   * `nums.length` is divisible by `k`
#   * `1 <= nums[i] <= nums.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
type BartProbability #Records the proposal probabilities for moves
  node_grow_prune::Float64
  change_decision_rule::Float64
  swap_decision_rule::Float64

  function BartProbability(n,c,s)
    assert(n+c+s==1)
    new(n,c,s)
  end
end

BartProbability() = BartProbability(0.5,0.4,0.1)

type BartOptions
  num_trees::Int #number of additive trees
  num_draws::Int #number of posterior draws
  num_burn_in::Int
  num_thinning::Int
  alpha::Float64 #calculating nonterminal probability
  beta::Float64
  k::Int
  proposal_probabilities::BartProbability
end

function bart_options(;num_trees::Int=200,num_draws::Int=1000,
                     num_burn_in::Int=200,num_thinning::Int=10,alpha::Float64=0.95,beta::Float64=2.0,
                     k::Int=2,proposal_probabilities=BartProbability())
   BartOptions(num_trees,num_draws,num_burn_in,num_thinning,alpha,beta,k,proposal_probabilities)
end

type BartStateParameters
  sigma::Float64
  sigma_mu::Float64
  nu::Int64
  lambda::Float64
end

type BartTree
  tree::DecisionTree
end

@make_tree_type(BartTree,DecisionNode)

type BartState
  trees::Vector{BartTree}
  parameters::BartStateParameters
end

type BartAdditiveTree
  trees::Vector{BartTree}
end

type Bart
  bart_additive_trees::Vector{BartAdditiveTree} #store the posterior draw of additive tree
  y_min::Float64
  y_max::Float64
  options::BartOptions
end

type BartLeaf <: DecisionLeaf
  value::Float64 #\mu_{ij}
  residual_mean::Float64
  residual_sigma::Float64
  leaf_data_indices::Vector{Int} #all the training data contained in this leaf

  function BartLeaf(value::Float64,residual::Vector{Float64},leaf_data_indices::Vector{Int})
        #here, we don't allow leaves with no training data
        residual_mean = mean(residual[leaf_data_indices])
        residual_sigma = sqrt(mean((residual[leaf_data_indices].-residual_mean).^2))
        new(value,residual_mean,residual_sigma,leaf_data_indices)
  end
end

function node_nonterminal(depth::Int,alpha::Float64=0.95,beta::Float64=2.0)
  alpha*depth^(-beta)
end

node_nonterminal(depth::Int,bartoptions::BartOptions) = node_nonterminal(depth,bartoptions.alpha,bartoptions.beta)

function log_tree_prior(branch::DecisionBranch,depth::Int,bartoptions::BartOptions)
    prior = log(node_nonterminal(depth,bartoptions))
    #Here, we don't consider the probability of selecting available features
    #since we always assume that all the variables are available
    #We only consider the probability of selecting available values
    #Here, it's possible that length(train_data_indices(branch))==0,then prior = Inf
    prior -= log(length(train_data_indices(branch)))
    prior += log_tree_prior(branch.left,depth+1,bartoptions)
    prior += log_tree_prior(branch.right,depth+1,bartoptions)
    prior
end

log_tree_prior(leaf::BartLeaf,depth::Int,bartoptions::BartOptions) = log(1-node_nonterminal(depth,bartoptions))

node_loglikelihood(branch::DecisionBranch,bart_state::BartState) = node_loglikelihood(branch.left,bart_state)+node_loglikelihood(branch.right,bart_state)

function node_loglikelihood(leaf::BartLeaf,bart_state::BartState)
    n = length(leaf.leaf_data_indices)
    #Here, we don't allow leaf with no train_data_indices
    if n==0
        return Inf
    end
    loglikelihood =
      0.5*(log(bart_state.parameters.sigma^2)-log(bart_state.parameters.sigma^2+n*bart_state.parameters.sigma_mu^2))
    loglikelihood -= 0.5*n*leaf.residual_sigma^2/bart_state.parameters.sigma^2
    loglikelihood -= 0.5*n*leaf.residual_mean^2/(n*bart_state.parameters.sigma_mu^2+bart_state.parameters.sigma^2)
    loglikelihood
end

function update_tree!(bart_state::BartState,tree::BartTree,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
    update_probability = rand()
    prob = bartoptions.proposal_probabilities #probabilities for posterior proposal moves

    if update_probability<prob.node_grow_prune
      updated = node_grow_prune!(bart_state,tree,x,residual,bartoptions)
    elseif update_probability<prob.node_grow_prune+prob.change_decision_rule
      updated = change_decision_rule!(bart_state,tree,x,residual,bartoptions)
    else
      updated = swap_decision_rule!(bart_state,tree,x,residual,bartoptions)
    end
    #after we sample T_j, now we sample M_j given T_j
    #no matter we update updated or not, I think we should update M_j(leaf_values)
    update_leaf_values!(tree,bart_state.parameters)
    updated
end

function node_grow_prune!(bart_state::BartState,tree::BartTree,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
    #if the tree only has one root,then we must grow
    probability_grow = typeof(tree.tree.root)==BartLeaf?1.0:0.5
    u = rand()
    if u<probability_grow
      updated = node_grow!(bart_state,tree,probability_grow,x,residual,bartoptions)
    else
      probability_prune = 1-probability_grow
      updated = node_prune!(bart_state,tree,probability_prune,x,residual,bartoptions)
    end
    updated
end

function node_grow!(bart_state::BartState,tree::BartTree,probability_grow::Float64,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
    #first we need to sample a leave to grow new leaves
    leaf_nodes = leaves(tree)
    number_leaves = length(leaf_nodes)
    i = rand(1:number_leaves)
    leaf = leaf_nodes[i]
    leaf_indices = leaf.leaf_data_indices

    #it's impossible to split a leaf with no data points
    #it's feasible to split a leaf with one data points
    if length(leaf_indices)==0
      return false
    end

    leaf_depth = depth(tree,leaf)
    leaf_loglikelihood = node_loglikelihood(leaf,bart_state)
    leaf_nonterminal = node_nonterminal(leaf_depth,bartoptions)

    #uniformly sample a splitting feature from the available predictors

    #split_feature = rand(1:size(x,2))
     split_feature = rand(1:size(x,1))
    #uniformly sample a splitting location from the available locations
    #feature = x[leaf_indices,split_feature]
    feature = vec(x[split_feature,leaf_indices])
    split_location = rand(1:length(feature))
    split_value = sort(feature)[split_location]

    #It's possible that length(left_indices) =0 or length(right_indices)=0
    left_indices = leaf_indices[feature.<=split_value]
    right_indices = leaf_indices[feature.>split_value]
    if length(left_indices)*length(right_indices)==0
       return false
    end
    left_leaf_nonterminal = length(left_indices)>0?node_nonterminal(leaf_depth+1,bartoptions):0.0
    right_leaf_nonterminal = length(right_indices)>0?node_nonterminal(leaf_depth+1,bartoptions):0.0

    #construct the new branch and new leaf
    left_leaf = BartLeaf(0.0,residual,left_indices)
    right_leaf = BartLeaf(0.0,residual,right_indices)
    branch = DecisionBranch(split_feature,split_value,left_leaf,right_leaf)
    left_leaf_loglikelihood = node_loglikelihood(left_leaf,bart_state)
    right_leaf_loglikelihood = node_loglikelihood(right_leaf,bart_state)

    if left_leaf_loglikelihood==Inf || right_leaf_loglikelihood==Inf
      return false
    end

    #calculate the number of internal nodes with two terminal nodes of the original tree, which is the
    #not_grand_branches
    not_grand_branch_nodes = not_grand_branches(tree)
    number_not_grand_branch_nodes = length(not_grand_branch_nodes)

    parent_leaf = parent(tree,leaf) #it's posssible that there is no parent
    if parent_leaf!=nothing
      if typeof(parent_leaf.left)==BartLeaf && typeof(parent_leaf.right)==BartLeaf
         number_not_grand_branch_nodes = number_not_grand_branch_nodes
      else
      number_not_grand_branch_nodes += 1
      end
    else
      number_not_grand_branch_nodes += 1
    end

    #the probability of pruning a node from the proposal tree(after grow)
    probability_prune = 0.5

    alpha = (1-left_leaf_nonterminal)*(1-right_leaf_nonterminal)*leaf_nonterminal*number_leaves*probability_prune
    alpha /= (1-leaf_nonterminal)*number_not_grand_branch_nodes*probability_grow
    alpha *= exp(left_leaf_loglikelihood+right_leaf_loglikelihood-leaf_loglikelihood)

    #decide whether or not to grow two leaves
    if rand()<alpha
       if parent_leaf==nothing
         tree.tree.root = branch
       elseif leaf==parent_leaf.left
         parent_leaf.left = branch
       else
         parent_leaf.right = branch
       end
       updated = true
    else
       updated = false
    end
    updated
end

function node_prune!(bart_state::BartState,tree::BartTree,probability_prune::Float64,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
  if typeof(tree.tree.root.left)==BartLeaf && typeof(tree.tree.root.right)==BartLeaf
     probability_grow = 1.0
  else
     probability_grow = 0.5
  end

  not_grand_branch_nodes = not_grand_branches(tree)
  number_not_grand_branch_nodes = length(not_grand_branch_nodes)
  i = rand(1:number_not_grand_branch_nodes)
  not_grand_branch = not_grand_branch_nodes[i]

  leaf_nodes = leaves(tree)
  number_leaves = length(leaf_nodes)

  not_grand_branch_depth = depth(not_grand_branch)
  not_grand_branch_nonterminal = node_nonterminal(not_grand_branch_depth,bartoptions)
  left_leaf_nonterminal = length(not_grand_branch.left.leaf_data_indices)>0?node_nonterminal(not_grand_branch_depth+1,bartoptions):0.0
  right_leaf_nonterminal = length(not_grand_branch.right.leaf_data_indices)>0?node_nonterminal(not_grand_branch_depth+1,bartoptions):0.0

  not_grand_branch_loglikelihood = node_loglikelihood(not_grand_branch,bart_state)
  #here, I construct a new leaf by combining the information in the left and right leaves of the not_grand_branch
  not_grand_branch_data_indices = vcat(not_grand_branch.left.leaf_data_indices,not_grand_branch.right.leaf_data_indices)
  new_leaf = BartLeaf(0.0,residual,not_grand_branch_data_indices)
  proposal_new_leaf_loglikelihood= node_loglikelihood(new_leaf,bart_state)

  alpha = (1-not_grand_branch_nonterminal)*number_not_grand_branch_nodes*probability_grow
  alpha /= (1-left_leaf_nonterminal)*(1-right_leaf_nonterminal)*not_grand_branch_nonterminal*(number_leaves-1)*probability_prune
  alpha *= exp(proposal_new_leaf_loglikelihood-not_grand_branch_loglikelihood)

  parent_not_grand_branch = parent(tree,not_grand_branch)

  if rand() < alpha
    if parent_not_grand_branch == nothing
       tree.tree.root = new_leaf
    elseif not_grand_branch==parent_not_grand_branch.left
          parent_not_grand_branch.left = new_leaf
    else
          parent_not_grand_branch.right = new_leaf
    end
    updated = true
  else
    updated = false
  end
  updated
end

train_data_indices(leaf::BartLeaf) = leaf.leaf_data_indices

function train_data_indices(branch::DecisionBranch)
  function train_data_indices!(branch::DecisionBranch,indices::Vector{Int})
    train_data_indices!(branch.left,indices)
    train_data_indices!(branch.right,indices)
  end

  function train_data_indices!(leaf::BartLeaf,indices::Vector{Int})
    for i in 1:length(leaf.leaf_data_indices)
        push!(indices,leaf.leaf_data_indices[i])
    end
  end

  indices = Int[]
  train_data_indices!(branch,indices)
  indices
end

function tree_adjust!(parent::DecisionBranch,leaf::BartLeaf,x::Matrix{Float64},residual::Vector{Float64},branch_indices::Vector{Int},left::Bool)
   valid_split=true
   if left
     parent.left = BartLeaf(leaf.value,residual,branch_indices)
   else
     parent.right = BartLeaf(leaf.value,residual,branch_indices)
   end
   if length(branch_indices)==0
     valid_split = false
   end
   valid_split
end
#it might be possible that length(branch_indices)==0
function tree_adjust!(parent::DecisionBranch,branch::DecisionBranch,x::Matrix{Float64},residual::Vector{Float64},branch_indices::Vector{Int},left::Bool)
  tree_adjust!(branch,x,residual,branch_indices)
end

function tree_adjust!(branch::DecisionBranch,x::Matrix{Float64},residual::Vector{Float64},branch_indices::Vector{Int64})
    valid_split = true
    if length(branch_indices)==0
       left_data_indices = branch_indices
       right_data_indices = branch_indices
       valid_split = false
    else
        #feature = x[branch_indices,branch.feature]
        feature = vec(x[branch.feature,branch_indices])
        left_data_indices = branch_indices[feature.<=branch.value]
        right_data_indices = branch_indices[feature.>branch.value]
    end

    valid_split = tree_adjust!(branch,branch.left,x,residual,left_data_indices,true)
    valid_split = tree_adjust!(branch,branch.right,x,residual,right_data_indices,false)
    return valid_split
end

function change_decision_rule!(bart_state::BartState,tree::BartTree,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
  branch_nodes = branches(tree)

  number_branches = length(branch_nodes)
  if number_branches ==0
     return false
  end

  i = rand(1:number_branches)
  branch = branch_nodes[i]
  branch_depth = depth(branch)

  branch_indices = train_data_indices(branch)
  #because here we allow leaves with no training data, it's possible that
  #branch_indices is empty
  if length(branch_indices)==0
      return false
  end

  #split_feature = rand(1:size(x,2))
  #feature = x[branch_indices,split_feature]
  split_feature = rand(1:size(x,1))
  feature = vec(x[split_feature,branch_indices])

  split_location = rand(1:length(feature))
  split_value = feature[split_location]

  #calculate the prior tree loglikelihood
  branch_loglikelihood = node_loglikelihood(branch,bart_state)
  branch_prior = log_tree_prior(branch,branch_depth,bartoptions)

  #calculate the proposal tree loglikelihood
  old_feature = branch.feature
  old_value = branch.value

  branch.feature = split_feature
  branch.value = split_value
  valid_split = tree_adjust!(branch,x,residual,branch_indices)
  if !valid_split
    branch.feature = old_feature
    branch.value = old_value
    tree_adjust!(branch,x,residual,branch_indices)
    updated = false
    return updated
  end
  proposal_branch_loglikelihood = node_loglikelihood(branch,bart_state)
  #when the following value equals to Inf, it means that there is a branch that has no training data indices
  proposal_branch_prior = log_tree_prior(branch,branch_depth,bartoptions) #value might be Inf

  #calculate alpha
  if proposal_branch_loglikelihood==Inf
    alpha = 0.0
  else
    alpha = proposal_branch_prior==Inf?0.0:exp(proposal_branch_loglikelihood+proposal_branch_prior-branch_loglikelihood-branch_prior)
  end

  if rand()<alpha
    updated = true
  else
    branch.feature = old_feature
    branch.value = old_value
    tree_adjust!(branch,x,residual,branch_indices)
    updated = false #actually I think here should also set updated = true, since now the leaves values are all zero, thus we need to sample again
  end
  updated
end

function swap_decision_rule!(bart_state::BartState,tree::BartTree,x::Matrix{Float64},residual::Vector{Float64},bartoptions::BartOptions)
  left_or_right(branch::DecisionBranch) = typeof(branch.left)==DecisionBranch?(typeof(branch.right)==DecisionBranch?rand(1:2):1):2
  grand_branch_nodes = grand_branches(tree)
  if length(grand_branch_nodes)==0
    return false
  end

  #pick a parent
  grand_branch = grand_branch_nodes[rand(1:length(grand_branch_nodes))]
  grand_branch_depth = depth(grand_branch)
  grand_branch_indices = train_data_indices(grand_branch)
  type_left = typeof(grand_branch.left)
  type_right = typeof(grand_branch.right)
  swap_both = false

  if type_left==DecisionBranch && type_right==DecisionBranch
    if grand_branch.left.feature==grand_branch.right.feature && grand_branch.left.value==grand_branch.right.value
      swap_both = true
    else
      i = rand(1:2)
      swap_branch = i==1?grand_branch.left:grand_branch.right
    end
  else
      swap_branch = left_or_right(grand_branch)==1?grand_branch.left:grand_branch.right
  end

  grand_branch_loglikelihood = node_loglikelihood(grand_branch,bart_state)
  grand_branch_prior = log_tree_prior(grand_branch,grand_branch_depth,bartoptions)
  old_grand_branch_feature = grand_branch.feature
  old_grand_branch_value = grand_branch.value

  if swap_both==true
     old_branch_feature = grand_branch.left.feature
     old_branch_value = grand_branch.left.value

     grand_branch.feature = old_branch_feature
     grand_branch.value = old_branch_value
     grand_branch.left.feature = grand_branch.right.feature = old_grand_branch_feature
     grand_branch.left.value = grand_branch.right.value = old_grand_branch_value
   else
     old_branch_feature = swap_branch.feature
     old_branch_value = swap_branch.value

     grand_branch.feature = old_branch_feature
     grand_branch.value = old_branch_value
     swap_branch.feature = old_grand_branch_feature
     swap_branch.value = old_grand_branch_value
   end

   valid_split = tree_adjust!(grand_branch,x,residual,grand_branch_indices)
   if !valid_split
     grand_branch.feature = old_grand_branch_feature
     grand_branch.value = old_grand_branch_value
     if swap_both
       grand_branch.left.feature = grand_branch.right.feature = old_branch_feature
       grand_branch.left.value = grand_branch.right.value = old_branch_value
     else
       swap_branch.feature = old_branch_feature
       swap_branch.value = old_branch_value
     end
     tree_adjust!(grand_branch,x,residual,grand_branch_indices)
     updated = false
     return updated
   end
   proposal_grand_branch_loglikelihood = node_loglikelihood(grand_branch,bart_state)
   proposal_grand_branch_prior = log_tree_prior(grand_branch,grand_branch_depth,bartoptions)

   #calculate alpha
   if proposal_grand_branch_loglikelihood==Inf
     alpha = 0.0
   else
     alpha = proposal_grand_branch_prior==Inf?0.0:exp(proposal_grand_branch_loglikelihood+proposal_grand_branch_prior-grand_branch_loglikelihood-grand_branch_prior)
   end

   if rand()<alpha
     updated = true
   else
     grand_branch.feature = old_grand_branch_feature
     grand_branch.value = old_grand_branch_value
     if swap_both
       grand_branch.left.feature = grand_branch.right.feature = old_branch_feature
       grand_branch.left.value = grand_branch.right.value = old_branch_value
     else
       swap_branch.feature = old_branch_feature
       swap_branch.value = old_branch_value
     end
     tree_adjust!(grand_branch,x,residual,grand_branch_indices)
     updated = false #actually I think here should also set updated = true, since now the leaves values are all zero, thus we need to sample again
   end
   updated
end

function update_leaf_values!(tree::BartTree,parameters::BartStateParameters)
  for leaf = leaves(tree) #need to write a function to get all the leaves
    update_leaf_values!(leaf,parameters);
  end
end

function update_leaf_values!(leaf::BartLeaf,parameters::BartStateParameters)
  n = length(leaf.leaf_data_indices);
  post_mean = n*leaf.residual_mean*parameters.sigma_mu^2/(n*parameters.sigma_mu^2+parameters.sigma^2);
  post_sigma = sqrt(parameters.sigma^2*parameters.sigma_mu^2/(n*parameters.sigma_mu^2+parameters.sigma^2));
  leaf.value = post_mean + post_sigma*randn();
end

function update_sigma!(bart_state::BartState,residuals::Vector{Float64})
  alpha = (bart_state.parameters.nu+length(residuals))/2;
  beta = (bart_state.parameters.nu*bart_state.parameters.lambda+sum(residuals.^2))/2
  x = Gamma(alpha,1/beta)  #using package Distributions
  bart_state.parameters.sigma = sqrt(1/rand(x,1)[1])
end

#This two functions perform similar to predict,while significantly faster
function predict_value(bart_tree::BartTree,x::Matrix{Float64})
  leaf_nodes = leaves(bart_tree)
  y_predict = zeros(size(x,2))
  for leaf in leaf_nodes
    y_predict[leaf.leaf_data_indices]=leaf.value
  end
  y_predict
end

function predict_value(bart_state::BartState,x::Matrix{Float64})
   y_predict = zeros(size(x,2))
   for tree in bart_state.trees
      y_predict = y_predict.+predict_value(tree,x)
   end
   y_predict
end

function initialize_bart_state(x::Matrix{Float64},y_normalized::Vector{Float64},bartoptions::BartOptions)
  #number_observations = size(x,1)
  number_observations = size(x,2)
  #initialize a bayesian additive regression tree
  trees = BartTree[]
  #here, assume we estimate y_normalized by its mean, thus estimates of each single tree should be
  #mean(y_normalized)/bartoption.num_trees; then residual of this tree should be
  initial_residuals = y_normalized.-(bartoptions.num_trees-1)*mean(y_normalized)/bartoptions.num_trees

  for i in 1:bartoptions.num_trees
     push!(trees,BartTree(DecisionTree(BartLeaf(0.0,initial_residuals,collect(1:number_observations)))))
  end

  sigma = std(y_normalized);
  sigma_mu = 0.5/(sqrt(bartoptions.num_trees)*bartoptions.k);
  nu = 3.0;
  q = 0.90;
  lambda = sigma^2*quantile(Chisq(nu),0.1)/nu; #using package Distributions
  bartparameters= BartStateParameters(sigma,sigma_mu,nu,lambda);
  bart_state = BartState(trees,bartparameters);
  #given all the tree are one single node, update the leaf values
  for tree = bart_state.trees
      update_leaf_values!(tree,bart_state.parameters)
  end
  yhat = predict_value(bart_state,x)

  for (tree in bart_state.trees)
      y_tree_hat = predict_value(tree,x)

      residuals = y_normalized-(yhat-y_tree_hat)
      tree.tree.root.residual_mean = mean(residuals)
      tree.tree.root.residual_sigma = sqrt(mean((residuals.-mean(residuals)).^2))
      update_leaf_values!(tree,bart_state.parameters)
      #after we update one tree, we update estimate of y
      y_tree_hat_new = predict_value(tree,x)
      yhat+=y_tree_hat_new-y_tree_hat
 end
  bart_state
end

function bart_state_to_additivetree(bart_state::BartState)
  trees = [BartTree(DecisionTree(bart_node_regression_node(tree.tree.root))) for tree in bart_state.trees]
  BartAdditiveTree(trees)
end

function bart_node_regression_node(branch::DecisionBranch)
  DecisionBranch(branch.feature,
                 branch.value,
                 bart_node_regression_node(branch.left),
                 bart_node_regression_node(branch.right))
end

function bart_node_regression_node(leaf::BartLeaf)
  RegressionLeaf(leaf.value)
end

function normalize(response::Vector{Float64},y_min::Float64,y_max::Float64)
   k = 1/(y_max-y_min) #whatif y_max = y_min
   b = -0.5*(y_max+y_min)/(y_max-y_min)
   y_normalized = k*response+b;
   y_normalized
end

function denormalize(predict::Vector{Float64},y_min::Float64,y_max::Float64)
  y_denormalized = predict.*(y_max-y_min)+0.5*(y_max+y_min)
  y_denormalized
end

function StatsBase.fit(x::Vector{Float64},y::Vector{Float64},bartoptions::BartOptions)
    x = reshape(x,length(x),1)
    fit(x,y,bartoptions);
end

function StatsBase.fit(x::Matrix{Float64},y::Vector{Float64},bartoptions::BartOptions)
  println("Running BART with numeric y\n")
  println("Number of trees: ",bartoptions.num_trees)
  println("Prior:")
  println("     k: ",bartoptions.k)
  y_min = minimum(y)
  y_max = maximum(y)
  y_normalized = normalize(y,y_min,y_max)
  number_observations = length(y);
  number_predictors = size(x,2);
  x = x'

  bart_state = initialize_bart_state(x,y_normalized,bartoptions)
  println("     degrees of freedom in sigma prior: ",bart_state.parameters.nu)
  println("     quantile in sigma prior: 0.9")
  println("     power and base for tree prior: ",bartoptions.alpha," ",bartoptions.beta)
  println("     use quantile for rule cut points ",0)
  println("data: ")
  println("     number of training observations: ",number_observations)
  println("     number of explanatory variables: ",number_predictors)

  bart_additive_trees = Array(BartAdditiveTree,0)
  y_hat = predict_value(bart_state,x)

  println("\n")
  println("Running mcmc loop:")

  @time for i in range(1,bartoptions.num_draws+bartoptions.num_burn_in)
           if i % 100 ==0
            println("iteration: ",i," (of ",bartoptions.num_draws+bartoptions.num_burn_in,")")
           end
           updates = 0
           for  j = 1:bartoptions.num_trees
              y_tree_hat = predict_value(bart_state.trees[j],x)

              residual= y_normalized - (y_hat-y_tree_hat)
              updated = update_tree!(bart_state,bart_state.trees[j],x,residual,bartoptions)
              updates+=updated?1:0
              y_tree_hat_new = predict_value(bart_state.trees[j],x)
              y_hat+=y_tree_hat_new-y_tree_hat
              #y_hat += predict(bart_state.trees[j],x)-y_tree_hat
           end
           update_sigma!(bart_state,y_normalized-y_hat)
           #println("there is",updates, "in this iteration")
           if i>bartoptions.num_burn_in
              if i % bartoptions.num_thinning==0
                bart_additive_tree = bart_state_to_additivetree(bart_state)
                #bart_additive_tree = BartAdditiveTree(bart_state.trees)
                push!(bart_additive_trees,bart_additive_tree)
              end
           end
        end
   Bart(bart_additive_trees,y_min,y_max,bartoptions)
end

function StatsBase.predict(tree::BartTree,sample::Vector{Float64}) #import package StatsBase
  node = tree.tree.root
  while typeof(node)==DecisionBranch
    if sample[node.feature]<=node.value
      node = node.left
    else
      node = node.right
    end
  end
  node.value
end

function StatsBase.predict(bart_state::BartState,sample::Vector{Float64})
  Float64(sum([predict(tree,sample) for tree in bart_state.trees]))
end

function StatsBase.predict(bart_state::BartState,x::Matrix{Float64})
  [Float64(predict(bart_state,vec(x[:,i]))) for i in range(1,size(x,2))]
end

function StatsBase.predict(tree::BartTree,x::Matrix{Float64})
   [Float64(predict(tree,vec(x[:,i]))) for i in range(1,size(x,2))]
end

function StatsBase.predict(bart_additive_tree::BartAdditiveTree,x::Matrix{Float64})
  y_predict = zeros(size(x,2))
  for i in range(1,size(x,2))
    y_predict[i] = Float64(sum([predict(tree,vec(x[:,i])) for tree in bart_additive_tree.trees]))
  end
  y_predict
end

function StatsBase.predict(bart::Bart,x::Vector{Float64},confidence_interval=false::Bool)
   x = reshape(x,length(x),1)
   StatsBase.predict(bart,x,confidence_interval)
end

function StatsBase.predict(bart::Bart,x::Matrix{Float64},confidence_interval=false::Bool)
   x = x'
   y_predict = zeros(size(x,2),length(bart.bart_additive_trees))
   if confidence_interval
      y_CI = zeros(size(x,2),2)
      for i in range(1,length(bart.bart_additive_trees))
      y_predict[:,i] = predict(bart.bart_additive_trees[i],x)
      y_predict[:,i] = denormalize(y_predict[:,i],bart.y_min,bart.y_max)
      end

      for i in 1:size(x,2)
        y_CI[i,:] = quantile(vec(y_predict[i,:]),[0.025,0.975])
      end
      #point estimate
      y_hat = vec(mean(y_predict,2))
      y_hat,y_CI
   else
      for i in range(1,length(bart.bart_additive_trees))
        y_predict[:,i] = predict(bart.bart_additive_trees[i],x)
        y_predict[:,i] = denormalize(y_predict[:,i],bart.y_min,bart.y_max)
      end
      #point estimate
      y_hat = vec(mean(y_predict,2))
      y_hat
   end
end

function model_selection(x::Vector{Float64},y::Vector{Float64},bartoptions::BartOptions,num_trees=5::Int)
   x = reshape(x,length(x),1)
   model_selection(x,y,bartoptions,num_trees);
end

function model_selection(x::Matrix{Float64},y::Vector{Float64},bartoptions::BartOptions,num_trees=5::Int)
  println("Running BART with numeric y\n")
  bartoptions.num_trees = num_trees;
  println("Number of trees: ",bartoptions.num_trees)
  println("Prior:")
  println("     k: ",bartoptions.k)
  y_min = minimum(y)
  y_max = maximum(y)
  y_normalized = normalize(y,y_min,y_max)
  number_observations = length(y);
  number_predictors = size(x,2);
  x = x'

  bart_state = initialize_bart_state(x,y_normalized,bartoptions)
  println("     degrees of freedom in sigma prior: ",bart_state.parameters.nu)
  println("     quantile in sigma prior: 0.9")
  println("     power and base for tree prior: ",bartoptions.alpha," ",bartoptions.beta)
  println("     use quantile for rule cut points ",0)
  println("data: ")
  println("     number of training observations: ",number_observations)
  println("     number of explanatory variables: ",number_predictors)

  bart_additive_trees = Array(BartAdditiveTree,0)
  y_hat = predict_value(bart_state,x)
  println("\n")
  count = zeros(number_predictors);
  println("Running mcmc loop:")

  @time for i in range(1,bartoptions.num_draws+bartoptions.num_burn_in)
           if i % 100 ==0
            println("iteration: ",i," (of ",bartoptions.num_draws+bartoptions.num_burn_in,")")
           end
           updates = 0
           for  j = 1:bartoptions.num_trees
              y_tree_hat = predict_value(bart_state.trees[j],x)
              residual= y_normalized - (y_hat-y_tree_hat)
              updated = update_tree!(bart_state,bart_state.trees[j],x,residual,bartoptions)
              updates+=updated?1:0
              y_hat += predict_value(bart_state.trees[j],x)-y_tree_hat
           end
           update_sigma!(bart_state,y_normalized-y_hat)
           #println("there is",updates, "in this iteration")
           if i>bartoptions.num_burn_in
              if i % bartoptions.num_thinning==0 && i>1000
                 for j in 1:bartoptions.num_trees
                    branch_nodes = branches(bart_state.trees[j])
                    for branch in branch_nodes
                        count[branch.feature]+=1
                    end
                 end
              end
           end
        end
   count
end

max_stars_repo_name: mathcg/BayesTree

max_issues_repo_name: mathcg/BayesTree

max_forks_repo_name: mathcg/BayesTree


text:
using Parquet
using Test
using Dates

function test_load(file::String)
    p = Parquet.File(file)
    @debug("load file", file)
    @test isa(p.meta, Parquet.FileMetaData)

    rgs = rowgroups(p)
    @test length(rgs) > 0

    cnames = colnames(p)
    @test length(cnames) > 0
    @debug("data columns", file, cnames)

    for rg in rgs
        ccs = columns(p, rg)
        @debug("reading row group", file, ncolumnchunks=length(ccs))

        for cc in ccs
            npages = 0
            ccp = Parquet.ColumnChunkPages(p, cc)
            result = iterate(ccp)
            npages = 0
            ncompressedbytes = 0
            nuncompressedbytes = 0
            while result !== nothing
                page,nextpos = result
                result = iterate(ccp, nextpos)
                npages += 1
                ncompressedbytes += Parquet.page_size(page.hdr)
                nuncompressedbytes += page.hdr.uncompressed_page_size
            end
            @test npages > 0
            @test ncompressedbytes > 0
            @test nuncompressedbytes >= 0
            @debug("read column chunk", file, npages, ncompressedbytes, nuncompressedbytes)
        end
    end

    @debug("done loading file", file, showbuffer=String(sb))

    iob = IOBuffer()
    show(iob, p)
    sb = take!(iob)
    @test !isempty(sb)
    @debug("parquet file show", file, showbuffer=String(sb))

    iob = IOBuffer()
    show(iob, p.meta)
    sb = take!(iob)
    @test !isempty(sb)
    @debug("parquet file metadata show", file, showbuffer=String(sb))

    show(iob, schema(p))
    sb = take!(iob)
    @test !isempty(sb)
    @debug("schema show", file, showbuffer=String(sb))

    nothing
end

function test_decode(file)
    p = Parquet.File(file)
    @debug("load file", file)
    @test isa(p.meta, Parquet.FileMetaData)

    rgs = rowgroups(p)
    @test length(rgs) > 0

    for rg in rgs
        ccs = columns(p, rg)
        @debug("reading row group", file, ncolumnchunks=length(ccs))

        for cc in ccs
            jtype = Parquet.elemtype(Parquet.elem(schema(p), colname(p,cc)))
            npages = 0
            ccpv = Parquet.ColumnChunkPageValues(p, cc, jtype)
            result = iterate(ccpv)
            valcount = repncount = defncount = 0
            while result !== nothing
                resultdata,nextpos = result
                valcount += resultdata.value.offset
                repncount += resultdata.repn_level.offset
                defncount += resultdata.defn_level.offset
                result = iterate(ccpv, nextpos)
                npages += 1
            end

            @debug("read", file, npages, valcount, defncount, repncount)
        end
    end
    nothing
end

function test_decode_all_pages()
    @testset "decode parquet-compatibility test files" begin
        testfolder = joinpath(@__DIR__, "parquet-compatibility")
        for encformat in ("SNAPPY", "GZIP", "NONE")
            for fname in ("nation", "customer")
                testfile = joinpath(testfolder, "parquet-testdata", "impala", "1.1.1-$encformat", "$fname.impala.parquet")
                test_decode(testfile)
            end
        end
    end

    @testset "decode julia-parquet-compatibility test files" begin
        testfolder = joinpath(@__DIR__, "julia-parquet-compatibility")
        for encformat in ("ZSTD", "SNAPPY", "GZIP", "NONE")
            for fname in ("nation", "customer")
                testfile = joinpath(testfolder, "Parquet_Files", "$(encformat)_pandas_pyarrow_$(fname).parquet")
                test_decode(testfile)
            end
        end
    end

    @testset "decode missingvalues test file" begin
        test_decode(joinpath(@__DIR__, "missingvalues", "synthetic_data.parquet"))
    end
end

function test_load_all_pages()
    @testset "load parquet-compatibility test files" begin
        testfolder = joinpath(@__DIR__, "parquet-compatibility")
        for encformat in ("SNAPPY", "GZIP", "NONE")
            for fname in ("nation", "customer")
                testfile = joinpath(testfolder, "parquet-testdata", "impala", "1.1.1-$encformat", "$fname.impala.parquet")
                test_load(testfile)
            end
        end
    end

    @testset "load julia-parquet-compatibility test files" begin
        testfolder = joinpath(@__DIR__, "julia-parquet-compatibility")
        for encformat in ("ZSTD", "SNAPPY", "GZIP", "NONE")
            for fname in ("nation", "customer")
                testfile = joinpath(testfolder, "Parquet_Files", "$(encformat)_pandas_pyarrow_$(fname).parquet")
                test_load(testfile)
            end
        end
    end
end

function test_load_boolean_and_ts()
    @testset "load booleans and timestamps" begin
        @debug("load booleans and timestamps")
        p = Parquet.File(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"))

        rg = rowgroups(p)
        @test length(rg) == 1
        cc = columns(p, 1)
        @test length(cc) == 11
        cnames = colnames(p)
        @test length(cnames) == length(cc)
        @test cnames[2] == ["bool_col"]

        rc = RecordCursor(p; rows=1:2, colnames=colnames(p))
        @test length(rc) == 2
        @test eltype(rc) == NamedTuple{(:id, :bool_col, :tinyint_col, :smallint_col, :int_col, :bigint_col, :float_col, :double_col, :date_string_col, :string_col, :timestamp_col),Tuple{Union{Missing, Int32},Union{Missing, Bool},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int32},Union{Missing, Int64},Union{Missing, Float32},Union{Missing, Float64},Union{Missing, Array{UInt8,1}},Union{Missing, Array{UInt8,1}},Union{Missing, DateTime}}}

        values = collect(rc)
        @test [v.bool_col for v in values] == [true,false]
        @test [v.timestamp_col for v in values] == [DateTime("2009-04-01T12:00:00"), DateTime("2009-04-01T12:01:00")]

        cc = BatchedColumnsCursor(p)
        values, _state = iterate(cc)
        @test values.timestamp_col == [DateTime("2009-04-01T12:00:00"), DateTime("2009-04-01T12:01:00")]

        p = Parquet.File(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"); map_logical_types=Dict(["date_string_col"]=>(String,logical_string)))
        rc = RecordCursor(p; rows=1:2, colnames=colnames(p))
        values = collect(rc)
        @test [v.date_string_col for v in values] == ["04/01/09", "04/01/09"]

        cc = BatchedColumnsCursor(p)
        values, _state = iterate(cc)
        @test values.date_string_col == ["04/01/09", "04/01/09"]

        p = Parquet.File(joinpath(@__DIR__, "booltest", "alltypes_plain.snappy.parquet"); map_logical_types=Dict(["timestamp_col"]=>(DateTime,(v)->logical_timestamp(v; offset=Dates.Second(30)))))
        rc = RecordCursor(p; rows=1:2, colnames=colnames(p))
        values = collect(rc)
        @test [v.timestamp_col for v in values] == [DateTime("2009-04-01T12:00:30"), DateTime("2009-04-01T12:01:30")]

        cc = BatchedColumnsCursor(p)
        values, _state = iterate(cc)
        @test values.timestamp_col == [DateTime("2009-04-01T12:00:30"), DateTime("2009-04-01T12:01:30")]
        #dlm,headers=readdlm("booltest/alltypes.csv", ','; header=true)
        #@test [v.bool_col for v in values] == dlm[:,2]  # skipping for now as this needs additional dependency on DelimitedFiles
    end
end

function test_load_nested()
    @testset "load nested columns" begin
        @debug("load nested columns")
        p = Parquet.File(joinpath(@__DIR__, "nested", "nested1.parquet"))

        @test nrows(p) == 100
        @test ncols(p) == 5

        rc = RecordCursor(p)
        @test length(rc) == 100
        @test eltype(rc) == NamedTuple{(:_adobe_corpnew,),Tuple{NamedTuple{(:id, :vocab, :frequency, :max_len, :reduced_max_len),Tuple{Union{Missing, Int32},Union{Missing, String},Union{Missing, Int32},Union{Missing, Float64},Union{Missing, Int32}}}}}

        values = Any[]
        for rec in rc
            push!(values, rec)
        end

        v = values[1]._adobe_corpnew
        @test v.frequency == 3
        @test v.id == 1375
        @test v.max_len == 64192.0
        @test v.reduced_max_len == 64
        @test v.vocab == "10385911_a"

        v = values[100]._adobe_corpnew
        @test v.frequency == 61322
        @test v.id == 724
        @test v.max_len == 64192.0
        @test v.reduced_max_len == 64
        @test v.vocab == "12400277_a"

        p = Parquet.File(joinpath(@__DIR__, "nested", "nested.parq"))

        @test nrows(p) == 10
        @test ncols(p) == 1

        rc = RecordCursor(p)
        @test length(rc) == 10
        @test eltype(rc) == NamedTuple{(:nest,),Tuple{Union{Missing, NamedTuple{(:thing,),Tuple{Union{Missing, NamedTuple{(:list,),Tuple{Array{NamedTuple{(:element,),Tuple{Union{Missing, String}}},1}}}}}}}}}

        values = collect(rc)
        v = first(values)
        @test length(v.nest.thing.list) == 2
        @test v.nest.thing.list[1].element == "hi"
        v = last(values)
        @test length(v.nest.thing.list) == 2
        @test v.nest.thing.list[1].element == "world"
    end
end

function test_load_multiple_rowgroups()
    @testset "testing multiple rowgroups..." begin
        @debug("load multiple rowgroups")
        p = Parquet.File(joinpath(@__DIR__, "rowgroups", "multiple_rowgroups.parquet"))

        @test nrows(p) == 100
        @test ncols(p) == 12

        rc = RecordCursor(p)
        @test length(rc) == 100
        vals = collect(rc)
        @test length(vals) == 100
        @test vals[1].int64 == vals[51].int64
        @test vals[1].int32 == vals[51].int32

        cc = BatchedColumnsCursor(p)
        @test length(cc) == 2
        colvals = collect(cc)
        @test length(colvals) == 2
        @test length(colvals[1].int32) == 50
    end
end

function test_load_file()
    @testset "load a file" begin
        df = read_parquet(joinpath(@__DIR__, "rowgroups", "multiple_rowgroups.parquet"))

        #all columns must be 100 rows long
        @test all([length(val)==100 for (_, val) in df])

        # 12 columns
        @test length(df) == 12
    end
end

function test_load_at_offset()
    @testset "load file at offset" begin
        testfolder = joinpath(@__DIR__, "parquet-compatibility")
        testfile = joinpath(testfolder, "parquet-testdata", "impala", "1.1.1-NONE", "customer.impala.parquet")
        parquet_file = Parquet.File(testfile)

        vals_20000_40000 = first(collect(Parquet.BatchedColumnsCursor(parquet_file; rows=20000:40000))).c_custkey
        vals_1_40000 = first(collect(Parquet.BatchedColumnsCursor(parquet_file; rows=1:40000))).c_custkey
        vals_2_40001 = first(collect(Parquet.BatchedColumnsCursor(parquet_file; rows=2:40001))).c_custkey

        @test vals_20000_40000 == vals_1_40000[20000:40000]
        @test vals_20000_40000 != vals_1_40000[1:20000]
        @test vals_2_40001[1:10000] == vals_1_40000[2:10001]
        @test vals_2_40001[1:10000] != vals_1_40000[1:10000]
    end
end

@testset "load files" begin
    test_load_all_pages()
    test_decode_all_pages()
    test_load_boolean_and_ts()
    test_load_nested()
    test_load_multiple_rowgroups()
    test_load_file()
    test_load_at_offset()
end


text:
# ------------------------------------------------------------------------------------------
# # Representing data in a computer
#
# The core of data science and machine learning is **data**: we are interested in extracting
# knowledge from data.
#
# But how exactly do computers represent data? Let's find out exactly what an "artificial
# intelligence" has at its disposal to learn from.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Data is represented as arrays
#
# Let's take a look at some fruit. Using the `Images.jl` library, we can load in some
# images:
# ------------------------------------------------------------------------------------------

using Images

apple = load("data/10_100.jpg")

banana = load("data/104_100.jpg")

# ------------------------------------------------------------------------------------------
# Here we have images of apples and bananas. We would eventually like to build a program
# that can automatically distinguish between the two. However, the computer doesn't "see" an
# apple or a banana; instead, it just sees numbers.
#
# An image is encoded in something called an **array**, which is like a container that has
# boxes or slots for individual pieces of data:
#
# ![attachment:array_cartoon.png](data/array_cartoon.png)
#
# An array is a bunch of numbers in connected boxes; the figure above shows a 1-dimensional
# array. Our images are instead 2-dimensional arrays, or matrices, of numbers, arranged
# something like this:
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![attachment:array2d.png](data/array2d.png)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# For example, `apple` is an image, consisting of a 100x100 array of numbers:
# ------------------------------------------------------------------------------------------

typeof(apple)

typeof(a)

a = [ 1 2 3;4 5 6]

size(a)

# ------------------------------------------------------------------------------------------
# We can grab the datum stored in the box at row `i` and column `j` by *indexing* using
# square brackets: `[i, j]`. For example, let's get the pixel (piece of the image) in box
# $(40, 60)$, i.e. in the 40th row and 60th column of the image:
# ------------------------------------------------------------------------------------------

apple

dump(typeof(apple[40, 60]))

apple[18:20,29:31]

# ------------------------------------------------------------------------------------------
# We see that Julia displays a coloured box! Julia, via the `Colors.jl` package, is clever
# enough to display colours in a way that is useful to us humans!
#
# So, in fact, an image is a 2D array, in which each element of the array is an object (a
# collection of numbers) describing a coloured pixel.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Colors as numbers
#
# How, then, are these colors actually stored? Computers store colors in RGB format, that is
# they store a value between 0 and 1 for each of three "channels": red, green, and blue.
# Here, 0 means none of that color and 1 means the brightest form of that color. The overall
# color is a combination of those three colors.
#
# For example, we can pull out the `red` value using the function `red` applied to the
# color. Since internally the actual value is stored in a special format, we choose to
# convert it to a standard floating-point number using the `Float64` function:
# ------------------------------------------------------------------------------------------

Float64(red(apple[40, 60]))

[ mean(float.(c.(img))) for c = [red,green,blue], img = [apple,banana] ]

using Plots

histogram(float.(green.(apple[:])),color="red",label="apple", normalize=true, nbins=25)
histogram!(float.(green.(banana[:])),color="yellow",label="banana",normalize=true, nbins=25)

apple

float(red(banana[50,20]))

banana[50,20]

pixel = apple[40, 60]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

# ------------------------------------------------------------------------------------------
# Since the red value is high while the others are low, this means that at pixel `(40, 60)`,
# the picture of the apple is very red. If we do the same at one of the corners of the
# picture, we get the following:
# ------------------------------------------------------------------------------------------

pixel = apple[1, 1]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

apple

# ------------------------------------------------------------------------------------------
# We see that every color is bright, which corresponds to white.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Working on an image as a whole
#
# In Julia, to apply a function to the whole of an array, we place a `.` between the
# function name and the left parenthesis (`(`), so the following gives us the `red` value of
# every pixel in the image:
# ------------------------------------------------------------------------------------------

redpartofapple = Float64.(red.(apple))
mean(redpartofapple)

using Plots

gr()

histogram(redpartofapple[:],color=:red,label="redness in the apple")

# ------------------------------------------------------------------------------------------
# Note that we get a 2D array (matrix) back.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia's [mathematical standard
# library](https://docs.julialang.org/en/stable/stdlib/math/#Mathematics-1) has many
# mathematical functions built in. One of them is the `mean` function, which computes the
# average value. If we apply this to our apple:
# ------------------------------------------------------------------------------------------

mean(Float64.(red.(apple)))

# ------------------------------------------------------------------------------------------
# we see that the value indicates that the average amount of red in the image is a value
# between the amount of red in the apple and the amount of red in the white background.
#
# *Somehow we need to teach a computer to use this information about a picture to recognize
# that there's an apple there!*
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## A quick riddle
#
# Here's a quick riddle. Let's check the average value of red in the image of the banana.
# ------------------------------------------------------------------------------------------

mean(Float64.(red.(banana)))

# ------------------------------------------------------------------------------------------
# Oh no, that's more red than our apple? This isn't a mistake and is actually true! Before
# you move onto the next exercise, examine the images of the apple and the banana very
# carefully and see if you can explain why this is expected.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# What is the average value of blue in the banana?
#
# (To open a new box use <ESC>+b (b is for "below", what do you think a does?))
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Does the banana have more blue or more green?
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/01. ML - Representing data in a computer.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials