text:
module CRT

using JFVM
using PyPlot

export CRT_SNIA, CRT_SIA, CRT_GIA_eps, CRT_GIA_kin_eq, CRT_GIA_kin_fb, plot_concentrations

#Define functions
function epsilon(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
    # A + B <-> C
    return 1/(2*nu_A*nu_B*K_eq)*(nu_A*C_B*K_eq+nu_B*C_A*K_eq+nu_C-sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2))
end

function epsilon_A(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return 1/(2*nu_A*nu_B*K_eq)*(nu_B*K_eq-(C_A*K_eq^2*nu_B^2-C_B*K_eq^2*nu_A*nu_B+K_eq*nu_B*nu_C)/sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2))
end

function epsilon_AA(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return 1/(2*nu_A*nu_B*K_eq)*(-(K_eq^2*nu_B^2)/sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)+((C_A*K_eq^2*nu_B^2-C_B*K_eq^2*nu_A*nu_B+K_eq*nu_B*nu_C)*(2*C_A*K_eq^2*nu_B^2-2*C_B*K_eq^2*nu_A*nu_B+2*K_eq*nu_B*nu_C))/(2*(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)^(3/2)))
end

function epsilon_B(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return 1/(2*nu_A*nu_B*K_eq)*(nu_A*K_eq-(C_B*K_eq^2*nu_A^2-C_A*K_eq^2*nu_A*nu_B+K_eq*nu_A*nu_C)/sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2))
end

function epsilon_BB(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return 1/(2*nu_A*nu_B*K_eq)*(-(K_eq^2*nu_A^2)/sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)+((-C_A*K_eq^2*nu_A*nu_B+C_B*K_eq^2*nu_A+K_eq*nu_A*nu_C)*(-2*C_A*K_eq^2*nu_A*nu_B+2*C_B*K_eq^2*nu_A^2+2*K_eq*nu_A*nu_C))/(2*(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)^(3/2)))
end

function epsilon_C(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return -1/sqrt(nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)
end

function epsilon_CC(C_A, C_B, C_C, nu_A, nu_B, nu_C, K_eq)
  return (2*K_eq*nu_A*nu_B)/((nu_A^2*C_B^2*K_eq^2-2*nu_A*nu_B*C_A*C_B*K_eq^2+nu_B^2*C_A^2*K_eq^2+4*nu_A*nu_B*C_C*K_eq+2*nu_A*nu_C*C_B*K_eq+2*nu_B*nu_C*C_A*K_eq+nu_C^2)^(3/2))
end

include("CRT_SNIA.jl")
include("CRT_SIA.jl")
include("CRT_GIA_eps.jl")
include("CRT_GIA_kin_eq.jl")
include("CRT_GIA_kin_fb.jl")
include("plot_concentrations.jl")

end

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
using Symbolics
using ForwardDiff
using LinearAlgebra
using Plots
using ForwardDiff

struct Segment{T}
    domain::Tuple{T,T}
    func
    f
    df
    ddf
    r
    rx
    rp
    δp
    dir_norm
    δ_dir_norm
end

function Segment(domain, func)
    @variables x p[1:2]

    # surface
    f = func(x)
    df = Symbolics.derivative(f, x)
    ddf = Symbolics.derivative(df, x)

    # distance (squared)
    r = x - p[1] + (f - p[2]) * df
    rx = 1.0 + df * df + (f - p[2]) * ddf
    rp = [-1.0 -df]
    δp = [-rp ./ rx; -df * rp ./ rx]

    # fast methods
    f_func = eval(Symbolics.build_function(f, x))
    df_func = eval(Symbolics.build_function(df, x))
    ddf_func = eval(Symbolics.build_function(ddf, x))

    r_func = eval(Symbolics.build_function(r, x, p))
    rx_func = eval(Symbolics.build_function(rx, x, p))
    fp_func = eval(Symbolics.build_function(rp, x, p)[1])
    δp_func = eval(Symbolics.build_function(δp, x, p)[1])

    # surface norm
    f2 = func(p[1])
    dir = [-Symbolics.derivative(f2, p[1]); 1.0]
    dir_norm = dir ./ norm(dir)
    dir_norm_func = eval(Symbolics.build_function(dir_norm, p)[1])
    δ_dir_norm = Symbolics.jacobian(dir_norm, p)
    δ_dir_norm_func = eval(Symbolics.build_function(δ_dir_norm, p)[1])

    Segment(domain, func,
        f_func, df_func, ddf_func,
        r_func, rx_func, fp_func,
        δp_func,
        dir_norm_func, δ_dir_norm_func)
end
# seg1 = Segment((-10.0, 10.0), x -> 0.5 * x^2.0)

seg1 = Segment((-2.0, 0.5), x -> 0.0 * x)
seg2 = Segment((0.5, 1.5), x -> x - 0.5)
# seg3 = Segment((1.5, 2.5), x -> 1.0)
# seg4 = Segment((2.5, 2.51), x -> 100.0 * x - 249.0)
# seg5 = Segment((2.51, 3.51), x -> 2.0)
# seg6 = Segment((3.51, 3.52), x -> -100.0 * x + 353.0)
# seg7 = Segment((3.52, 4.52), x -> 1.0)

surf = [seg1, seg2]#, seg3]#, seg4]#, seg5, seg6, seg7]

get_point(s::Segment, a) = [a; s.f(a)]

p = [0.0; 1.0]

function minimum_distance_point(s::Segment, p)
    # initialization
    if p[1] < s.domain[1]
        a = s.domain[1]
    elseif p[1] > s.domain[2]
        a = s.domain[2]
    else
        a = p[1]
    end

    # Newton
    for i = 1:100
        abs(s.r(a, p)) < 1.0e-8 && break
        # step
        Δ = -s.rx(a, p) \ s.r(a, p)
        α = 1.0
        a_cand = a + α * Δ
        # linesearch
        iter = 0
        while abs(s.r(a_cand, p)) >= abs(s.r(a, p)) && iter < 100
            α *= 0.5
            a_cand = a + α * Δ
            iter += 1
        end
        a = a_cand
    end

    # solution
    if a < s.domain[1]
        a = s.domain[1]
    elseif a > s.domain[2]
        a = s.domain[2]
    end

    # point and sensitivity
    return get_point(s, a), s.δp(a, p)
end

get_surf_dir(s::Segment, a) = [-s.df(a[1]); 1.0]

function surface_normal(s::Segment, p)
    dir = get_surf_dir(s, p)
    dir_norm = dir ./ norm(dir)
    # dir_norm = s.dir_norm(p)

    θ = -0.5 * π
    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]
    tan_norm = R * dir_norm

    return dir_norm, tan_norm
end

function rotation(ns)
	# world-frame normal
	nw = [0.0, 1.0]

	angle = atan(nw[2], nw[1]) - atan(ns[2], ns[1])

	# Rw->s
	[cos(angle) -sin(angle); sin(angle) cos(angle)]
end

function check_domain(surf::Vector{Segment{T}}, p; r=0.0) where T
    status = Bool[]
    for s in surf
        if p[1] >= s.domain[1] - r && p[1] <= s.domain[2] + r
            push!(status, true)
        else
            push!(status, false)
        end
    end
    return status
end

function contact_point(s::Segment, p)
    pmin, δp = minimum_distance_point(s, p)
    d = norm(p - pmin)
    dir_norm, tan_norm = surface_normal(s, pmin)
    rot = rotation(dir_norm)
    δ_dir_norm = s.δ_dir_norm(pmin)
    return pmin, d, dir_norm, tan_norm, rot, δp, δ_dir_norm
end

function contact_point(surf::Vector{Segment{T}}, p; r=0.0) where T
    domain = check_domain(surf, p, r=r)
    p_eval = []
    for (i, status) in enumerate(domain)
        if status
            push!(p_eval, contact_point(surf[i], p))
        end
    end
    min_idx = findmin([pt[2] for pt in p_eval])[2]
    return p_eval[min_idx]
end

function plot_segment(s::Segment; N=100)
    x = range(s.domain[1], stop=s.domain[2], length=N)
    y = s.f.(x)
    plot!(x, y, width=2.0, color=:black, label="")
end

function visualize_segment!(vis, s::Segment;
    idx = 1, color=RGBA(0.0, 0.0, 0.0, 1.0), linewidth=50, num_points=100)

    line_mat = LineBasicMaterial(color=color, linewidth=linewidth)
    domain = range(s.domain[1], stop=s.domain[2], length=num_points)
    points = [Point(x, 0.0, s.func(x)) for x in domain]
    setobject!(vis["traj_$idx"], MeshCat.Line(points, line_mat))
end

function visualize!(vis, surf::Vector{Segment{T}};
    colors=[RGBA(0.0, 0.0, 0.0, 1.0) for i = 1:length(surf)], linewidth=50, num_points=100) where T

    for (i, s) in enumerate(surf)
        visualize_segment!(vis, s, idx = i, color=colors[i], linewidth=linewidth, num_points=num_points)
    end
end

plt = plot(xlims=(-0.5, 4.0), ylims=(-1.0, 3.0), aspect_ratio=:equal)
for s in surf
    plt = plot_segment(s)
end

# p = [2.0; 1.5]
# p = [0.0, 1.0]
# p = [0.0, 1.4]
p = [0.5, 1.0]
p = [1.5, 1.25]

plt = scatter!([p[1]], [p[2]], color=:red, label="")

p_eval = [contact_point(s, p) for s in surf]
for x in p_eval
    plt = plot!([p[1]; x[1][1]], [p[2]; x[1][2]], width=1.5, color=:cyan, label="")
end
display(plt)

for (i, pc) in enumerate(p_eval)
    @show x = pc[1][1]
    @show y = pc[1][2]

    dir_norm = pc[3]
    tan_norm = pc[4]

    scale = 0.25
    plt = plot!([x; x + scale * dir_norm[1]], [y; y + scale * dir_norm[2]], width=1.5, color=:magenta, label="")
    plt = plot!([x; x + scale * tan_norm[1]], [y; y + scale * tan_norm[2]], width=1.5, color=:orange, label="")

end

display(plt)

function generate_polynomial(seg1::Segment{T}, seg2::Segment{T}; ϵ=0.01) where {T}
    # segments
    seg1.domain[2]
    seg2.domain[1]

    # offsets
    x⁺ = seg1.domain[2] + ϵ
    x⁻ = seg1.domain[2] - ϵ

    # end points and their derivatives
    p1x = seg1.func(x⁻)
    dp1x = ForwardDiff.derivative(seg1.func, x⁻)
    p2x = seg2.func(x⁺)
    dp2x = ForwardDiff.derivative(seg2.func, x⁺)

    # fit polynomial parameters
    b = [p1x; dp1x; p2x; dp2x]
    A = [1.0 x⁻ x⁻^2.0 x⁻^3.0;
        0.0 1.0 2.0 * x⁻ 3.0 * x⁻^2.0;
        1.0 x⁺ x⁺^2.0 x⁺^3.0;
        0.0 1.0 2.0 * x⁺ 3.0 * x⁺^2.0]
    param = A \ b

    # polynomial
    function poly(y)
        param[1] + param[2] * y + param[3] * y^2.0 + param[4] * y^3.0
    end
    @assert abs(poly(x⁻) - seg1.func(x⁻)) < 1.0e-8
    @assert abs(poly(x⁺) - seg2.func(x⁺)) < 1.0e-8

    # polynomial segment
    Segment((x⁻, x⁺), poly)
end

# seg_poly = generate_polynomial(seg1, seg2)

function polynomial_blend(seg1::Segment{T}, seg2::Segment{T}; ϵ=0.01) where {T}
    # polynominal segment
    seg_poly = generate_polynomial(seg1, seg2, ϵ=ϵ)

    # update segments
    seg1_update = Segment((seg1.domain[1], seg_poly.domain[1]), seg1.func)
    seg2_update = Segment((seg_poly.domain[2], seg2.domain[2]), seg2.func)

    return seg1_update, seg_poly, seg2_update
end

seg1_update, seg_poly, seg2_update = polynomial_blend(seg1, seg2)

function polynomial_blend(surf::Vector{Segment{T}}; ϵ=0.01) where {T}
    @assert length(surf) > 1
    seg1_update, seg_poly, seg2_update = polynomial_blend(surf[1], surf[2], ϵ=ϵ)
    surf_update = [seg1_update, seg_poly, seg2_update]
    if length(surf) > 2
        for i = 2:length(surf)-1
            seg1_update, seg_poly, seg2_update = polynomial_blend(surf_update[end], surf[i+1], ϵ=ϵ)
            pop!(surf_update)
            push!(surf_update, seg1_update)
            push!(surf_update, seg_poly)
            push!(surf_update, seg2_update)
        end
    end
    return surf_update
end

surf_update = polynomial_blend(surf, ϵ=0.1)

# plt = plot(xlims=(1.45, 1.55), ylims=(0.95, 1.05), aspect_ratio=:equal)
plt = plot(xlims=(-0.5, 4.0), ylims=(-1.0, 2.5), aspect_ratio=:equal)
# plt = plot(xlims=(3.45, 3.55), ylims=(0.95, 1.05), aspect_ratio=:equal)
plt = plot(xlims=(0.0, 1.0), ylims=(-0.1, 0.25), aspect_ratio=:equal)

for s in surf_update
    plt = plot_segment(s)
end
display(plt)

p = [1.4; 1.25]
plt = scatter!([p[1]], [p[2]], color=:red, label="")

p_eval = [contact_point(s, p) for s in surf_update]
for x in p_eval
    plt = plot!([p[1]; x[1][1]], [p[2]; x[1][2]], width=1.5, color=:cyan, label="")
end
display(plt)

for (i, pc) in enumerate(p_eval)
    @show x = pc[1][1]
    @show y = pc[1][2]

    dir_norm = pc[3]
    tan_norm = pc[4]

    scale = 0.25
    plt = plot!([x; x + scale * dir_norm[1]], [y; y + scale * dir_norm[2]], width=2.0, color=:magenta, label="")
    plt = plot!([x; x + scale * tan_norm[1]], [y; y + scale * tan_norm[2]], width=2.0, color=:orange, label="")
end
display(plt)

# # test rotation
# p_eval[3][1]
# rotation(p_eval[6][3])' * [1.0; 0.0] # map vector in segment to world frame

ext: jl


text:
using Distributions, DataFrames, Plots

mutable struct DistributionalQAgent
    αp::Array{Real, 1}
    αn::Array{Real, 1}
    q::Array{Real, 1}
    function DistributionalQAgent(αp::Array{Real, 1}, αn::Array{Real, 1}, n::Real)
        q = Array{Real, 1}(zeros(n))
        new(αp, αn, q)
    end
end

function generate_lrs(min::Real, max::Real, n::Real)
    quantile(Uniform(min, max), 0.:(1 / (n - 1)):1.)
end

function update(agent::DistributionalQAgent, reward::Real)
    δ = ifelse.(reward .- agent.q .<= 0, -1., 1.)
    agent.q .+= δ .* ifelse.(δ .<= 0, agent.αn, agent.αp)
end

struct Environment
    p::Array{Real, 1}
    value::Array{Real, 1}
end

function step(env::Environment)
    i = rand(Categorical(env.p))
    env.value[i]
end

function step(env::Environment, agent::DistributionalQAgent)
    reward = step(env)
    update(agent, reward)
    DataFrame(agent.q', :auto)
end

N = 50
T = 1000
WARMUP = 100

αp = Array{Real, 1}(generate_lrs(0.01, 0.1, 50))
αn = reverse(αp)
agent = DistributionalQAgent(αp, αn, N)
env = Environment([0.5, 0.5], [1., 2.])

result = Array{Real, 2}(reduce(vcat, map(_ -> step(env, agent), 1:T)))
result_without_warmup = reshape(result[WARMUP+1:end, :], ((T - WARMUP) * N, 1))
p1 = plot(result, label="", xlabel="Trial", ylabel="Q-value")
p2 = histogram(result_without_warmup, label="", xlabel="Q-value", ylabel="counts")
plot(p1, p2)


text:
using TropicalYao.Reversible.LogLikeNumbers
using Test

@testset "base" begin
    include("base.jl")
end

@testset "blas" begin
    include("blas.jl")
end

@testset "einsum" begin
    include("einsum.jl")
end


text:
using ProgressMeter
using StatsBase
using LinearAlgebra
using Printf
using DelimitedFiles

include("utils.jl")

seconds_from_days(days) = days*24*60*60

#solitons for different experiments
abstract type ContactData end
abstract type Sociopatterns <: ContactData end
struct InVS15 <: Sociopatterns end
struct InVS13 <: Sociopatterns end
struct LyonSchool <: Sociopatterns end
struct Copenhagen <: ContactData end

label(experiment::InVS15) = "InVS15"
label(experiment::InVS13) = "InVS13"
label(experiment::LyonSchool) = "LyonSchool"
label(experiment::Copenhagen) = "Copenhagen"

timestep(experiment::InVS15) = 20
timestep(experiment::Copenhagen) = 5*60

duration(experiment::InVS15)     = 14*24*60*60
duration(experiment::Copenhagen) = 28*24*60*60

# A sketch
###############################################################################
#
# Interactions
# ID 1 .........|||.....|||..................|||||......
# ID 2 .........|||..........................|||........
# ID 3 .................|||....................|||......
#
# Interaction Activity
# ID 1 .........111.....111..................11211......
# ID 2 .........111..........................112........
# ID 3 .................111....................211......
#
# Contact Train
# ID 1 .........|.......|....................|.|........
# ID 2 .........|............................|..........
# ID 3 .................|......................|........
#
###############################################################################

###############################################################################
###############################################################################
###############################################################################
# functions that load objects from filtered data

"""
    load_processed_data(experiment, minimum_duraton, path [, filter_out_incomplete])

loads contact data for `experiment` (default Copenhagen()) from `path`, filters
for co-location times that persist a `minimum_duration` (in units of
experiment) and are close enough (specified for each experiment), and returns
the objects that carry the contact activity, the encoutners trains and the list
of contacts.

For the Copenhagen Network Study, one can choose to filter out those
individuals for which no signal was recording on neither first nor last day of
experiment, indicating that they did not participate the full duration.
"""
function load_processed_data(
        experiment,                   # soliton which experiment, e.g. `Coppenhagen()`
        minimum_duration,             # absolute time (in units of experiment)
        path::String;                 # filepath to load data from
        filter_out_incomplete=false   # if `true`, trains that do not have any signal in the first and last 24h are discarded
    )

    # ids is a list of ids, colocation_times is a 2d array with columns: timestamp, id, jd
    # @CodeRev: this should be clear now from documentation of the function
    selected_ids, colocation_times = load_data(experiment, path, filter_out_incomplete=filter_out_incomplete);
    @printf("Duration: %.2f days\n", (colocation_times[end, 1] - colocation_times[1, 1])/60/60/24.0)

    # get co-location times sorted per tuple
    # `ids_raw` has same format as `ids`, `pairwise_colocation_times` is a dict with tuple (id, jd) as keys and a list of colocation_times as values.
    # call as e.g. colocation_time_pairwise[(0,15)]
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times);

    # get contacts and filter out those that are too short
    list_contacts = get_contacts(pairwise_colocation_times, selected_ids, timestep(experiment))
    filter_out_short_contacts!(list_contacts, 15*60)

    # calculate encoutner trains and contact activity from list of contacts
    ets = encounter_trains(list_contacts, selected_ids, duration(experiment), timestep(experiment));
    cas = contact_activities(list_contacts, selected_ids, 0:timestep(experiment):duration(experiment));

    return cas, ets, list_contacts
end
load_processed_data() = load_processed_data(Copenhagen(), 15*60, "./dat")

###############################################################################
###############################################################################
###############################################################################
# load data

"""
    load_data(experiment, path [, filter_out_incomplete])

loads interaction data for selected `experiment` from `path` and returns
participant ids and 2D array of colocation_times with column (time, id A, id B).

In case of Copenhagen(), the raw data is from bluetooth signals such that raw
data is filtered for close proximity (only bluetooth signals with RSSI>=-80, cf
Fig. 2 in [1]) and only indlcudes colocation_times with valid id B >=0 (discarding
signals from empty scans or with devices outside of study).

In the case of Copenhange(), there is the option to filter out incomplete
participant participation by returning only ids that have any bluetooth signal
recording on both first and last day of the experiment.

# References
[1] https://doi.org/10.1371/journal.pone.0100915).
"""
# we know that data is in integer format
function load_data(experiment::Copenhagen, path; filter_out_incomplete=false)
    @printf("Load dataset '%s' \n", label(experiment))
    # timestamp, user_a, user_b, rssi
    data = readdlm(@sprintf("%s/bt_symmetric.csv", path), ',', Int64, skipstart=1)

    # filter out rows where bluetooth signal strength (rssi) is below -80
    # Reason: * rssi value in data is maximum that occured within 5min window
    #         * for the devices used in this study, the authors studied the
    #         relation between rssi and distance and find that signal for 1m is
    #         distributed around -75 and for 2m around -80 (perfect scenario
    #         Fig.2 in [1]
    #         * if rssi < -80 this means that for 5min no single encoutner closer
    #         than 1.5m
    colocation_times = data[data[:,4].>=-80,:]

    # make dimension match sociopatterns format
    colocation_times = colocation_times[:,1:3]

    # filter out rows where user_b is not a valid target (>=0)
    colocation_times = colocation_times[colocation_times[:,3].>=0,:]

    # find uniqe ids (from user a and user b)
    unique_a = unique!(colocation_times[:,2])
    unique_b = unique!(colocation_times[:,3])
    ids = sort!(unique!(vcat(unique_a, unique_b)))

    # filter out ids for which signals occur neither on first nor on last day (minority)
    # TODO@CodeRev: should this maybe be in an extra function? If so, how?
    if filter_out_incomplete
        print("... filter out ids of participants that had no signal on both first and last day of study\n")
        time_first_occurence = [data[findfirst(x->x==id, data[:,2]),1] for id in ids]
        time_last_occurence = [data[findlast(x->x==id, data[:,2]),1] for id in ids]
        print("... reduced from ", length(ids))
        ids = ids[(time_first_occurence .< 1*24*60*60) .== (time_last_occurence .> 27*24*60*60)]
        print("... to ", length(ids), " ids\n")
    end

    return ids, colocation_times
end
function load_data(experiment::Sociopatterns, path; filter_out_incomplete=false)
    if filter_out_incomplete
        throw(error("filter_out_incomplete is not implemented for experiments from sociopatterns because they do not rely on bluetooth"))
    end
    @printf("Load dataset '%s' \n", label(experiment))
    colocation_times = readdlm(@sprintf("./dat/tij_pres_%s.dat", label(experiment)), ' ', Int64)
    #meta = readdlm(@sprintf("./dat/metadata_%s.dat", label(experiment)), '\t', Any)

    # get ids as array of Ints
    unique_a = unique!(colocation_times[:,2])
    unique_b = unique!(colocation_times[:,3])
    ids = sort!(unique!(vcat(unique_a, unique_b)))
    #ids = meta[:,1]
    #ids = convert(Array{typeof(ids[1])}, ids)

    return ids, colocation_times
end
#default
load_data(;path="./dat/") = load_data(experiment_Copenhagen(), path)

"""
    sort_colocation_times_pairwise(colocation_times)

sorts 2D array of colocation_times with colums (time, id A, id B) into a dictionary
of pairwise colocation_time times and returns a dictionary of
`pairwise_colocation_times`.

access pairwise_colocation_times[(id A, id B)] yields list of times

Because for all our data sets the ids are integer, this implementation yields a
keys of type Tuple{Int64,Int64}. For other data sets where ids may be
non-integer, this function may have to be implemented differenly (depending on
the format of `colocation_times`, which could also be Vector{Tuple{T, I, I}}).
"""
function sort_colocation_times_pairwise(colocation_times::Array{T,2}) where T
    unique_i = unique!(colocation_times[:,2])
    unique_j = unique!(colocation_times[:,3])
    ids = sort!(unique!(vcat(unique_i, unique_j)))

    # initiate dictionary as empty lists
    pairwise_colocation_times = Dict{Tuple{T,T},Vector{T}}()
    for id in ids
        for jd in ids
            pairwise_colocation_times[(id,jd)] = T[]
        end
    end

    # go through data and append times to lists
    # TODO@CodeRev we go through twice here because id1, id2 could be in prnciple also "A", "B", ...
    for (t, id1, id2) in eachrow(colocation_times)
        push!(pairwise_colocation_times[(id1,id2)], t)
        push!(pairwise_colocation_times[(id2,id1)], t)
    end

    return pairwise_colocation_times
end

"""
    filter_out_short_contacts!(list_contacts, minimum_duration)

removes contacts that are shorter than `minimum_duration` from `list_contacts`.

When used with list_contacts from `get_contacts` then the convention of contact
duration is such that `minimum_duration` sets the duration below which
colocation_times get removed:
* if minimum_duration=300  and  times=[300,]        -> duration 0   -> removed
* if minimum_duration=600  and  times=[300,600,900] -> duration 600 -> kept

Operates either directly on a vector of colocation_time times or on a dictionary of
pairwise colocation_time times.
"""
function filter_out_short_contacts!(list_contacts::Vector{Vector{Tuple{T,T}}}, minimum_duration::T ) where {T}
    for contacts in list_contacts
        filter_out_short_contacts!(contacts, minimum_duration)
    end
end
function filter_out_short_contacts!(contacts::Vector{Tuple{T,T}}, minimum_duration::T) where T
    filter!(x->getindex(x,2)>=minimum_duration, contacts)
end

###############################################################################
###############################################################################
###############################################################################
# Interaction activity

"""
   activity

activity objects stores a timeseries of the number of colocation_times of `id` in a
`time_range` (typically 0:timestep:duration) together with a sum of
colocation_times for some of the sampling algorithms.

# Constructors
1) directly with activity(timeseries)
2) from the contacts of one individual with any other (vector of Tuples (time contact start, contact duration))
"""
struct activity
    timeseries::Vector{Int64}
    sum::Int64

    function activity(timeseries::Vector{Int64})
        new(copy(timeseries), sum(timeseries))
    end
end
function activity(contacts, time_range)
    # store encoutner list
    timeseries = zeros(Int64, length(time_range))

    # iterate over all contacts
    for contact in contacts
        contact_start = contact[1]
        contact_end   = contact[1] + contact[2]
        timeseries[searchsortedlast(time_range, contact_start):searchsortedlast(time_range, contact_end)] .+= 1
    end
    return activity(timeseries)
end
timeseries(a::activity) = a.timeseries
Base.sum(a::activity) = a.sum
Base.length(a::activity) = length(a.timeseries)
Base.getindex(a::activity, idx::Real) = a.timeseries[idx]
Base.getindex(a::activity, idxs::AbstractRange) = a.timeseries[idxs]
Base.lastindex(a::activity) = lastindex(a.timeseries)
Base.iterate(a::activity) = iterate(a.timeseries)
Base.iterate(a::activity, idx::Real) = iterate(a.timeseries, idx)

"""
    contact_activities{I,T}

includes a vector of `ids`, the `time_range` of the experiment, a vector of
`activities`

* `activity' is defined as a timeseries of the number of co-locations (per
timestep defined in time_range) for each id.

# Constructors:
1) Directly from a list of `activities`, a list of `ids`, and a `time_range`
2) From a dictionary of `pairwise_colocation_times`, a lust of `ids` that has to be
a subset of the ids defined in `pairwise_colocation_times`, and a `time_range`,
which additionally calculates the list of activities by counting the number of
colocation_times for each id in the specified time_range (e.g. 0:timestep:duration)

# Access
* activities(cas) returns list of activities
* random_colocation_time(rng, cas) returns (index_id, index_activity) to access activities(cas)[index_train][index_time]
* time_range(cas) returns time_range of the experiemnt
* ids(cas) returns ids of the experiemnt
* cas[1:3] returns a copy of the object with only the first three activitiess and ids
* cas[BitArray] returns a copy of the object with BitArray applied to activities an ids
"""
struct contact_activities{I,T}
    _activities::Vector{activity}     # outer: id, inner: activity (time series)
    _ids::Vector{I}
    _time_range::AbstractRange{T}     # e.g. 1:3:100

    # minimal constructor that processes already calculated activities
    function contact_activities(activities::Vector{activity}, ids::Vector{I}, time_range::AbstractRange{T}) where {I,T}
        @assert length(activities) == length(ids)
        @assert length(activities[1]) == length(time_range)
        @assert length(activities[end]) == length(time_range)
        new{I,T}(activities, ids, time_range)
    end
end
# constructor that creates activites from pairwise colocation times and passes
# those to the inner constructor
function contact_activities(list_contacts::Vector{Vector{Tuple{T,T}}}, ids::Vector{I}, time_range::AbstractRange{T}) where {I,T}
    #this is used because `ids` can be a selection of all potential ids
    activities = activity[]
    @showprogress 1 for contacts in list_contacts
        push!(activities, activity(contacts, time_range))
    end

    contact_activities(activities, ids, time_range)
end

###############################################################################
# access functions
# TODO: bound checks! -> see how this is done in Vectors??
Base.getindex(cas::contact_activities{I,T}, idxs::Union{AbstractRange, BitArray}) where {I,T} = contact_activities(cas._activities[idxs], cas._ids[idxs], cas._time_range)
#
ids(cas::contact_activities{I,T}) where {I,T} = cas._ids
#
activities(cas::contact_activities) = cas._activities
#
times(cas::contact_activities) where {I,T} = cas._time_range
#
Base.length(cas::contact_activities{I,T}) where {I,T} = length(cas._activities)
#
Base.lastindex(cas::contact_activities{I,T}) where {I,T} = lastindex(cas._activities)

# get random contact (e.g. for disease start) -> replaces contact_pointer
random_contact(rng::AbstractRNG, cas::contact_activities{I,T}) where {I,T} = random_contact(rng::AbstractRNG, activities(cas))
function random_contact(rng::AbstractRNG, activities::Vector{activity})
    sum_contacts_target = rand(rng, 1:sum(sum.(activities)))

    # find index_a, index_t such that
    # sum(sum.(activities[1:index_a-1])) + sum(activitites[index_a][1:index_t] >= sum_contacts_target
    sum_contacts = 0
    for (index_a, activity) in enumerate(activities)
        # find index_a such that sum_contacts_target is reached inside
        # @CodeRev: sum(activity) is not calculating the sum each time, it is
        # accessing a value in memory
        if sum_contacts + sum(activity) >= sum_contacts_target
            for (index_t, num_contacts) in enumerate(activity)
                sum_contacts += num_contacts
                # check for >= because num_contacts can be > 1 (thereby
                # adding all degenerate interactiosn can increase
                # sum_contacts above target)
                if  sum_contacts >= sum_contacts_target
                    return (index_a, index_t)
                end
            end
        end
        sum_contacts += sum(activity)
    end

    throw(error("something went wrong when locating a random contact"))
end

###############################################################################
###############################################################################
###############################################################################
# encounter (for encounter duration distribution)

"""
    get_contacts(pairwise_colocation_times, id, separation_threshold)

calculate and returns a list of `contacts` as a vectors of contacts (time
start, duration) for each specific `id` with any other `jd`.  The contacts are
calculated by merging co-location times between id and jd that are separated
apart shorter than separation_threshold into a single encounter.

If instead of a single id a vector of `ids` is provided, this function returns
Vector{Vector{encounter}}

The duration is here defined as the time-difference between the first and the
last time, such that isolated (single timestep) colocation_time times have duration 0.
Examples for Copenhagen() with timestep=300s:
* times=[300,]        -> duration 0
* times=[300,600,900] -> duration 600

"""
function get_contacts(pairwise_colocation_times::Dict{Tuple{I,I},Vector{T}}, ids::Vector{I}, separation_threshold::T) where {I,T}
    # this is needed because `ids` can be a selection of all ids
    all_ids = unique(getindex.(keys(pairwise_colocation_times), 1))
    list_contacts = Vector{Tuple{T, T}}[]
    @showprogress 1 for id in ids
        encounter = get_contacts(pairwise_colocation_times, all_ids, id, separation_threshold)
        push!(list_contacts, encounter)
    end
    return list_contacts
end
function get_contacts(pairwise_colocation_times::Dict{Tuple{I,I},Vector{T}}, ids::Vector{I}, id::I, separation_threshold::T) where {I,T}
    @assert id in ids

    # store encoutner list
    c_list = Tuple{T,T}[]

    # iterate over all possible targets (jd)
    for jd in ids
        # exclude self from targets
        if (jd==id)
            continue
        end

        # get encoutner list for tuple from the co-location times
        colocation_times = pairwise_colocation_times[(id,jd)]
        #TODO@CodeRev: this is called twice! once for id and once for jd. okay or not?
        _contacts_per_tuple!(c_list, colocation_times, separation_threshold)
    end

    #sort and return
    sort!(c_list)
    return c_list
end
#internal use only
#
#fill contact list `c_list` [list of (begin, duration)] from the
#colocation times of pairs (id, jd). A contact last as long as subsequent
#colocation times are seperated at most by `separation_threshold`.
#
#Remark
#The duration for contacts with n time steps is (n-1)*timestep, such that for
#isolated interaction time the duration is 0. This convention considers that an
#interaction in a single timestep has vanishing probability to last exactly for
#1 time step.
function _contacts_per_tuple!(c_list::Vector{Tuple{T, T}}, colocation_times::Vector{T}, separation_threshold::T) where {I,T}
    if (length(colocation_times) == 0)
        return
    end

    # variables in time units of data
    dt = 0

    e_beg = colocation_times[1]

    for i in 1:length(colocation_times)
        c_end = colocation_times[i]
        # end if last colocation time
        if i==length(colocation_times)
            push!(c_list, (e_beg, c_end-e_beg))
        else
            # check for dt between colocation times
            dt = colocation_times[i+1] - colocation_times[i]
            # if two encounter separated long enough then save current encounter
            # to list and reset beginning to the next encounter time
            if dt > separation_threshold
                push!(c_list, (e_beg, c_end-e_beg))
                # reset beginning
                e_beg = colocation_times[i+1]
            end
        end
    end
end

###############################################################################
###############################################################################
###############################################################################
# Contact trains

"""
    encounter_trains{T}

includes a vector of `ids`, a list of encoutner `trains`, the `duration` and
`timestep` of the experiment

* a single `train' for user id is defined as the list of all encoutner times
(start times of encounter, no information about duration) that user id had
during experiment.

# Constructors:
1) Directly from a list of `trains`{list of lists), a list of `ids`(list), a
`duration`(T), and a `timestep` (T)
2) From a list of contacts
2) From a dictionary of `pairwise_colocation_times` (Dict), a list of `ids` that has to be
a subset of the ids defined in `pairwise_colocation_times`, and a `time_range`,
which additionally calculates the list of activities by counting the number of
co-location times for each id in the specified time_range (e.g. 0:timestep:duration)

# Access
* trains(ets) returns list of encoutner trains
* random_encoutner(ets) returns (index_train, index_time) to access trains(ets)[index_train][index_time)
* duration(ets) returns duration of experiment
* timestep(ets) returns timestep of experiment
* ids(ets) returns ids of experiment
* ets[1:3] returns a copy of the object with only the first three contat trains and ids
* ets[BitArray] returns a copy of the object with BitArray applied to trains an ids

"""
struct encounter_trains{I,T}
    _trains::Vector{Vector{T}}  # nested list: outer for participants, inner for the encoutner times of the participant
    _ids::Vector{I}             # 1d list of the ids of participants
    _duration::T                # duration of the experiment
    _timestep::T                # time resolution of the recording
    _cumsum_encounter_for_random_selection::Vector{Int}

    # minimal constructor that processes already calculated trains
    function encounter_trains(trains::Vector{Vector{T}}, ids::Vector{I}, duration, timestep) where {I,T}
        cumsum_encounter = cumsum(length.(trains))
        new{I,T}(trains, ids, duration, timestep, cumsum_encounter)
    end
end
# constructor that extracts trains from list of encounters and passes those to
# the inner constructor
function encounter_trains(list_contacts::Vector{Vector{Tuple{T,T}}}, ids::Vector{I}, duration::T, timestep::T) where {I,T}
    trains = Vector{T}[]
    for contact in list_contacts
        # encounter is a list of tuples (time start, duration) and this
        # extracts only time
        push!(trains, getindex.(contact,1))
    end
    encounter_trains(trains, ids, duration, timestep)
end
# constructor that creates lists of encounters from pairwise co-location times,
# extracts trains, and passes those to the next constructor
function encounter_trains(pairwise_colocation_times::Dict{Tuple{I,I},Vector{T}}, ids::Vector{I}, duration::T, timestep::T) where {I,T}
    @assert length(pairwise_colocation_times) == length(ids)^2

    trains = Vector{T}[]
    @showprogress 1 for id in ids
        train = getindex.(get_contacts(pairwise_colocation_times, ids, id, timestep),1)
        push!(trains, train)
    end

    encounter_trains(trains, ids, duration, timestep)
end

"""
    trains_from_dts(list_dts::Vector{Vector{T}}) where T

convenience function to map a vector (over ids) of vectors of inter-encoutner
times to a vector (over ids) of encoutner times. Because time of first encoutner is
unknown, the most simple convention of encoutner at time zero is assumed.
"""
function trains_from_dts(list_dts::Vector{Vector{T}}) where T
    trains = Vector{T}[]
    for dts in list_dts
        push!(trains, pushfirst!(cumsum(dts),0))
    end
    return trains
end

###############################################################################
# access functions
Base.getindex(ets::encounter_trains{I,T}, idxs::Union{AbstractRange, BitArray}) where {I,T} = encounter_trains(ets._trains[idxs], ets._ids[idxs], ets._duration, ets._timestep)
#
trains(ets::encounter_trains{I,T}) where {I,T} = ets._trains
#
ids(ets::encounter_trains{I,T}) where {I,T} = ets._ids
#
duration(ets::encounter_trains{I,T}) where {I,T} = ets._duration
timestep(ets::encounter_trains{I,T}) where {I,T} = ets._timestep
#
Base.length(ets::encounter_trains{I,T}) where {I,T} = length(trains(ets))
#
Base.lastindex(ets::encounter_trains{I,T}) where {I,T} = lastindex(trains(ets))
# this line makes broadcasting possible by rediricting it to the trains
Broadcast.broadcastable(ets::encounter_trains{I,T}) where {I,T} = trains(ets)

# get random interaction (e.g. for disease start)
#checked to yield identical result as old version
function random_encoutner(rng::AbstractRNG, ets::encounter_trains{I,T}) where {I,T}
    num_encounters_target = rand(rng, 1:ets._cumsum_encounter_for_random_selection[end])

    # find index_train from cumulative number encounters over trains
    index_train = searchsortedfirst(ets._cumsum_encounter_for_random_selection, num_encounters_target)
    if index_train == 1
        num_encounters_past = 0
    else
        num_encounters_past = ets._cumsum_encounter_for_random_selection[index_train-1]
    end

    return (index_train, num_encounters_target - num_encounters_past)
end
#function random_encoutner(rng::AbstractRNG, trains::Vector{Vector{T}}) where T
#    num_encoutners_target = rand(rng, 1:sum(length.(trains)))
#
#    # find index_train, index_time such that
#    num_encoutners = 0
#    for (index_train, train) in enumerate(trains)
#        # find the train in which the num_encoutners_target will be reached (index = [1:length])
#        if num_encoutners + length(train) >= num_encoutners_target
#            return (index_train, num_encoutners_target - num_encoutners)
#        end
#        num_encoutners += length(train)
#    end
#
#    throw(error("something went wrong when locating a random encoutner"))
#end

###############################################################################
###############################################################################
###############################################################################
# tests

function test_encounter_train_vs_contact_activity(;
        experiment=Copenhagen(),
        path="./dat",
    )
    # contact activity
    selected_ids, colocation_times = load_data(experiment, path)
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)
    filter_out_short_contacts!(pairwise_colocation_times, timestep(Copenhagen()), 15*60)
    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    # encoutner train
    ets = encounter_trains(pairwise_colocation_times, selected_ids, duration(experiment), timestep(experiment))

    valid = true

    # test if same number of trains are zero
    println("number of zero trains: ", sum(sum.(activities(cas)).== 0) == sum(length.(trains(ets)).== 0))
    valid &= (sum(sum.(activities(cas)).== 0) == sum(length.(trains(ets)).== 0))

    return valid
end

function run_tests()
end

###############################################################################
###############################################################################
###############################################################################
# pauls tinkering

"""
ensure that for every id, in the time_bins where activity is > 0, we also have an entry
in the raw (and filtered) co-location times.

* repeat for activity on the actual raw data
* repeat for encoutner trains
* repeat for encounters ->  is there an easy way to get our encounter representation in code?
    for id in ids:
        for encounter in get_encounters_for_id:
            * check that start time has signal in raw data
            * check that end time has signal in raw data
            * check that _some time between start and end_ has signal in raw data
"""
function test_activity_vs_colocation(;
        experiment=Copenhagen(),
        path="./dat",
    )
    selected_ids, colocation_times = load_data(experiment, path)
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)

    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    filter_out_short_contacts!(pairwise_colocation_times, timestep(Copenhagen()), 15*60)
    cas_filtered = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    # some trivial sanity checks
    @assert length(ids(cas)) >= length(ids(cas_filtered))
    # @assert that for every id, the sum of activity is <= after filtering
    for (index, id) in enumerate(ids(cas))
        jndex = findall(x->x==id, ids(cas_filtered))[1]
        @assert sum(activities(cas)[index]) >= sum(activities(cas_filtered)[jndex])
    end

    # go over every timestep and check that id is (not) in the raw data
    @showprogress for (tdx, ts) in enumerate(0:timestep(experiment):duration(experiment))
        to_check = colocation_times[colocation_times[:,1].==ts,:]
        # check the unfiltered case
        for (index, id) in enumerate(ids(cas))
            if activities(cas)[index][tdx] == 0
                @assert !(id in to_check[:,2] || id in to_check[:,3])
            else
                @assert (id in to_check[:,2] || id in to_check[:,3])
            end
        end

        # check after filtering, here we cannot check for 0 activity
        for (index, id) in enumerate(ids(cas_filtered))
            if activities(cas_filtered)[index][tdx] > 0
                @assert (id in to_check[:,2] || id in to_check[:,3])
            end
        end
    end
end

"""
simple, artificial cases to test api works as expected
"""
struct Dummy <: ContactData end
label(experiment::Dummy)     = "Dummy"
timestep(experiment::Dummy)  = 500
duration(experiment::Dummy)  = 20*500           # leads to 21 timesteps

function test_dummy_cases()

    experiment = Dummy()
    ts = timestep(experiment)
    passed = true

    selected_ids = [1,2,3]
    colocation_times = [
        1*ts 1 2;
        2*ts 1 2;
        3*ts 1 2;
        5*ts 2 1;   # below separation threshold, but long enough
        6*ts 2 1;
        11*ts 1 2;  # above threshold
        12*ts 1 2;
        12*ts 1 3;  # two encoutners -> activity 2
        12*ts 2 3;
        13*ts 1 2;
        13*ts 1 3;
        13*ts 2 3;
        14*ts 1 2;
        14*ts 1 3;
        14*ts 2 3;
        15*ts 1 2;
        15*ts 1 3;
        15*ts 2 3;
        19*ts 2 3;  # filtered because too short
    ]

    println("\nBefore Filtering:")
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)
    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    expected = Dict()
    expected[1] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    [0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,0,0]
    expected[2] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    [0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,1,0]
    expected[3] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    [0,0,0,0,0,0,0,0,0,0,0,0,2,2,2,2,0,0,0,1,0]

    for (idx, i) in enumerate(ids(cas))
        println("id: ", i, "\nExpected | Found")
        for t in 1:length(expected[idx])
            println(expected[idx][t], " | ", activities(cas)[idx][t])
        end
        if timeseries(activities(cas)[idx]) == expected[idx]
            println("check")
        else
            println("failed")
            passed = false
        end
    end

    println("\nFiltering Nothing:")
    # dummy filter, filtering for zero should not change data
    # minimum duration of 0 -> no filtering
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)
    filter_out_short_encounter!(pairwise_colocation_times, 1*ts, 0*ts)
    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    for (idx, i) in enumerate(ids(cas))
        println("id: ", i, "\nExpected | Found")
        for t in 1:length(expected[idx])
            println(expected[idx][t], " | ", activities(cas)[idx][t])
        end
        if timeseries(activities(cas)[idx]) == expected[idx]
            println("check")
        else
            println("failed")
            passed = false
        end
    end

    println("\nFiltering duration 1step:")
    # now filter > 1 timesteps length -> only interactions longer or equal to 2 timesteps are kept
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)
    filter_out_short_encounter!(pairwise_colocation_times, 1*ts, 1*ts)
    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    expected = Dict()
    expected[1] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,0,0 # unfiltered
    [0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,0,0]
    expected[2] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,1,0] # unfiltered
    [0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,0,0] # filter the length-1 event at the end
    expected[3] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,0,0,0,0,0,0,0,0,0,0,0,2,2,2,2,0,0,0,1,0] # unfiltered
    [0,0,0,0,0,0,0,0,0,0,0,0,2,2,2,2,0,0,0,0,0] # filter the length-1 event at the end

    for (idx, i) in enumerate(ids(cas))
        println("id: ", i, "\nExpected | Found")
        for t in 1:length(expected[idx])
            println(expected[idx][t], " | ", activities(cas)[idx][t])
        end
        if timeseries(activities(cas)[idx]) == expected[idx]
            println("check")
        else
            println("failed")
            passed = false
        end
    end

    println("\nFiltering duration 2step:")
    # now filter > 2 timesteps length -> only interactions longer or equal to 3 timesteps are kept
    pairwise_colocation_times = sort_colocation_times_pairwise(colocation_times)
    filter_out_short_encounter!(pairwise_colocation_times, 1*ts, 2*ts)
    cas = contact_activities(pairwise_colocation_times, selected_ids, 0:timestep(experiment):duration(experiment))

    expected = Dict()
    expected[1] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,0,0 # unfiltered
    [0,1,1,1,0,0,0,0,0,0,0,1,2,2,2,2,0,0,0,0,0]
    expected[2] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,1,1,1,0,1,1,0,0,0,0,1,2,2,2,2,0,0,0,1,0] # unfiltered
    [0,1,1,1,0,0,0,0,0,0,0,1,2,2,2,2,0,0,0,0,0] # also filter the length-2 event in the beginning
    expected[3] =
    #0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
    #0,0,0,0,0,0,0,0,0,0,0,0,2,2,2,2,0,0,0,1,0] # unfiltered
    [0,0,0,0,0,0,0,0,0,0,0,0,2,2,2,2,0,0,0,0,0] # also filter the length-2 event in the beginning

    for (idx, i) in enumerate(ids(cas))
        println("id: ", i, "\nExpected | Found")
        for t in 1:length(expected[idx])
            println(expected[idx][t], " | ", activities(cas)[idx][t])
        end
        if timeseries(activities(cas)[idx]) == expected[idx]
            println("check")
        else
            println("failed")
            passed = false
        end
    end

    if passed
        println("\nAll checks out")
    else
        println("\nTest failed")
    end

    return passed

end


text:
__precompile__()

module KVS
# Module Initializer
# function __init__()
# end

# Depending modules

# Exported
export KVSLevelDb, kvs_open, kvs_close, kvs_put, kvs_get, kvs_list_channels, kvs_delete, kvs_delete_range
export kvs_atomic_update, kvs_atomic_batch_increment, kvs_wait_for_updaters
export kvs_atomic_try_update, kvs_put_chunk, byte_array, array_to_type
export AbstractKVS, kvs_clobber!, kvs_path_append, kvs_get_many
export atomic_op, ATOMIC_NONE, ATOMIC_INC, ATOMIC_DEC
export KVSRocksDB, kvs_put_sync
export kvs_create_checkpoint, kvs_get_db_root, kvs_write_batch
export KVSRiak, KVSConsul, kvs_put_atomic, kvs_get_index

# Declarations
abstract type AbstractKVS end
@enum atomic_op ATOMIC_NONE ATOMIC_INC ATOMIC_DEC # Passed to atomic_update()

# Source Files
include("KVSDeclarations.jl")
include("KVSLevelDb.jl")
include("KVSRocksDB.jl")
include("KVSRiak.jl")
include("KVSConsul.jl")
# Deprecated

end # module KVS


text:
module TiffImages

using ColorTypes
using DocStringExtensions
using FileIO
using FixedPointNumbers
using IndirectArrays
using OffsetArrays
using OrderedCollections
using PkgVersion
using ProgressMeter
using Base.Iterators

const PKGVERSION = @PkgVersion.Version 0

include("enum.jl")
include("utils.jl")
include("files.jl")
include("compression.jl")
include("tags.jl")
include("ifds.jl")
include("layout.jl")
include(joinpath("types", "common.jl"))
include(joinpath("types", "dense.jl"))
include(joinpath("types", "mmap.jl"))
include("load.jl")

end # module


text:
using EMIRT
using EMIRT.SegmentMSTs
using EMIRT.IOs
using Watershed
using Test

@testset "test segmentation with mst" begin
aff = readaff(joinpath(dirname(@__FILE__),"../assets/piriform.aff.h5"))

seg, rg = watershed(aff; is_threshold_relative=true);
segmentPairs, segmentPairAffinities = rg2segmentPairs(rg)
sgm = SegMST(seg, segmentPairs, segmentPairAffinities)

println("merge mst: ")
@time merge!(sgm, 0.5);

# include(joinpath(Pkg.dir(), "EMIRT/plugins/emshow.jl"))
# show(sgm.seg)
#save(joinpath(dirname(@__FILE__),"../assets/sgm.h5", sgm)

end # end of test set


text:
abstract type AbstractLayer end

struct WeightedLayer{T} <: AbstractLayer
    scaling::T
end

"""
    WeightedInput(scaling)
    WeightedInput(;scaling=0.1)

Returns a weighted layer initializer object, that will produce a weighted input matrix with
a with random non-zero elements drawn from [-```scaling```, ```scaling```], as described
in [1]. The ```scaling``` factor can be given as arg or kwarg.

[1] Lu, Zhixin, et al. "_Reservoir observers: Model-free inference of unmeasured variables
in chaotic systems._"
Chaos: An Interdisciplinary Journal of Nonlinear Science 27.4 (2017): 041102.
"""
function WeightedLayer(; scaling=0.1)
    return WeightedLayer(scaling)
end

function create_layer(
    input_layer::WeightedLayer,
    approx_res_size,
    in_size;
    matrix_type = Matrix{Float64}
)
    scaling = input_layer.scaling
    res_size = Int(floor(approx_res_size/in_size)*in_size)
    layer_matrix = zeros(res_size, in_size)
    q = floor(Int, res_size/in_size)

    for i in 1:in_size
        layer_matrix[(i-1)*q+1:(i)*q, i] = rand(Uniform(-scaling, scaling), 1, q)
    end

    return Adapt.adapt(matrix_type, layer_matrix)
end

function create_layer(layer, args...; kwargs...)
    return layer
end

"""
    DenseLayer(scaling)
    DenseLayer(;scaling=0.1)

Returns a fully connected layer initializer object, that will produce a weighted input
matrix with a with random non-zero elements drawn from [-```scaling```, ```scaling```].
The ```scaling``` factor can be given as arg or kwarg. This is the default choice in the
```ESN``` construction.
"""
struct DenseLayer{T} <: AbstractLayer
    scaling::T
end

function DenseLayer(; scaling=0.1)
    return DenseLayer(scaling)
end

"""
    create_layer(input_layer::AbstractLayer, res_size, in_size)

Returns a ```res_size``` times ```in_size``` matrix layer, built accordingly to the
```input_layer``` constructor.
"""
function create_layer(
    input_layer::DenseLayer,
    res_size,
    in_size;
    matrix_type = Matrix{Float64}
)
    scaling = input_layer.scaling
    layer_matrix = rand(Uniform(-scaling, scaling), res_size, in_size)
    return Adapt.adapt(matrix_type, layer_matrix)
end

"""
    SparseLayer(scaling, sparsity)
    SparseLayer(scaling; sparsity=0.1)
    SparseLayer(;scaling=0.1, sparsity=0.1)

Returns a sparsely connected layer initializer object, that will produce a random sparse
input matrix with random non-zero elements drawn from [-```scaling```, ```scaling```] and
given sparsity. The ```scaling``` and ```sparsity``` factors can be given as args or kwargs.
"""
struct SparseLayer{T} <: AbstractLayer
    scaling::T
    sparsity::T
end

function SparseLayer(; scaling=0.1, sparsity=0.1)
    return SparseLayer(scaling, sparsity)
end

function SparseLayer(scaling_arg; scaling=scaling_arg, sparsity=0.1)
    return SparseLayer(scaling, sparsity)
end

function create_layer(
    input_layer::SparseLayer,
    res_size,
    in_size;
    matrix_type = Matrix{Float64}
)
    layer_matrix = Matrix(sprand(res_size, in_size, input_layer.sparsity))
    layer_matrix = 2.0 .*(layer_matrix.-0.5)
    replace!(layer_matrix, -1.0=>0.0)
    layer_matrix = input_layer.scaling .*layer_matrix
    return Adapt.adapt(matrix_type, layer_matrix)
end

#from "minimum complexity echo state network" Rodan
#and "simple deterministically constructed cycle reservoirs with regular jumps"
#by Rodan and Tino
struct BernoulliSample{T}
    p::T
end

"""
    BernoulliSample(p)
    BernoulliSample(;p=0.5)

Returns a Bernoulli sign constructor for the ```MinimumLayer``` call. The ```p``` factor
determines the probability of the result as in the Distributions call. The value can be
passed as an arg or kwarg. This sign weight determination for input layers is introduced
in [1].

[1] Rodan, Ali, and Peter Tino. "_Minimum complexity echo state network._"
IEEE transactions on neural networks 22.1 (2010): 131-144.
"""
function BernoulliSample(;p=0.5)
    return BernoulliSample(p)
end

struct IrrationalSample{K}
    irrational::Irrational
    start::K
end

"""
    IrrationalSample(irrational, start)
    IrrationalSample(;irrational=pi, start=1)

Returns an irrational sign contructor for the '''MinimumLayer''' call. The values can be
passed as args or kwargs. The sign of the weight are decided from the decimal expansion of
the given ```irrational```. The first ```start``` decimal digits are thresholded at 4.5,
then the n-th input sign will be + and - respectively.

[1] Rodan, Ali, and Peter Tiňo. "_Simple deterministically constructed cycle reservoirs
with regular jumps._" Neural computation 24.7 (2012): 1822-1852.
"""
function IrrationalSample(;irrational=pi, start=1)
    return IrrationalSample(irrational, start)
end

struct MinimumLayer{T,K} <: AbstractLayer
    weight::T
    sampling::K
end

"""
    MinimumLayer(weight, sampling)
    MinimumLayer(weight; sampling=BernoulliSample(0.5))
    MinimumLayer(;weight=0.1, sampling=BernoulliSample(0.5))

Returns a fully connected layer initializer object. The matrix constructed with this
initializer presents the same absolute weight value, decided by the ```weight``` factor.
The sign of each entry is decided by the ```sampling``` struct. Construction detailed
in [1] and [2].

[1] Rodan, Ali, and Peter Tino. "_Minimum complexity echo state network._"
IEEE transactions on neural networks 22.1 (2010): 131-144.
[2] Rodan, Ali, and Peter Tiňo. "_Simple deterministically constructed cycle reservoirs
with regular jumps._" Neural computation 24.7 (2012): 1822-1852.
"""
function MinimumLayer(weight; sampling=BernoulliSample(0.5))
    return MinimumLayer(weight, sampling)
end

function MinimumLayer(; weight=0.1, sampling=BernoulliSample(0.5))
    return MinimumLayer(weight, sampling)
end

function create_layer(
    input_layer::MinimumLayer,
    res_size,
    in_size;
    matrix_type = Matrix{Float64}
)
    sampling = input_layer.sampling
    weight = input_layer.weight
    layer_matrix = create_minimum_input(sampling, res_size, in_size, weight)
    return Adapt.adapt(matrix_type, layer_matrix)
end

function create_minimum_input(sampling::BernoulliSample, res_size, in_size, weight)
    p = sampling.p
    input_matrix = zeros(res_size, in_size)
    for i in 1:res_size
        for j in 1:in_size
            rand(Bernoulli(p)) ? input_matrix[i, j] = weight : input_matrix[i, j] = -weight
        end
    end

    return input_matrix
end

function create_minimum_input(sampling::IrrationalSample, res_size, in_size, weight)
    setprecision(BigFloat, Int(ceil(log2(10)*(res_size*in_size+sampling.start+1))))
    ir_string = string(BigFloat(sampling.irrational)) |> collect
    deleteat!(ir_string, findall(x -> x == '.', ir_string))
    ir_array = zeros(length(ir_string))
    input_matrix = zeros(res_size, in_size)

    for i in 1:length(ir_string)
        ir_array[i] = parse(Int, ir_string[i])
    end

    co = sampling.start

    for i in 1:res_size
        for j in 1:in_size
            ir_array[c] < 5 ? input_matrix[i, j] = -weight : input_matrix[i, j] = weight
            counter += 1
        end
    end

    return input_matrix
end

struct InformedLayer{T,K,M} <: AbstractLayer
    scaling::T
    gamma::K
    model_in_size::M
end

"""
    InformedLayer(model_in_size; scaling=0.1, gamma=0.5)

Returns a weighted input layer matrix, with random non-zero elements drawn from
[-```scaling```, ```scaling```], where some γ of reservoir nodes are connected exclusively
to the raw inputs, and the rest to the outputs of the prior knowledge model,
as described in [1].

[1] Jaideep Pathak et al. "_Hybrid Forecasting of Chaotic Processes: Using Machine Learning
in Conjunction with a Knowledge-Based Model_" (2018)
"""
function InformedLayer(model_in_size; scaling=0.1, gamma=0.5)
    return InformedLayer(scaling, gamma, model_in_size)
end

function create_layer(
    input_layer::InformedLayer,
    res_size,
    in_size;
    matrix_type = Matrix{Float64}
)
    scaling = input_layer.scaling
    state_size = in_size - input_layer.model_in_size

    if state_size <= 0
        throw(DimensionMismatch("in_size must be greater than model_in_size"))
    end

    input_matrix = zeros(res_size, in_size)
    #Vector used to find res nodes not yet connected
    zero_connections = zeros(in_size)
    #Num of res nodes allotted for raw states
    num_for_state = floor(Int, res_size*input_layer.gamma)
    #Num of res nodes allotted for prior model input
    num_for_model = floor(Int, (res_size*(1-input_layer.gamma)))

    for i in 1:num_for_state
        #find res nodes with no connections
        idxs = findall(Bool[zero_connections == input_matrix[i,:]
            for i in 1:size(input_matrix,1)])
        random_row_idx = idxs[rand(1:end)]
        random_clm_idx = range(1, state_size, step = 1)[rand(1:end)]
        input_matrix[random_row_idx,random_clm_idx] = rand(Uniform(-scaling, scaling))
    end

    for i in 1:num_for_model
        idxs = findall(Bool[zero_connections == input_matrix[i,:]
            for i in 1:size(input_matrix,1)])
        random_row_idx = idxs[rand(1:end)]
        random_clm_idx = range(state_size+1, in_size, step = 1)[rand(1:end)]
        input_matrix[random_row_idx,random_clm_idx] = rand(Uniform(-scaling, scaling))
    end

    return Adapt.adapt(matrix_type, input_matrix)
end

"""
    NullLayer(model_in_size; scaling=0.1, gamma=0.5)

Returns a vector of zeros.
"""
struct NullLayer <: AbstractLayer end

function create_layer(
    input_layer::NullLayer,
    res_size,
    in_size;
    matrix_type = Matrix{Float64})
    return Adapt.adapt(matrix_type, zeros(res_size, in_size))
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

module ValidatedNumerics

using Reexport

@reexport using IntervalArithmetic
@reexport using IntervalRootFinding
@reexport using IntervalContractors
@reexport using IntervalConstraintProgramming
@reexport using IntervalOptimisation
@reexport using TaylorModels

end # module ValidatedNumerics


text:
using Phylo
using Phylo: Rootedness, Rooted, TreeType
using Phylo: AbstractNode, AbstractBranch, AbstractTree
using Compat

# Need to be able to create new node and branch labels
function _newlabel(ids::Vector{Label}) where Label <: Integer
    return isempty(ids) ? 1 : maximum(ids) + 1
end
function _newlabel(ids::Vector{Vector{Label}}) where Label <: Integer
    return mapreduce(_newlabel, max, 1, ids)
end
_newlabel(names::Vector{String}, prefix::String) =
    prefix * "$(_newlabelfree(names, prefix))"
_newlabel(names::Vector{Vector{String}}, prefix::String) =
    prefix * "$(mapreduce(_newlabelfree, max, 1, names))"

function _newlabelfree(names::Vector{String}, prefix)
    num = length(names) + 1
    if "$prefix$num" ∉ names
        return num
    end
    len = length(prefix)
    goodnames = [name[(len - 1):end] for name in names
             if length(name) > len && name[1:len] == prefix]
    num = length(goodnames) + 1
    while ("$num" ∈ goodnames)
        num += 1
    end
    return num
end

"""
    _newbranchlabel(tree::AbstractTree)

Returns a new unique branch name for a tree.
"""
function _newbranchlabel end
_newbranchlabel(tree::AbstractTree) = _newlabel(collect(getbranchnames(tree)))

"""
    _newnodelabel(tree::AbstractTree)

Returns a new unique node name for a tree.
"""
function _newnodelabel end
_newnodelabel(tree::AbstractTree{TT, RT, String, N, B}) where
    {TT, RT, N, B} = _newlabel(collect(getnodenames(tree)), "Node ")
_newnodelabel(tree::AbstractTree{TT, RT, I, N, B}) where
    {TT, RT, I <: Integer, N, B} = _newlabel(collect(getnodenames(tree)))

# AbstractTree type methods
"""
    _treenametype(::Type{AbstractTree})

Returns the label type for a tree type. Must be implemented for any tree type.
"""
function _treenametype end
_treenametype(::Type{<: AbstractTree{OneTree}}) = Int

# AbstractTree methods
"""
    _ntrees(::AbstractTree)

Returns the number of trees in an object. Must be implemented for any ManyTrees
type.
"""
function _ntrees end
_ntrees(tree::AbstractTree{OneTree}) = 1
_ntrees(tree::AbstractTree{ManyTrees}) = length(_gettreenames(tree))

"""
    _gettrees(::AbstractTree)

Returns the trees in an object. Must be implemented for any ManyTrees type.
"""
function _gettrees end
_gettrees(tree::AbstractTree{OneTree}) = [tree]

"""
    _gettreenames(::AbstractTree)

Returns the names for the trees. Can be implemented for any ManyTrees type.
"""
function _gettreenames end
_gettreenames(tree::AbstractTree{ManyTrees}) = _gettreename.(_gettrees(tree))

"""
    _gettreename(::AbstractTree)

Returns the name for a single tree. Should be implemented for any OneTree
type where they have names.
"""
function _gettreename end
_gettreename(tree::AbstractTree{OneTree}) = 1

"""
    _gettree(::Pair{Label, AbstractTree})
    _gettree(::AbstractTree, id)

Returns a tree - either itself if it is a single tree, or the single tree
in a set with label id. Must be implemented for any ManyTrees type.
"""
function _gettree end
_gettree(treepair::Pair{Label, <: AbstractTree{OneTree}}) where Label =
    treepair[2]
function _gettree(tree::T, id) where {T <: AbstractTree{OneTree}}
    id isa _treenametype(T) || throw(TypeError(:_gettree, "index argument",
                                               _treenametype(T), id))
    id == _gettreename(tree) || throw(BoundsError(tree, id))
    return tree
end

"""
    _getnodes(tree::AbstractTree)

Returns a vector of nodes for a OneTree tree. Either _getnodes() or
_getnodenames() must be implemented for any OneTree tree type.
"""
function _getnodes end
_getnodes(tree::AbstractTree{OneTree}) =
    [_getnode(tree, name) for name in _getnodenames(tree)]

"""
    _getnodenames(tree::AbstractTree)

Returns a vector of node names for a OneTree tree. Either _getnodes() or
_getnodenames() must be implemented for any OneTree tree type.
"""
function _getnodenames end
_getnodenames(tree::AbstractTree{OneTree}) =
    [_getnodename(tree, node) for node in _getnodes(tree)]

"""
    _nnodes(::AbstractTree)

Returns the number of nodes (internal nodes and leaves) in a single tree. May
be implemented for any OneTree tree type (otherwise infers from _getnodes()).
"""
_nnodes(tree::AbstractTree{OneTree}) = length(_getnodes(tree))

"""
    _nbranches(::AbstractTree)

Returns the number of branches in a single tree. May be implemented for any OneTree tree type (otherwise infers from _getbranches()).
"""
_nbranches(tree::AbstractTree{OneTree}) = length(_getbranches(tree))

"""
    _getleafnames(::AbstractTree)

Returns the leaf names of a tree. May be implemented for any
tree type (otherwise determined from _getnodenames() and _isleaf() functions).
"""
function _getleafnames end
_getleafnames(tree::AbstractTree{OneTree}) =
    [node for node in _getnodenames(tree) if _isleaf(tree, node)]
_getleafnames(tree::AbstractTree{ManyTrees}) =
    _getleafnames(_gettree(tree, first(_gettreenames(tree))))

"""
    _getleaves(::AbstractTree)

Returns the leaves (tips) of a single tree. May be implemented for any
OneTree type (otherwise determined from _getnodes() and _isleaf() functions).
"""
function _getleaves end
_getleaves(tree::AbstractTree{OneTree}) =
    [node for node in _getnodes(tree) if _isleaf(tree, node)]

"""
    _nleaves(::AbstractTree)

Returns the number of leaves (tips) in a tree. May be implemented for any
tree type (otherwise determined from the _getleafnames() function).
"""
_nleaves(tree::AbstractTree) = length(_getleafnames(tree))

"""
    _getroots(::AbstractTree)

Returns the root(s) of a tree. May be implemented for any
OneTree type (otherwise determined from _getnodes() and _isroot() functions).
"""
_getroots(tree::AbstractTree{OneTree, <: Rooted}) =
    [node for node in _getnodes(tree) if _isroot(tree, node)]

"""
    _getroot(::AbstractTree)

Returns the unique root of a rooted tree. May be implemented for any
OneTree type (otherwise determined from _getroots()).
"""
function _getroot end
function _getroot(tree::AbstractTree{OneTree, <: Rooted})
    @assert nroots(tree) == 1 "More than one root for tree " *
    "(found $(length(roots)))"
    return first(_getroots(tree))
end

"""
    _nroots(::AbstractTree)

Returns the number of roots (subtrees) in a OneTree tree. May be implemented
for any ManyRoots type (otherwise infers from _getroots()).
"""
function _nroots end
_nroots(tree::AbstractTree{OneTree, <: Rooted}) = length(_getroots(tree))
_nroots(tree::AbstractTree{OneTree, Unrooted}) = 0

"""
    _getnode(::AbstractTree, id)
    _getnode(id)

Returns the node associated with id (which could be a name, a node or a pair)
from a tree. For some id types, it will be possible to extract the node
without reference to the tree. Must be implemented for any tree and node label
type.
"""
function _getnode end
_getnode(::AbstractTree{OneTree, RT, NL, N}, node::N) where {RT, NL, N} =
    node
_getnode(::AbstractTree{OneTree, RT, NL, N}, pair::Pair{NL, N}) where
    {RT, NL, N} = pair[2]
_getnode(pair::Pair{NL, N}) where {RT, NL, N <: AbstractNode{RT, NL}} = pair[2]

"""
    _getnodename(::AbstractTree, id)
    _getnodename(id)

Returns the name of a node associated with id (which could be a name,
a node or a pair) from a tree. For some id and node types, it will be
able to extract the node name without reference to the tree.
"""
function _getnodename end
_getnodename(::AbstractTree{OneTree, RT, NL, N, B}, pair::Pair{NL, N}) where
    {RT, NL, N, B} = pair[1]
_getnodename(::AbstractTree{OneTree, RT, NL, N, B}, nodename::NL) where
    {RT, NL, N, B} = nodename
_getnodename(pair::Pair{NL, N}) where {RT, NL, N <: AbstractNode{RT, NL}} =
    pair[1]

"""
    _getbranch(::AbstractTree, id)
    _getbranch(id)

Returns the branch associated with id (which could be a name, a branch or a
pair) from a tree. For some id types, it will be possible to extract the branch
without reference to the tree. Must be implemented for any OneTree tree type.
"""
function _getbranch end
_getbranch(::AbstractTree{OneTree, RT, NL, N, B}, branch::B) where
    {RT, NL, N, B} = branch
_getbranch(::AbstractTree{OneTree, RT, NL, N, B}, pair::Pair{Int, B}) where
    {RT, NL, N, B} = pair[2]
_getbranch(pair::Pair{Int, B}) where {RT, NL, B <: AbstractBranch{RT, NL}} =
    pair[2]

"""
    _getbranchname(::AbstractTree, id)
    _getbranchname(id)

Returns the branch associated with id (which could be a name, a branch or a
pair) from a tree. For some id and branch types, it will be able to extract
the branch name without reference to the tree.
"""
function _getbranchname end
_getbranchname(::AbstractTree{OneTree, RT, NL, N, B}, branchname::Int) where
    {RT, NL, N, B} = branchname
_getbranchname(::AbstractTree{OneTree, RT, NL, N, B}, pair::Pair{Int, B}) where
    {RT, NL, N, B} = pair[1]
_getbranchname(pair::Pair{Int, B}) where
    {RT, NL, B <: AbstractBranch{RT, NL}} = pair[1]

"""
    _hasnode(tree::AbstractTree, node[name])

May be implemented for any OneTree tree type.
"""
function _hasnode end
_hasnode(tree::AbstractTree{OneTree, RT, NL, N, B}, nodename::NL) where
    {RT, NL, N, B} = nodename ∈ _getnodenames(tree)
_hasnode(tree::AbstractTree{OneTree, RT, NL, N, B}, node::N) where
    {RT, NL, N, B} = node ∈ _getnodes(tree)

"""
    _getbranches(tree::AbstractTree)

Returns a vector of branches for a OneTree tree. Either _getbranches() or
_getbranchnames() must be implemented for any OneTree tree type.
"""
function _getbranches end
_getbranches(tree::AbstractTree{OneTree, RT, NL, N, B}) where
    {RT, NL, N, B} =
    [_getbranch(tree, name) for name in _getbranchnames(tree)]

"""
    _getbranchnames(tree::AbstractTree)

Returns a vector of branch names for a OneTree tree. Either _getbranches() or
_getbranchnames() must be implemented for any OneTree tree type.
"""
function _getbranchnames end
_getbranchnames(tree::AbstractTree{OneTree}) =
    [_getbranchname(tree, branch) for branch in _getbranches(tree)]

"""
    _hasbranch(tree::AbstractTree, branch::AbstractBranch)
    _hasbranch(tree::AbstractTree, source::AbstractNode, dest::AbstractNode)

Tests whether a branch is present in a tree. source and dest method must and
branch method may be implemented for any OneTree tree type.
"""
function _hasbranch end
_hasbranch(::AbstractTree, _) = false
_hasbranch(tree::AbstractTree{OneTree}, branch::AbstractBranch) =
    branch ∈ _getbranches(tree)
_hasbranch(tree::AbstractTree{OneTree}, branchname::Int) =
    branchname ∈ _getbranchnames(tree)

"""
    _createbranch!(tree::AbstractTree, source, destination,
                   length::Float64, data)

Create a new branch and add it to a tree. Must be implemented for any
AbstractTree subtype.
"""
function _createbranch! end

"""
    _deletebranch!(tree::AbstractTree, branch)

Delete a branch, reoving it from a tree. Must be implemented for any
AbstractTree subtype.
"""
function _deletebranch! end

"""
    _createnode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _createnode! end

"""
    _deletenode!(tree::AbstractTree, nodename)

Must be implemented for any AbstractTree subtype.
"""
function _deletenode! end

"""
    _hasrootheight(::AbstractTree)

"""
function _hasrootheight(::AbstractTree)
    return false
end

"""
    _getrootheight(::AbstractTree)

"""
function _getrootheight end

"""
    _setrootheight!(::AbstractTree, value)

"""
function _setrootheight! end

"""
    _clearrootheight!(::AbstractTree)

"""
function _clearrootheight! end

"""
    _validate(::AbstractTree)

"""
function _validate(::AbstractTree)
    return true
end

"""
    _traversal(tree::AbstractTree, order::TraversalOrder, todo, sofar)

Return an iterable object containing nodes in given order - preorder, inorder,
postorder or breadthfirst
"""
function _traversal end
function _traversal(tree::AbstractTree{OneTree, <: Rooted},
                    order::TraversalOrder, todo = collect(_getroots(tree)),
                    sofar = eltype(todo)[])
    while !isempty(todo)
        if order == Phylo.breadthfirst
            append!(sofar, todo)
            children = eltype(todo)[]
            for node in todo
                append!(children, getchildren(tree, node))
            end
            todo = children
        else
            now = popfirst!(todo)
            children = getchildren(tree, now)
            if isempty(children)
                push!(sofar, now)
            else
                if order == Phylo.preorder
                    push!(sofar, now)
                    sofar = _traversal(tree, order, children, sofar)
                elseif order == Phylo.postorder
                    sofar = _traversal(tree, order, children, sofar)
                    push!(sofar, now)
                elseif order == Phylo.inorder
                    child = popfirst!(children)
                    sofar = _traversal(tree, order, [child], sofar)
                    push!(sofar, now)
                    sofar = _traversal(tree, order, children, sofar)
                end
            end
        end
    end
    return sofar
end

# AbstractNode methods
"""
    _isleaf(tree::AbstractTree, node)

Is the node a leaf? Does not need to be implemented for any node type –
inferred from _outdegree or _degree - unless tree knows which nodes are
leaves and not nodes.
"""
function _isleaf end
_isleaf(tree::AbstractTree{OneTree, <: Rooted}, node) =
    _outdegree(tree, _getnode(tree, node)) == 0
function _isleaf(tree::AbstractTree{OneTree, Unrooted}, node)
    _degree(tree, node) == 0 && return true
    _degree(tree, node) > 1 && return false
    return missing
end

"""
    _isroot(tree::AbstractTree, node)

Is the node a root node of the tree?
"""
function _isroot end
_isroot(tree::AbstractTree{OneTree, <: Rooted, NL, N}, node::N) where {NL, N} =
    !_hasinbound(tree, node)
_isroot(tree::AbstractTree{OneTree, <: Rooted, NL}, node::NL) where NL =
    _isroot(tree, _getnode(tree, node))
_isroot(tree::AbstractTree{OneTree, Unrooted}, node) = false

"""
    _isinternal(tree::AbstractTree, node::AbstractNode)

Is the node internal to the tree?
"""
function _isinternal end
_isinternal(tree::AbstractTree{OneTree, RT, NL}, node::NL) where {RT, NL} =
    _isinternal(tree, _getnode(tree, node))
_isinternal(tree::AbstractTree, node::AbstractNode{<: Rooted}) =
    _outdegree(tree, node) > 0 && _hasinbound(tree, node)
_isinternal(tree::AbstractTree, node::AbstractNode{Unrooted}) =
    _degree(tree, node) < 2 ? false : missing

"""
    _isunattached(tree::AbstractTree, node::AbstractNode)

Does the node currently form its own (sub)tree?
"""
_isunattached(tree::AbstractTree, node) =
    _degree(tree, _getnode(tree, node)) == 0

"""
    _indegree(tree::AbstractTree, node::AbstractNode)

In degree of node. Can be implemented for rooted nodes, otherwise inferred
from _hasinbound.
"""
function _indegree end
_indegree(tree::AbstractTree, node::AbstractNode) = Int(_hasinbound(tree, node))
_indegree(tree::AbstractTree, node::AbstractNode{Unrooted}) =
    _degree(tree, node) == 0 ? 0 : missing

"""
    _hasinboundspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new inbound connection on a node?
"""
function _hasinboundspace end
_hasinboundspace(tree::AbstractTree, node::AbstractNode) =
    !_hasinbound(tree, node)

"""
    _outdegree(tree::AbstractTree, node::AbstractNode)

Out degree of node.
"""
function _outdegree end
_outdegree(tree::AbstractTree{OneTree, <: Rooted}, node::AbstractNode) =
    length(_getoutbounds(tree, node))
_outdegree(tree::AbstractTree{OneTree, Unrooted},
           node::AbstractNode{Unrooted}) =
    _degree(tree, node) == 0 ? 0 : missing

"""
    _hasoutboundspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new outbound connection on a node? Must be implemented if
a node has a limit on the number of outbound connections (eg for a binary tree)
"""
function _hasoutboundspace end

"""
    _hasspace(tree::AbstractTree, node::AbstractNode)

Is there space for a new connection on a node? Must be implemented if
a node has a limit on the number of connections (eg for a binary tree)
"""
function _hasspace end
_hasspace(::AbstractTree{OneTree}, ::AbstractNode) =
    _hasinboundspace(tree, node) | _hasoutboundspace(tree, node)
_hasspace(::AbstractTree{OneTree, Unrooted}, ::AbstractNode) = true

"""
    _degree(tree::AbstractTree, node::AbstractNode)

Degree of node. Must be implemented for Unrooted nodes, otherwise can be
inferred from indegree and outdegree.
"""
function _degree end
_degree(tree::AbstractTree, node::AbstractNode) =
    _indegree(tree, node) + _outdegree(tree, node)

"""
    _hasinbound(tree::AbstractTree, node::AbstractNode)

Must be implemented for any AbstractNode subtype.
"""
function _hasinbound end
_hasinbound(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    _hasinbound(tree, _getnode(tree, name))
_hasinbound(tree::AbstractTree{OneTree, Unrooted},
            node::AbstractNode{Unrooted}) =
    _degree(tree, node) == 0 ? false : missing

"""
    _getinbound(tree::AbstractTree, node::AbstractNode)

Get the inbound connection. Must be implemented for any rooted AbstractNode subtype.
"""
function _getinbound end
_getinbound(::AbstractTree{OneTree, Unrooted}, ::AbstractNode{Unrooted}) =
    error("Can't ask for inbound connection on Unrooted trees")
_getinbound(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    _getinbound(tree, _getnode(tree, name))

"""
    _getparent(tree::AbstractTree, node)

Return the parent node for this node. Can be implemented for Rooted node types.
"""
function _getparent end
_getparent(tree::AbstractTree{OneTree, Unrooted}, _) =
    error("Can't ask for parent of Unrooted tree")
_getparent(tree::AbstractTree{OneTree, RT, NL, N, B}, node::N) where
    {RT <: Rooted, NL, N, B} = _getnode(tree,
                                        _src(tree, _getinbound(tree, node)))
_getparent(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    _getnodename(tree, _getparent(tree, _getnode(tree, name)))

"""
    _addinbound!(tree::AbstractTree, node::AbstractNode, inbound)

Adds a branch to the input of a rooted node. Must be implemented for any rooted
AbstractNode subtype unless this happens when a branch is created.
"""
function _addinbound! end
_addinbound!(::AbstractTree{OneTree, <: Rooted}, ::AbstractNode{Unrooted}, _) =
    error("Unrooted trees don't have inbound connections")

"""
    _removeinbound!(tree::AbstractTree, node::AbstractNode, inbound)

Removes a branch from the input of a rooted node. Must be implemented for any
rooted AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeinbound! end
_removeinbound!(::AbstractTree{OneTree, Unrooted}, ::AbstractNode, _) =
    error("Unrooted trees don't have inbound connections")

"""
    _getoutbounds(tree::AbstractTree, node::AbstractNode)

Returns the outbound connections of a rooted node. Must be implemented for any
rooted AbstractNode subtype.
"""
function _getoutbounds end
_getoutbounds(::AbstractTree{OneTree, <: Rooted}, ::AbstractNode{Unrooted}) =
    error("Unrooted trees don't have outbound connections")

"""
    _addoutbound!(tree::AbstractTree, node::AbstractNode, branch)

Add an outbound branch to a rooted node. Must be implemented for any Rooted
AbstractNode subtype unless this happens when a branch is created.
"""
function _addoutbound! end
_addoutbound!(::AbstractTree{OneTree, Unrooted}, ::AbstractNode{Unrooted}, _) =
    error("Unrooted trees don't have inbound and outbound connections")

"""
    _removeoutbound!(tree::AbstractTree, node::AbstractNode, branch)

Remove an outbound branch from a rooted node. Must be implemented for any
AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeoutbound! end
_removeoutbound!(::AbstractTree{OneTree, Unrooted},
                 ::AbstractNode{Unrooted}, _) =
    error("Unrooted trees don't have outbound connections")

"""
    _getchildren(tree::AbstractTree, node)
    _getchildren(tree::AbstractTree, nodename)

Return the child node(s) for this node. May be implemented for any rooted
AbstractNode subtype.
"""
function _getchildren end
_getchildren(::AbstractTree{OneTree, Unrooted}, _) =
    error("Can't ask for children of Unrooted tree")
_getchildren(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    [_getnodename(tree, _dst(tree, branch))
     for branch in _getoutbounds(tree, _getnode(tree, name))]

_getchildren(tree::AbstractTree{OneTree, RT, NL, N, B}, node::N) where
    {RT <: Rooted, NL, N, B} =
    [_getnode(tree, _dst(tree, branch)) for branch in _getoutbounds(tree, node)]

"""
    _getconnections(tree::AbstractTree, node::AbstractNode)

Returns all of the connections of a node. Must be implemented for any unrooted
AbstractNode subtype, can be inferred from _getinbound and _getoutbounds for
a rooted node.
"""
function _getconnections end
_getconnections(tree::AbstractTree{OneTree}, node::AbstractNode) =
    _hasinbound(tree, node) ?
        append!([_getinbound(tree, node)], _getoutbounds(tree, node)) :
        _getoutbounds(tree, node)
_getconnections(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    _hasinbound(tree, name) ?
        append!([_getinbound(tree, name)], _getoutbounds(tree, name)) :
        _getoutbounds(tree, name)

"""
    _getsiblings(tree::AbstractTree, node::AbstractNode)

Returns all of the siblings of a node. May be implemented for any
AbstractNode subtype, can be inferred from _getparent and _getchildren for
a rooted node or _getconnections for an unrooted node.
"""
function _getsiblings end
_getsiblings(tree::AbstractTree{OneTree, <: Rooted}, node::AbstractNode) =
    _hasinbound(tree, node) ?
        append!([_getparent(tree, node)], _getchildren(tree, node)) :
        _getchildren(tree, node)
_getsiblings(tree::AbstractTree{OneTree, Unrooted}, node::AbstractNode) =
    [_conn(tree, branch, node) for branch in _getconnections(tree, node)]
_getsiblings(tree::AbstractTree{OneTree, <: Rooted, NL}, name::NL) where NL =
    _hasinbound(tree, name) ?
        append!([_getparent(tree, name)], _getchildren(tree, name)) :
        _getchildren(tree, name)
_getsiblings(tree::AbstractTree{OneTree, Unrooted, NL}, name::NL) where NL =
    [_getnodename(tree, _conn(tree, branch, _getnode(tree, name)))
     for branch in _getconnections(tree, name)]

"""
    _addconnection!(tree::AbstractTree, node::AbstractNode, branch)

Add a connection to an unrooted node. Must be implemented for any unrooted
AbstractNode subtype unless this happens when a branch is added.
"""
function _addconnection! end
_addconnection!(::AbstractTree, ::AbstractNode{<: Rooted}, _) =
    error("Direction of connection must be specified for a rooted node")

"""
    _removeconnection!(tree::AbstractTree, node::AbstractNode, branch)

Remove a connection from an unrooted node. Must be implemented for any Unrooted
AbstractNode subtype unless this happens when a branch is deleted.
"""
function _removeconnection! end
_removeconnection!(::AbstractTree, ::AbstractNode{<: Rooted}, _) =
    error("Direction of connection must be specified for a rooted node")

"""
    _hasheight(tree::AbstractTree, nodename)

"""
function _hasheight(::AbstractTree, _)
    return false
end

"""
    _getheight(tree::AbstractTree, nodename)

"""
function _getheight(::AbstractTree, _)
    throw(NullException())
    return NaN
end

"""
    _setheight!(::AbstractTree, nodename, value)

"""
function _setheight!(::AbstractTree, _, value)
    throw(NullException())
    return value
end

# AbstractBranch methods
"""
    _src(branch::AbstractBranch)

Return source node for a branch. Must be implemented for any rooted
AbstractBranch subtype.
"""
function _src end
_src(::AbstractTree, ::AbstractBranch{Unrooted}) =
    error("Unrooted branches do not have in and out connections")

"""
    _dst(branch::AbstractBranch)

Return destination node for a branch. Must be implemented for any rooted
AbstractBranch subtype.
"""
function _dst end
_dst(::AbstractTree, ::AbstractBranch{Unrooted}) =
    error("Unrooted branches do not have in and out connections")

"""
    _conns(tree::AbstractTree, branch::AbstractBranch)

Return a vector of connections for a branch. Must be implemented for any
Unrooted AbstractBranch subtype, otherwise can combine _src and _dst.
"""
function _conns end
_conns(tree::AbstractTree{OneTree, <: Rooted}, branch::AbstractBranch) =
    [_src(tree, branch), _dst(tree, branch)]

"""
    _conn(branch::AbstractBranch, exclude::AbstractNode)

Return the connection for a branch that isn't the `exclude` node. May be
implemented for any Unrooted AbstractBranch subtype, otherwise will use
_conns.
"""
function _conn end
function _conn(tree::AbstractTree{OneTree, RT, NL, N, B},
               branch::B, exclude::N) where {RT <: Rooted, NL,
                                             B <: AbstractBranch{RT, NL},
                                             N <: AbstractNode{RT, NL}}
    other = [node for node in _conns(tree, branch) if node != exclude]
    @assert length(other) == 1 _getnodename(tree, exclude) *
        " is not connected to branch $(_getbranchname(tree, branch))"
    return first(other)
end

"""
    _getlength

Return length of a branch. May be implemented for any AbstractBranch
subtype.
"""
function _getlength end
_getlength(::AbstractTree, ::AbstractBranch) = NaN
_getlength(tree::AbstractTree, bn::Int) = _getlength(tree, _getbranch(tree, bn))

"""
    _leafinfotype(::Type{<:AbstractTree})

Returns the type of the leaf info data.
"""
function _leafinfotype end
_leafinfotype(::Type{<:AbstractTree}) = Nothing

"""
    _noderecordtype(::Type{<:AbstractTree})

Returns the type of the node info data.
"""
function _noderecordtype end
_noderecordtype(::Type{<:AbstractTree}) = Nothing

"""
    _branchrecordtype(::Type{<:AbstractTree})

Returns the type of the branch info data.
"""
function _branchrecordtype end
_branchrecordtype(::Type{<:AbstractTree}) = Nothing

function _getleafinfo end
function _setleafinfo! end
function _getnoderecord end
function _setnoderecord! end
function _getbranchinfo end
function _setbranchinfo! end
function _resetleaves! end


text:
using DashBootstrapComponents

navbar = dbc_navbarsimple(
    children=[
        dbc_navitem(dbc_navlink("Page 1", href="#")),
        dbc_dropdownmenu(
            children=[
                dbc_dropdownmenuitem("More pages", header=true),
                dbc_dropdownmenuitem("Page 2", href="#"),
                dbc_dropdownmenuitem("Page 3", href="#"),
            ],
            nav=true,
            in_navbar=true,
            label="More",
        ),
    ],
    brand="NavbarSimple",
    brand_href="#",
    color="primary",
    dark=true,
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Serialization
using LightGraphs
using SimpleWeightedGraphs
using MetaGraphs
using Bijections
using Statistics
import LightGraphs.Parallel
import Base.Threads.@threads

function save_graph_data(graph::GGraphData)
    serialize(CALCULATED_FOLDER * "GGraphData", graph)
end
function save_graph_data(graph::GGraphData, epsilon, name)
    serialize(CALCULATED_FOLDER * "GGraphData $epsilon $name", graph)
end
function load_graph_data()
    if isfile(CALCULATED_FOLDER * "GGraphData")
        deserialize(CALCULATED_FOLDER * "GGraphData")
    else
        GGraphData()
    end
end
function load_graph_data(epsilon)
    deserialize(CALCULATED_FOLDER * "GGraphData $epsilon")
end
function load_graph_data(epsilon, name)
    deserialize(CALCULATED_FOLDER * "GGraphData $epsilon $name")
end

function get_overall_genotype_distribution(graph_data::GGraphData)
    occurances = Int64[]

    for vertex in 1:nv(graph_data.genotype_graph)
        outgoing = outneighbors(graph_data.genotype_graph, vertex)

        if length(outgoing) == 0
            push!(occurances, 0)
        else
            push!(occurances, sum([graph_data.genotype_graph.weights[u, vertex] for u in outgoing]))
        end
    end

    occurances = occurances ./ sum(occurances)

    return occurances
end

function build_genotype_graph!(graph_data::GGraphData)
    parent_offsprings = deserialize(CALCULATED_FOLDER * "parent_offspring_ocurrances")

    genotype_mapping = Bijection()

    for ((parent, _), _) in parent_offsprings
        if haskey(genotype_mapping, parent) == false
            genotype_mapping[parent] = length(genotype_mapping) + 1
        end
    end

    @info "Mapping created. ($(length(genotype_mapping))) Building graph..."

    genotype_graph = SimpleWeightedDiGraph(length(genotype_mapping))

    count = 0

    for ((parent, offspring), num) in parent_offsprings
        if haskey(genotype_mapping, offspring) == false continue end   # we skip sinks without offsprings

        parent_index = genotype_mapping[parent]
        offspring_index = genotype_mapping[offspring]
        add_edge!(genotype_graph, parent_index, offspring_index, num)

        count += 1
        if count % 1000 == 0
            @info (count / length(parent_offsprings))
        end
    end

    @info "Graph built."

    graph_data.genotype_vertex_mapping = genotype_mapping
    graph_data.genotype_graph = genotype_graph
end

function build_phenotype_graph!(graph_data::GGraphData, min_occurances, tolerance, min_samples, max_samples)
    phenotype_graph = SimpleWeightedDiGraph(nv(graph_data.genotype_graph))

    @info "Collect relevant edges."

    edges_to_test = []

    for edge in edges(graph_data.genotype_graph)
        u, v = src(edge), dst(edge)

        if u == v continue end
        if (v, u) in edges_to_test continue end
        if has_edge(graph_data.genotype_graph, u, v) == false continue end
        if graph_data.genotype_graph.weights[v, u] < min_occurances continue end

        push!(edges_to_test, (u, v))
    end

    @info "Calculate Phenotype Similarities for $(length(edges_to_test)) edges."

    num_edges_to_test = length(edges_to_test)
    phenotype_similarities = SharedArrays.SharedArray{Float64}(num_edges_to_test)

    done = Threads.Atomic{Int}(0)

    @threads for (i, (u, v)) in unique(enumerate(edges_to_test))
        genotype_u, genotype_v = graph_data.genotype_vertex_mapping(u), graph_data.genotype_vertex_mapping(v)
        similarity = get_phenotype_similarity(genotype_u, genotype_v, tolerance, min_samples, max_samples)

        phenotype_similarities[i] = similarity

        Threads.atomic_add!(done, 1)
        @info (done[] / length(edges_to_test))
    end

    phenotype_similarities_dict = Dict()

    for i in 1:num_edges_to_test
        similarity = phenotype_similarities[i]
        u, v = edges_to_test[i]

        phenotype_similarities_dict[(u, v)] = similarity
        phenotype_similarities_dict[(v, u)] = similarity
    end

    @info "Build Phenotype Similarities Graph."

    for edge in edges(graph_data.genotype_graph)
        u, v = src(edge), dst(edge)

        if haskey(phenotype_similarities_dict, (u, v)) == false continue end

        similarity = max(eps(), phenotype_similarities_dict[(u, v)])
        add_edge!(phenotype_graph, u, v, similarity)
    end

    graph_data.phenotype_graph = phenotype_graph
end

function build_neutral_networks!(graph_data::GGraphData, epsilon)
    undirected_phenotype_graph = get_simple_weighted_graph(graph_data.phenotype_graph)

    neutral_networks = []
    already_assigned_nodes = []

    # find all neutral networks

    for (i, edge) in edges(undirected_phenotype_graph) |> enumerate
	    if has_edge(undirected_phenotype_graph, edge) == false continue end

	    node = src(edge)    # we ignore sinks
        if node in already_assigned_nodes continue end

	    neutral_network = _find_neutral_network(undirected_phenotype_graph, node, epsilon)
        push!(neutral_networks, neutral_network)

        append!(already_assigned_nodes, neutral_network)

        if i % 1000 == 0
            @info " - $i, $(length(neutral_networks)) NNs"
        end
    end

    graph_data.neutral_networks = neutral_networks
end

function _find_neutral_network(graph, starting_vertex, epsilon)
    neutral_network = [starting_vertex]

    to_test = [(starting_vertex, n) for n in outneighbors(graph, starting_vertex)]

    while 0 < length(to_test)
        u, v = pop!(to_test)

        if graph.weights[v, u] <= epsilon
            push!(neutral_network, v)

            for n in outneighbors(graph, v)
                if n in neutral_network continue end

                push!(to_test, (v, n))
            end
        end
    end

    return unique(neutral_network)
end

function get_neutral_networks_by_g_sampling(graph_data, epsilon, n, tolerance, min_samples, max_samples)
    num_genotypes = nv(graph_data.genotype_graph)

    neutral_networks = []

    for _ in 1:n
        starting_genotype = rand(1:num_genotypes)
        neutral_network = build_neutral_network(graph_data, starting_genotype, epsilon, tolerance, min_samples, max_samples)
        push!(neutral_networks, neutral_network)
        @info "$(length(neutral_networks))-th generated with $(length(neutral_network)) nodes"
    end

    return neutral_networks
end

function build_neutral_network(graph_data, starting_vertex, epsilon, tolerance, min_samples, max_samples)
    neutral_network = []
    _build_neutral_network(graph_data, starting_vertex, neutral_network, epsilon, tolerance, min_samples, max_samples)
    return neutral_network
end

function _build_neutral_network(graph_data, vertex, already_built_nn, epsilon, tolerance, min_samples, max_samples)
    push!(already_built_nn, vertex)

    neighbors = outneighbors(graph_data.genotype_graph, vertex)

    for n in neighbors
        if n in already_built_nn continue end

        if has_edge(graph_data.phenotype_graph, vertex, n)
            phenotype_similarity = graph_data.phenotype_graph.weights[n, vertex]
        else
            genotype_u, genotype_v = graph_data.genotype_vertex_mapping(vertex), graph_data.genotype_vertex_mapping(n)

            phenotype_similarity = get_phenotype_similarity(genotype_u, genotype_v, tolerance, min_samples, max_samples)
            phenotype_similarity = max(eps(), phenotype_similarity)

            add_edge!(graph_data.phenotype_graph, vertex, n, phenotype_similarity)
        end

        if phenotype_similarity <= epsilon
            _build_neutral_network(graph_data, n, already_built_nn, epsilon, tolerance, min_samples, max_samples)
        end
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module Sand_01
using LinearAlgebra, Random
#using StaticArrays
"""
    Sand 01: simple uniform advection of a scalar
"""

mutable struct ring
    N::Integer              # number of cells
    Δt::Float64             # time step
    X_::Array{Float64,1}    # positions (nodes)
    ΔX_::Array{Float64,1}   # cell sizes (cell)
    U_::Array{Float64,1}    # velocity (on nodes)
    Σ_::Array{Float64,1}    # surface density
    S_::Array{Float64,1}    # integral of the above
    t::Float64              # time
end
function ring(N=100; U = 1.0, xi=0.0, xf=10.0)
    xtemp_ = collect(range(xi,xf,length=N+1))
    xtemp_[2:end-1] += (Random.rand(N-1).-0.5) * (xf-xi)/N * 0.9 # 0.8 is somewhat arbitrary
    @assert prod(diff(xtemp_).>0)  # checks to make sure xtemp_ monotonic
    X_ = xtemp_#[2:end]
    ΔX_ = diff(xtemp_)
    U_ = U * ones(N+1)
    """ U_[3*N÷4:end] .*= 0.5 # <--- this does NOT work properly; next step is to fix this. """
    # init w/top-hat density:
    Σ_ = ones(N); Σ_[1] = 0.0;
    Σ_[1:(2*N÷5)] .= 0.0; Σ_[(3*N÷5):end] .= 0.0
    S_ = Σ_ .* ΔX_
    #
    Δt = 0.5 / max(maximum(U_[1:end-1]./ΔX_),-minimum(U_[2:end]./ΔX_))
    #
    return ring(N,Δt,X_,ΔX_,U_,Σ_,S_,0.0)
end

function step!(r::ring; Δt=r.Δt)
    Σ_=r.Σ_; N = r.N; ΔX_ = r.ΔX_; U_ = r.U_; S_ = r.S_
    #
    ΔΣ_ = cat(Σ_[1]-Σ_[end], diff(Σ_), Σ_[1]-Σ_[end], dims=1)
    S_ = Σ_ .* ΔX_
    U⁺_ = U_[1:end-1] .* (U_[1:end-1].>0) # right-traveling
    U⁻_ = U_[2:end  ] .* (U_[2:end  ].<0) # left-traveling
    S_ -= Δt * U⁺_ .* ΔΣ_[1:end-1]
    S_ -= Δt * U⁻_ .* ΔΣ_[2:end]
    #
    r.S_ = S_
    r.Σ_ = S_ ./ ΔX_
    r.t += Δt
#    print(sum(S_),"\n")
    return sum(S_)
end

function tstep!(r::ring,Δt)
    M = floor(Δt / r.Δt)
    Δt_last = Δt - M*r.Δt
    @assert Δt_last ≥ 0
    for i in 1:M
        step!(r)
    end
    step!(r, Δt = Δt_last)
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
export IsMeasure

@trait IsMeasure{M,X} where {X = eltype(M)}

###################################################

export HasDensity
export logdensity

@trait IsMeasure{M, X} >: HasDensity{M, X} where {X = eltype(M)} begin
    logdensity :: [M, X] => Real
end

###################################################

export IsMeasurable
export measure

@trait IsMeasure{M,X} >: IsMeasurable{M,S,X} where {X = eltype(S)} begin
    measure :: [M, S] => Real
end

###################################################

export HasRand

@trait IsMeasure{M,X} >: HasRand{M,X} where {X = eltype(M)} begin
    rand :: [M] => eltype(M)
end

###################################################

# @trait Measure{M,X} >: FiniteMeasure{M,X} where {X = eltype{M}} begin

# end


text:
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "lateral_driving_models",  "lateral_driving_models.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "longitudinal_driving_models",  "longitudinal_driving_models.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "lane_change_models",  "lane_change_models.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "lane_change_models",  "MOBIL.jl"))

include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "static_gaussian_drivers.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "prerecorded_drivers.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "lat_lon_separable_drivers.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "lane_following_drivers.jl"))
include(Pkg.dir("AutomotiveDrivingModels", "src", "behaviors", "tim2d_drivers.jl"))