text:
#using Pkg
#Pkg.add(["DelimitedFiles", "ProgressMeter", "StatsBase", "Roots", "LinearAlgebra", "Printf", "HDF5", "Distributions", "SpecialFunctions", "Parameters"])
#Pkg.add(url="https://github.com/JuliaIO/HDF5.jl")
#
using HDF5
using Printf
using Statistics
using StatsBase
using LinearAlgebra
using Roots
using IterTools

using ProgressMeter

include("utils.jl")
include("load_data.jl")
include("write_files.jl")
include("temporal_features.jl")
include("disease_spread.jl")
include("error.jl")
include("sampling.jl")
include("surrogate.jl")

"""
# when using only intermediate functions:
    cas, ets, list_contacts = load_processed_data()
    filename="./test.h5"
    experiment=Copenhagen()
    skip_jackknife=true
"""
function analyse_all(
        experiment::ContactData;
        # optional
        minimum_duration = 15*60,
        path_dat = "./dat",
        path_out = "./out",
        filter_out_incomplete=false,
        ct_range_times = 0:60*60:10*24*60*60,
        cluster_minimum_separation = 15*60,
        disease_model = DiseaseDefault(),
        seed=1000,
        # 0: essentials, 1: essentials and fast extras, 2: most, 3: everything
        # to skip jackknife analysis set gloibal variable `skip_jackknife = true`
        level_of_details = 2,
    )
    lod = level_of_details
    if filter_out_incomplete
        filename = @sprintf("%s/results_%s_filtered_%dmin_filterOutIncomplete.h5", path_out, label(experiment), minimum_duration/60)
    else
        filename = @sprintf("%s/results_%s_filtered_%dmin.h5", path_out, label(experiment), minimum_duration/60)
    end

    ###########################################################################
    # data
    # experiment=Copenhagen(); minimum_duration=15*60; path_dat="./dat";filter_out_incomplete=false;
    cas, ets, list_contacts = load_processed_data(experiment, minimum_duration, path_dat, filter_out_incomplete=filter_out_incomplete);
    # write encounter train
    lod >= 1 && myh5write(filename,"/data/trains/",ets)

    if experiment == InVS15()
        support_crate=0:300:seconds_from_days(7*1.5)
    else
        support_crate=0:timestep(ets):seconds_from_days(7*1.5)
    end

    ###########################################################################
    # temporal features

    lod >= 2 && analyse_temporal_features_of_contact_activity(cas, filename, "/data/contact_activity")
    lod >= 1 && analyse_contact_duration(list_contacts, experiment, filename, "/data/contacts")
    lod >= 1 && analyse_temporal_features_of_encounter_train(ets, filename, "/data/encounter_train", support_crate=support_crate)

    ###########################################################################
    # surrogate data and sampling

    # disease periods
    range_latent = 0:0.05:8
    range_infectious = 0.05:0.05:8

    if lod >= 1
        root="/data_surrogate_randomize_per_train"
        sur_rand_train = surrogate_randomize_per_train(ets, seed)
        ref = sur_rand_train
        myh5write(filename,@sprintf("%s/trains/",root), sur_rand_train)
        analyse_temporal_features_of_encounter_train(sur_rand_train, filename, @sprintf("%s/encounter_train",root), support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, sur_rand_train, ref, filename, @sprintf("%s/disease/delta", root))
    end

    if lod >= 2
        # surrogate version 2
        root="/data_surrogate_randomize_all"
        sur_rand_all = surrogate_randomize_all(ets, seed)
        myh5write(filename,@sprintf("%s/trains/",root), sur_rand_all)
        analyse_temporal_features_of_encounter_train(sur_rand_all, filename, @sprintf("%s/encounter_train",root), support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, sur_rand_all, ref, filename, @sprintf("%s/disease/delta", root))
    end

    # surrogates from sampling
    if lod >= 3
        num_sample_trains = length(ets)
        # weights for trains
        train_weights = length.(ets) ./ mean(length.(ets))

        # (inhomogeneous) poisson processes
        contact_rate = rate(ets, 0:timestep(ets):seconds_from_days(7))
        mean_rate = mean(contact_rate.weights)

        root="sample/poisson_homogeneous/"
        samples = sample_encounter_trains_poisson(mean_rate, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets))
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))

        root="sample/poisson_homogeneous_weighted_trains/"
        samples = sample_encounter_trains_poisson(mean_rate, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets), weights=train_weights)
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))

        root="sample/poisson_inhomogeneous/"
        samples = sample_encounter_trains_poisson(contact_rate, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets))
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))

        root="sample/poisson_inhomogeneous_weighted_trains/"
        samples = sample_encounter_trains_poisson(contact_rate, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets), weights=train_weights)
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))

        # weibull renewal processes
        list_dts = inter_encounter_intervals(ets);
        list_dts = list_dts[length.(list_dts) .> 0];
        dist = distribution_durations(list_dts, timestep=timestep(ets));
        args_weibull = fit_Weibull(dist)

        # weibull renewal
        root="sample/weibul_renewal_process/"
        samples = sample_encounter_trains_weibull_renewal(args_weibull, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets))
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))

        # weibull renewal with weighted trains
        root="sample/weibul_renewal_process_weighted_trains/"
        samples = sample_encounter_trains_weibull_renewal(args_weibull, num_sample_trains, -seconds_from_days(7), (0, duration(ets)), MersenneTwister(seed), timestep=timestep(ets), weights=train_weights)
        analyse_temporal_features_of_encounter_train(samples, filename, root, support_crate=support_crate)
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, samples, ref, filename, @sprintf("%s/disease/delta", root))
    end

    ###########################################################################
    # disease spread
    if lod >= 1
        analyse_infectious_encounter_scan_delta(range_latent, range_infectious, ets, ref, filename, "/disease/delta")
        #details
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(2), seconds_from_days(3)), ets, filename, "/disease/delta_2_3")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(6), seconds_from_days(3)), ets, filename, "/disease/delta_6_3")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(2), seconds_from_days(3)), ref, filename, "/disease/delta_2_3_surrogate")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(6), seconds_from_days(3)), ref, filename, "/disease/delta_6_3_surrogate")
        # details at different points
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5)), ets, filename, "/disease/delta_1_0.5")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5)), ets, filename, "/disease/delta_1.5_0.5")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5)), ref, filename, "/disease/delta_1_0.5_surrogate")
        analyse_infectious_encounter_detail(DeltaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5)), ref, filename, "/disease/delta_1.5_0.5_surrogate")
    end

    #gamma - scan
    if lod >= 1
        range_k = 10 .^ collect(0:0.2:5)
        analyse_infectious_encounter_scan_gamma((2,3), range_k, ets, filename, "/disease/gamma_2_3")
        analyse_infectious_encounter_scan_gamma((6,3), range_k, ets, filename, "/disease/gamma_6_3")
        analyse_infectious_encounter_scan_gamma((2,3), range_k, ref, filename, "/disease/gamma_2_3_surrogate")
        analyse_infectious_encounter_scan_gamma((6,3), range_k, ref, filename, "/disease/gamma_6_3_surrogate")
        #gamma - detail
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(2), seconds_from_days(3), 1), ets, filename, "/disease/gamma_2_3/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(6), seconds_from_days(3), 1), ets, filename, "/disease/gamma_6_3/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(2), seconds_from_days(3), 10), ets, filename, "/disease/gamma_2_3/k_10.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(6), seconds_from_days(3), 10), ets, filename, "/disease/gamma_6_3/k_10.0")
        # gamma - detail, surrogate
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(2), seconds_from_days(3), 1), ref, filename, "/disease/gamma_2_3_surrogate/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(6), seconds_from_days(3), 1), ref, filename, "/disease/gamma_6_3_surrogate/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(2), seconds_from_days(3), 10), ref, filename, "/disease/gamma_2_3_surrogate/k_10.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(6), seconds_from_days(3), 10), ref, filename, "/disease/gamma_6_3_surrogate/k_10.0")
    end

    #gamma - scan, more unrealistic points
    if lod >= 1
        range_k = 10 .^ collect(0:0.2:5)
        analyse_infectious_encounter_scan_gamma((1,0.5), range_k, ets, filename, "/disease/gamma_1_0.5")
        analyse_infectious_encounter_scan_gamma((1.5,0.5), range_k, ets, filename, "/disease/gamma_1.5_0.5")
        analyse_infectious_encounter_scan_gamma((1,0.5), range_k, ref, filename, "/disease/gamma_1_0.5_surrogate")
        analyse_infectious_encounter_scan_gamma((1.5,0.5), range_k, ref, filename, "/disease/gamma_1.5_0.5_surrogate")
        #gamma - detail
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5), 1), ets, filename, "/disease/gamma_1_0.5/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5), 1), ets, filename, "/disease/gamma_1.5_0.5/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5), 10), ets, filename, "/disease/gamma_1_0.5/k_10.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5), 10), ets, filename, "/disease/gamma_1.5_0.5/k_10.0")
        # gamma - detail, surrogate
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5), 1), ref, filename, "/disease/gamma_1_0.5_surrogate/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5), 1), ref, filename, "/disease/gamma_1.5_0.5_surrogate/k_1.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1), seconds_from_days(0.5), 10), ref, filename, "/disease/gamma_1_0.5_surrogate/k_10.0")
        analyse_infectious_encounter_detail(GammaDiseaseModel(seconds_from_days(1.5), seconds_from_days(0.5), 10), ref, filename, "/disease/gamma_1.5_0.5_surrogate/k_10.0")
    end

    return
end

###############################################################################
###############################################################################
### analyse contact activitity
function analyse_temporal_features_of_contact_activity(
        cas::contact_activities{I,T},
        filename::String,
        root::String;
        ############ optional
        lags=0:1:Int(seconds_from_days(7)/step(times(cas))),    # 1-week autocorrelation function
        times=0:60*60:seconds_from_days(7)           # 1-week contact rate (taken directly from first argument, seems to work)

    ) where {I,T}
    println("writing to ", root)
    println("... coefficient of variation")
    cv, cvj, cverr = coefficient_variation(cas, jackknife_error=true)
    myh5write(filename, @sprintf("%s/coefficient_variation", root), hcat(cv, cvj, cverr))
    myh5desc(filename, @sprintf("%s/coefficient_variation", root),
             "Coefficient of variation (std/mean of contact per timestep), d1: cv(full), d2: cv(jackknife), d3: error")

    #println("... autocorrelation function")
    #c, cj, cerr = jackknife(x->autocorrelation_function(x, lags), cas)
    #myh5write(filename, @sprintf("%s/autocorrelation", root), hcat(lags, c, cj, cerr))
    #myh5desc(filename, @sprintf("%s/autocorrelation", root),
    #         "autocorrelation of individual contacts over time, d1: timetep-lag, d2: C(full), d3: C(jackknife), d4: error")

    println("... contact rate")
    r, rj, rerr = jackknife(x->rate(x, times).weights, cas)
    myh5write(filename, @sprintf("%s/rate", root), hcat(times[1:end-1], r, rj, rerr))
    myh5desc(filename, @sprintf("%s/rate", root),
             "contact rate averaged across indivdials and across repetition of `times` in experiment, d1: times, d2: rate(full), d3: rate(jackknife), d4: error")

    println("... contact rate (high precision)")
    times_highres = 0:300:seconds_from_days(7)
    r, rj, rerr = jackknife(x->rate(x, times_highres).weights, cas)
    myh5write(filename, @sprintf("%s/rate_5min", root), hcat(times_highres[1:end-1], r, rj, rerr))
    myh5desc(filename, @sprintf("%s/rate_5min", root),
             "same as rate but at 5 min timesteps")

    println("... distribution of cluster sizes")
    samples_cluster_sizes = sample_cluster_sizes(cas)
    samples_cluster_sizes = samples_cluster_sizes[length.(samples_cluster_sizes) .> 0]
    full_dist = distribution(samples_cluster_sizes);
    weights_cluster(x) = distribution(x, edges=full_dist.edges[1]).weights
    P, Pj, Perr = jackknife_mj(weights_cluster, samples_cluster_sizes, naive=full_dist.weights);
    myh5write(filename, @sprintf("%s/distribution_encounter_cluster", root), hcat(full_dist.edges[1][1:end-1], P, Pj, Perr))
    myh5desc(filename, @sprintf("%s/distribution_encounter_cluster", root),
             "distribution of contact cluster (total sum of successive activity separated by zero activity), d1: size, d2: P(full), d3: P(jackknife), d4: error")

    xbin, Pbin, = logbin(full_dist)
    array_cluster_logbin(x) = hcat(logbin(distribution(x, edges=full_dist.edges[1]))...)[:,1:2]
    xP, xP_J, xP_err = jackknife_mj(array_cluster_logbin, samples_cluster_sizes, naive=hcat(xbin,Pbin));
    myh5write(filename, @sprintf("%s/distribution_encounter_cluster_logbin", root), hcat(xP, xP_J, xP_err))
    myh5desc(filename, @sprintf("%s/distribution_encounter_cluster_logbin", root),
             "distribution of inter-encounter intervals with log binning, d1: <duration>, d2: <P>, d3: <duration>(jackknife), d4: <P>(jackknife), d5: err(<duration>), d6: err(<P>)")
end

###############################################################################
###############################################################################
### analyse contact stuff
function  analyse_contact_duration(list_contacts, experiment::ContactData, filename::String, root::String)
    println("writing to ", root)
    list_durations = [getindex.(contacts, 2) for contacts in list_contacts];
    # remove lists with zero elements because they 1) do not contribute to the
    # statistics of the durations and 2) thereby break the jackknife routine
    list_durations = list_durations[length.(list_durations) .> 0]

    println("... distribution of contact durations")
    # full distribution
    full_dist = distribution_durations(list_durations, timestep=timestep(experiment));
    weights(x) = distribution_durations(x, edges=full_dist.edges[1]).weights
    P, Pj, Perr = jackknife_mj(weights, list_durations, naive=full_dist.weights);
    myh5write(filename, @sprintf("%s/distribution_contact_durations", root), hcat(full_dist.edges[1][1:end-1], P, Pj, Perr))
    myh5desc(filename, @sprintf("%s/distribution_contact_durations", root),
             "distribution of contact durations, d1: duration, d2: P(full), d3: P(jackknife), d4: error")

    # log-binning: remove initial range of short durations which are trivially excluded by our filter
    mask = trues(length(full_dist.edges[1]))
    mask[1:findfirst(full_dist.weights .> 0)] .= false
    full_dist = distribution_durations(list_durations, edges=full_dist.edges[1][mask])
    xbin, Pbin = logbin(full_dist)
    array_xP_logbin(x) = hcat(logbin(distribution_durations(x, edges=full_dist.edges[1]))...)[:,1:2]
    xP, xP_J, xP_err = jackknife_mj(array_xP_logbin, list_durations, naive=hcat(xbin,Pbin));
    myh5write(filename, @sprintf("%s/distribution_contact_durations_logbin", root), hcat(xP, xP_J, xP_err))
    myh5desc(filename, @sprintf("%s/distribution_contact_durations_logbin", root),
             "distribution of contact durations with log binning, d1: <duration>, d2: <P>, d3: <duration>(jackknife), d4: <P>(jackknife), d5: err(<duration>), d6: err(<P>)")
end

###############################################################################
###############################################################################
### analyse contact trains
function analyse_temporal_features_of_encounter_train(
        ets::encounter_trains{I,T},
        filename::String,
        root::String;
        ############ optional
        lags=0:1:200,                                         # lags for autocorrelation
        times_rate=0:60*60:seconds_from_days(7),        # 1-week contact rate in steps of hours
        support_crate=0:timestep(ets):seconds_from_days(7*1.5),# 1.5-week conditional contact rate

    ) where {I,T}
    println("writing to ", root)

    println("... distribution total number encounter per train")
    edges=0:1:700
    f_weights(x) = normalize!(fit(Histogram{Float64}, length.(trains(x)), edges)).weights
    w, wj, werr = jackknife(f_weights, ets)
    myh5write(filename, @sprintf("%s/distribution_total_number_encounter", root), hcat(edges[1:end-1], w, wj, werr))
    myh5desc(filename, @sprintf("%s/distribution_total_number_encounter", root),
             "distribution of total number of encounter across trains, d1: num encounter, d2: P(full), d3: P(jack), d4: error")
    fit_exp = fit(Exponential, length.(trains(ets)))
    myh5write(filename, @sprintf("%s/distribution_total_number_encounter_fit_exp", root), [mean(fit_exp)])
    myh5desc(filename, @sprintf("%s/distribution_total_number_encounter_fit_exp", root),
             "fit result for exponential distribution: scale (mean)")
    xbin, Pbin = logbin(edges[1:end-1], w, increment_factor=1, bin_width_start=10)
    # linear binning due to increment factor 1
    array_xP_linbin(x) = hcat(logbin(edges[1:end-1], f_weights(x), increment_factor=1, bin_width_start=10)...)[:,1:2]
    xP, xP_J, xP_err = jackknife(array_xP_linbin, ets, naive=hcat(xbin,Pbin));
    myh5write(filename, @sprintf("%s/distribution_total_number_encounter_linbin", root), hcat(xP, xP_J, xP_err))
    myh5desc(filename, @sprintf("%s/distribution_total_number_encounter_linbin", root),
             "distribution of inter-encounter intervals with log binning, d1: <number>, d2: <P>, d3: <number>(jackknife), d4: <P>(jackknife), d5: err(<number>), d6: err(<P>)")

    println("... daily number of encounter")
    # error is simple jackknife because the basic statistics is the number of
    # encounter per time window, which is the same across all trains (result is
    # a clear zero for empty trains)
    dist_daily = distribution_number_encounter(ets, seconds_from_days(1))[1]
    w, wj, werr = jackknife(x->distribution_number_encounter(x, seconds_from_days(1), edges=dist_daily.edges[1])[1].weights, ets);
    myh5write(filename, @sprintf("%s/distribution_daily_number_encounters", root), hcat(dist_daily.edges[1][1:end-1], w, wj, werr))
    myh5desc(filename, @sprintf("%s/distribution_daily_number_encounters", root),
             "distribution of daily number of encounter (weekend AND weekday), d1: num encounter, d2: P(full) , d3: P(jack), d4: error")

    # partition experiment into windows of size 1 day and sort according to
    # pattern to groups 1 (weekend) and 2 (weekday)
    pattern = [1,2*ones(5)...,1]
    # we want the results of distribution_number_contacts(ets, seconds_from_days(1), pattern=pattern, edges=dist_daily.edges[1])
    # but this is a dictionary of distributions. Instead, for the jackknife
    # error function we need a 2D Array to compute all errors consistently with
    # a single run -> hence we transform to an array of log-weights
    dict_dist(x, edges) = distribution_number_encounter(x, seconds_from_days(1), pattern, edges=edges)
    array_weights(x, edges) = hcat([dist.weights for dist in values(dict_dist(x, edges))]...)
    w, wj, werr = jackknife(x->array_weights(x, dist_daily.edges[1]), ets);
    myh5write(filename, @sprintf("%s/distribution_daily_number_encounters_sep", root), hcat(dist_daily.edges[1][1:end-1], w, wj, werr))
    myh5desc(filename, @sprintf("%s/distribution_daily_number_encounters_sep", root),
             "distribution of daily number of encounter (weekend VS weekday), d1: num encounter, d2: P(full) weekend, d3: P(full) weekday, d4: P(jack) weekend, d5: P(jack) weekday, d6: err weekend, d7: err weekday")

    println("... encounter rate")
    # error is simple jackknife because the basic statistics is number of
    # coontacts per time step, which is the same across all trains
    r, rj, rerr = jackknife(x->rate(x, times_rate).weights, ets)
    myh5write(filename, @sprintf("%s/rate", root), hcat(times_rate[1:end-1], r, rj, rerr))
    myh5desc(filename, @sprintf("%s/rate", root),
             "encounter rate averaged across indivdials and across repetition of `times` in experiment, d1: times, d2: rate(full), d3: rate(jackknife), d4: error")

    # we also need to write the encounter rate at the higher time-resolution because
    # the convolution with the duration needs to match time discretization
    times_rate_highres = 0:300:seconds_from_days(7)
    r, rj, rerr = jackknife(x->rate(x, times_rate_highres).weights, ets)
    myh5write(filename, @sprintf("%s/rate_5min", root), hcat(times_rate_highres[1:end-1], r, rj, rerr))
    myh5desc(filename, @sprintf("%s/rate_5min", root),
             "same as `rate` but at 5 min timesteps")

    println("... conditional encounter rate")
    # error is delete-mj jackknife because basic statistics is the average over encounter that we condition on. Hence, we need to remove the emtpy trains
    ets_noempty = ets[length.(ets).>0]
    r, rj, rerr = jackknife_mj(x->conditional_encounter_rate(x, support_crate).weights, ets_noempty)
    myh5write(filename, @sprintf("%s/conditional_encounter_rate", root), hcat(support_crate[1:end-1], r, rj, rerr))
    myh5desc(filename, @sprintf("%s/conditional_encounter_rate", root),
             "average encounter rate conditioned on having an encounter at time 0, averaged across all encounter in experiment, d1: times, d2: rate(full), d3: rate(jackknife), d4: error")

    println("... distribution of inter-encounter intervals")
    # error is delete-mj jackknife because basic statistics is the
    # inter-encounter interval, a random variable that occurs differenlty often
    # in different trains
    list_dts = inter_encounter_intervals(ets);
    list_dts = list_dts[length.(list_dts) .> 0]
    full_dist = distribution_durations(list_dts, timestep=timestep(ets));
    weights(x) = distribution_durations(x, edges=full_dist.edges[1]).weights
    P, Pj, Perr = jackknife_mj(weights, list_dts, naive=full_dist.weights);
    myh5write(filename, @sprintf("%s/distribution_inter_encounter_intervals", root), hcat(full_dist.edges[1][1:end-1], P, Pj, Perr))
    myh5desc(filename, @sprintf("%s/distribution_inter_encounter_intervals", root),
             "distribution of inter-encounter intervals, d1: duration, d2: P(full), d3: P(jackknife), d4: error")

    xbin, Pbin, = logbin(full_dist)
    array_xP_logbin(x) = hcat(logbin(distribution_durations(x, edges=full_dist.edges[1]))...)[:,1:2]
    xP, xP_J, xP_err = jackknife_mj(array_xP_logbin, list_dts, naive=hcat(xbin,Pbin));
    myh5write(filename, @sprintf("%s/distribution_inter_encounter_intervals_logbin", root), hcat(xP, xP_J, xP_err))
    myh5desc(filename, @sprintf("%s/distribution_inter_encounter_intervals_logbin", root),
             "distribution of inter-encounter intervals with log binning, d1: <duration>, d2: <P>, d3: <duration>(jackknife), d4: <P>(jackknife), d5: err(<duration>), d6: err(<P>)")

    println("... autocorrelation function of inter-encounter interval")
    c, cj, cerr = jackknife_mj(x->autocorrelation_function(x, lags), list_dts)
    myh5write(filename, @sprintf("%s/autocorrelation_inter_encounter_intervals", root), hcat(lags, c, cj, cerr))
    myh5desc(filename, @sprintf("%s/autocorrelation_inter_encounter_intervals", root),
             "autocorrelation of inter-encounter intervals, d1: k, d2: C(full), d3: C(jackknife), d4: error")

    println("... distribution of cluster sizes")
    minimum_separation = 8*300 # corresponds to roughly mean duration (0.655277 h ~ 8*300 s)
    samples_cluster_sizes = sample_cluster_sizes(list_dts, minimum_separation)
    samples_cluster_sizes = samples_cluster_sizes[length.(samples_cluster_sizes) .> 0]
    full_dist = distribution(samples_cluster_sizes);
    weights_cluster(x) = distribution(x, edges=full_dist.edges[1]).weights
    P, Pj, Perr = jackknife_mj(weights_cluster, samples_cluster_sizes, naive=full_dist.weights);
    myh5write(filename, @sprintf("%s/distribution_encounter_cluster", root), hcat(full_dist.edges[1][1:end-1], P, Pj, Perr))
    myh5desc(filename, @sprintf("%s/distribution_encounter_cluster", root),
             "distribution of contact cluster (total sum of successive activity separated by zero activity), d1: size, d2: P(full), d3: P(jackknife), d4: error")

    xbin, Pbin, = logbin(full_dist)
    array_cluster_logbin(x) = hcat(logbin(distribution(x, edges=full_dist.edges[1]))...)[:,1:2]
    xP, xP_J, xP_err = jackknife_mj(array_cluster_logbin, samples_cluster_sizes, naive=hcat(xbin,Pbin));
    myh5write(filename, @sprintf("%s/distribution_encounter_cluster_logbin", root), hcat(xP, xP_J, xP_err))
    myh5desc(filename, @sprintf("%s/distribution_encounter_cluster_logbin", root),
             "distribution of inter-encounter intervals with log binning, d1: <duration>, d2: <P>, d3: <duration>(jackknife), d4: <P>(jackknife), d5: err(<duration>), d6: err(<P>)")
end

###############################################################################
###############################################################################
### disease spread

function analyse_infectious_encounter_detail(
        disease_model::AbstractDiseaseModel,
        ets::encounter_trains{I,T},
        filename::String,
        root::String;
       ) where {I,T}
    println("... detailed analysis of infectious encounter for model ", disease_model)
    ets_noempty = ets[length.(ets) .> 0]

    # viral load
    println("... ... disease progression")
    time_range = 0:60*60:seconds_from_days(21)
    load = viral_load(disease_model, time_range)
    myh5write(filename, @sprintf("%s/viral_load", root), load)
    myh5desc(filename, @sprintf("%s/viral_load", root),
        "probability to be infectious for current disease model, d1: time range, d2: probability")

    # distribution of infectious encounters
    println("... ... distribution infectious encounter")
    try
        data = samples_infectious_encounter(disease_model, ets)
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights

        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename, @sprintf("%s/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")

        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/mean_number_infectious_encounter", root), hcat(o,oj,oerr))
        myh5desc(filename, @sprintf("%s/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end

    # controls
    println("... ... controls")
    try
        println("... ... ... random disase onset")
        data = samples_infectious_encounter_random_onset(disease_model, ets)
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights
        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/control_random_disease_onset/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")
        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/control_random_disease_onset/mean_number_infectious_encounter", root), hcat(o,oj,oerr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end
    try
        println("... ... ... random disase onset (weighted train)")
        data = samples_infectious_encounter_random_onset(disease_model, ets, weighted_train=true)
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights
        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtrain/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtrain/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")
        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtrain/mean_number_infectious_encounter", root), hcat(o,oj,oerr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtrain/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end
    try
        println("... ... ... random disase onset (weighted time)")
        data = samples_infectious_encounter_random_onset(disease_model, ets, weighted_time=true)
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights
        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtime/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtime/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")
        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtime/mean_number_infectious_encounter", root), hcat(o,oj,oerr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtime/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end
    try
        println("... ... ... random disase onset (weighted train/time)")
        data = samples_infectious_encounter_random_onset(disease_model, ets, weighted_time=true, weighted_train=true)
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights
        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtrain_wtime/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtrain_wtime/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")
        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/control_random_disease_onset_wtrain_wtime/mean_number_infectious_encounter", root),  hcat(o,oj,oerr))
        myh5desc(filename,  @sprintf("%s/control_random_disease_onset_wtrain_wtime/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end
    try
        println("... ... ... disease inset from encounters of first week only")
        data = samples_infectious_encounter_with_infection_in_early_interval(disease_model, ets, seconds_from_days(7))
        data_noempty = data[length.(data) .> 0]
        dist = distribution_from_samples_infectious_encounter(data)
        # jackknife_mj because relevant statistics is the number of encounters where infection could start
        weights(x) = distribution_from_samples_infectious_encounter(x, edges=dist.edges[1]).weights
        w, wj, werr = jackknife_mj(weights, data_noempty, naive=dist.weights);
        myh5write(filename, @sprintf("%s/control_sample_first_week_only/distribution_infectious_encounter", root), hcat(dist.edges[1][1:end-1], w, wj, werr))
        myh5desc(filename,  @sprintf("%s/control_sample_first_week_only/distribution_infectious_encounter", root),
                 "distribution of number of encounters during infectious interval, d1: number of encounters, d2: P(full), d3: P(jack), d4:error")
        o,oj, oerr = jackknife_mj(x->mean_from_samples_infectious_encounter(x), data_noempty)
        myh5write(filename, @sprintf("%s/control_sample_first_week_only/mean_number_infectious_encounter", root),  hcat(o,oj,oerr))
        myh5desc(filename,  @sprintf("%s/control_sample_first_week_only/mean_number_infectious_encounter", root),
                 "mean number of infectious encounter , d1: full, d2: jack, d3: err")
    catch e
        println(e)
    end
end

function analyse_infectious_encounter_scan_delta(
        range_latent_in_days,
        range_infectious_in_days,
        ets::encounter_trains{I,T1},
        ref::encounter_trains{I,T2},
        filename::String,
        root::String;
        # optional
        seed=1000,
       ) where {I,T1, T2}
    print("... scan mean number infectious encounters for delta disease\n")

    # scan
    results  = Array{Float64,2}(undef,length(range_latent_in_days),length(range_infectious_in_days))
    relative = Array{Float64,2}(undef,length(range_latent_in_days),length(range_infectious_in_days))
    @showprogress 1 for (i, latent) in enumerate(range_latent_in_days)
        for (j, infectious) in enumerate(range_infectious_in_days)
            mean_ets = mean_infectious_encounter(DeltaDiseaseModel(seconds_from_days(latent), seconds_from_days(infectious)), ets)
            mean_ref = mean_infectious_encounter(DeltaDiseaseModel(seconds_from_days(latent), seconds_from_days(infectious)), ref)
            results[i,j] = mean_ets
            relative[i,j] = mean_ets/mean_ref
        end
    end

    myh5write(filename, @sprintf("%s/scan_mean_number_infectious_encounter/mean", root), results)
    myh5write(filename, @sprintf("%s/scan_mean_number_infectious_encounter/mean_relative_to_poisson", root), relative)
    myh5write(filename, @sprintf("%s/scan_mean_number_infectious_encounter/range_latent", root), collect(range_latent_in_days))
    myh5write(filename, @sprintf("%s/scan_mean_number_infectious_encounter/range_infectious", root), collect(range_infectious_in_days))
    myh5desc(filename, @sprintf("%s/scan_mean_number_infectious_encounter", root),
        "scan of the mean number of infectious encounters for different latent and infectious periods, d1: latent period, d2: infectious period")
end

function analyse_infectious_encounter_scan_gamma(
        params_in_days,
        range_k,
        ets::encounter_trains{I,T},
        filename::String,
        root::String;
        # optional
        samples = Int(1e5)
    ) where {I,T}
    println("... scan k-dependence for gamma disease model with ", params_in_days)
    latent=getindex(params_in_days, 1)
    infectious=getindex(params_in_days, 2)

    # scan
    o = Vector{Float64}(undef,length(range_k))
    oj = Vector{Float64}(undef,length(range_k))
    oerr  = Vector{Float64}(undef,length(range_k))
    p = Progress(length(range_k), 1, "k-values: ", offset=0)
    for (i, k) in enumerate(range_k)
        next!(p; showvalues = [(:i,i), (:k,k)])
        # next!(p)
        model = GammaDiseaseModel(seconds_from_days(latent), seconds_from_days(infectious), k)
        rng = MersenneTwister(1000)
        try
            #first need to generate data (not mean infectious encounter but all samples of infectious encounter)
            data = samples_infectious_encounter(model, ets)
            data_noempty = data[length.(data).>0]
            o[i], oj[i], oerr[i] = jackknife_mj(x->mean(vcat(x...)), data_noempty)
        catch
            o[i] = oj[i] = oerr[i] = NaN
        end
    end

    myh5write(filename, @sprintf("%s/scan_k", root), hcat(collect(range_k), o, oj, oerr))
    myh5desc(filename, @sprintf("%s/scan_k", root),
             "k dependence of the mean number of infectious encounters for gamma-distributied disease models (parameters in dsets), d1: k, d2: mean(full), d3: mean(jack), d4: mean(err)")
end

###############################################################################
###############################################################################
### sample disease spread with generative model that captures some features of data
### and analyse effective R from the time course of new cases.

function analyse_effective_R(filename::String;
    path_out = "./out",
    dsetname::String="measurements")
    infectious = 3
    h5open(filename, "r") do file
        encounter_rate = read(file[dsetname], "encounter_rate")
        cases = file[dsetname]["cases"]
        keys_latent = keys(cases)
        keys_latent = keys_latent[occursin.("latent", keys_latent)]
        list_latent = parse.(Float64, last.(split.(keys_latent, "_")))
        # go through hdf5 file
        R4_avg = zeros(length(keys_latent))
        R4_std = zeros(length(keys_latent))
        Rg_avg = zeros(length(keys_latent))
        Rg_std = zeros(length(keys_latent))
        R0_avg = zeros(length(keys_latent))
        R0_std = zeros(length(keys_latent))
        Tg = zeros(length(keys_latent))
        for (l, key_latent) in enumerate(keys_latent)
            println("\n", key_latent, " ", list_latent[l])
            #todo: string manipulation
            disease_model = DeltaDiseaseModel(seconds_from_days(list_latent[l]), seconds_from_days(infectious))

            # calculate generation time as the expectation value over infection times
            interval_infectious  = infectious_interval(disease_model)
            mask_infectious = interval_infectious[1] .< encounter_rate[:,1] .< interval_infectious[2];
            time_generation = sum(encounter_rate[mask_infectious, 1].*encounter_rate[mask_infectious,2])/sum(encounter_rate[mask_infectious, 2])
            time_generation /= seconds_from_days(1)
            keys_seed = keys(cases[key_latent])
            mean_R4 = zeros(length(keys_seed))
            mean_Rg = zeros(length(keys_seed))
            mean_R0 = zeros(Union{Float64, Missing}, length(keys_seed))
            @showprogress 1 for (s, seed) in enumerate(keys_seed)
                println("...", seed)
                measurement = read(file, @sprintf("%s/cases/%s/%s/", dsetname, key_latent, seed))
                measurement = measurement[:, 2]
                r0, stat = read(file, @sprintf("%s/R0/%s/%s/", dsetname, key_latent, seed))
                mean_R0[s] = r0

                # find left border where cases >= some threshold
                large_weights = findall(measurement[:].>1000)
                if length(large_weights) < 1
                    mean_Rg[s] = NaN
                    mean_R4[s] = NaN
                    mean_R0[s] = missing
                    print("no data qualified, skipping ", r0, "\n")
                    continue
                end

                # include inbetween points where fluctuations drive case numbers below thresh
                valid_range = large_weights[1]:large_weights[end]
                # average R estimage across 3 days
                a = 3
                # time shift new vs old, first use estimated generation time from data
                g = round(Int64, time_generation)
                println("average over: ", a, " range from, to: ", valid_range, " time shift: ", g)

                if length(valid_range) <= g+a
                    mean_Rg[s] = NaN
                else
                    # calculate R with local (`a`-long) sums for every time point
                    R_g = zeros(length(valid_range)-g-a)
                    for (i,t) in enumerate(valid_range[1]:valid_range[end]-g-a)
                        R_g[i] = sum(measurement[t+g : t+g+a]) / sum(measurement[t : t+a])
                    end
                    # and average across all time points
                    mean_Rg[s] = mean(R_g)
                end

                # repeat for conventional g
                g = 4
                if length(valid_range) <= g+a
                    mean_R4[s] = NaN
                else
                    R_4 = zeros(length(valid_range)-g-a)
                    for (i,t) in enumerate(valid_range[1]:valid_range[end]-g-a)
                        R_4[i] = sum(measurement[t+g : t+g+a]) / sum(measurement[t : t+a])
                    end
                    mean_R4[s] = mean(R_4)
                end

                println("R0, Rg, R4: ", mean_R0[s], " ", mean_Rg[s], " ", mean_R4[s])
            end
            # l is latent period
            # mean_and_std returns nan if any element is nan.
            R4_avg[l], R4_std[l] = mean_and_std(mean_R4)
            Rg_avg[l], Rg_std[l] = mean_and_std(mean_Rg)
            R0_avg[l], R0_std[l] = mean_and_std(skipmissing(mean_R0))
            Tg[l] = time_generation
        end

        # store results somewhere
        open(@sprintf("%s/analysis_mean-field_R_%s.dat", path_out, dsetname), "w") do io
            write(io, "#latent\t avg(R4)\t std(R4)\t avg(Rg)\t std(Rg)\t avg(R0)\t std(R0)\t g\n")
            writedlm(io, zip(list_latent, R4_avg, R4_std, Rg_avg, Rg_std, R0_avg, R0_std, Tg))
        end
    end
end

function sample_mean_field_for_effective_R(;
        #optional
        experiment=Copenhagen(),
        minimum_duration = 15*60,
        path_dat = "./dat",
        path_out = "./out",
    )

    function do_it(ets, filename, dsetname)
        # encounter rate
        encounter_rate = conditional_encounter_rate(ets, 0:timestep(ets):seconds_from_days(8+4))
        myh5write(filename, @sprintf("/%s/encounter_rate", dsetname),
            hcat(encounter_rate.edges[1][1:end-1], encounter_rate.weights))
        myh5desc(filename, @sprintf("/%s/encounter_rate", dsetname),
            "average encounter rate conditioned on having an encounter at time 0, averaged across all encounter in experiment, d1: times, d2: rate(full)")

        probability_infection = 0.12 # heuristically chosen to match R approx 3.3 for infectious=3 and latent=4 when g=4
        infectious = 3 # days
        range_latent = 0.0:0.5:8
        for (l, latent) in enumerate(range_latent)
            println(latent)
            disease_model = DeltaDiseaseModel(seconds_from_days(latent), seconds_from_days(infectious))
            #dist_offspring = distribution_from_samples_infectious_encounter(samples_infectious_encounter(disease_model, ets))

            # sample
            rng = MersenneTwister(1000)
            range_seeds = 1000:1010
            @showprogress 1 for (s, seed) in enumerate(range_seeds)
                println("...", seed)
                I0 = floor(Int, 1 + latent)
                measurement, sum_offsprings, sum_samples = spread_mean_field(disease_model, encounter_rate, probability_infection, seconds_from_days.(0:120.), max_cases=5e6, seed=seed, initial_infection_times=rand(rng,I0)*infectious)

                myh5write(filename, @sprintf("/%s/cases/latent_%.2f/%d", dsetname, latent, seed), hcat(measurement.edges[1][1:end-1], measurement.weights))
                myh5write(filename, @sprintf("/%s/R0/latent_%.2f/%d", dsetname, latent, seed), [sum_offsprings/sum_samples, sum_samples])
            end
        end
    end

    # use real data
    filename = @sprintf("%s/mean_field_samples_%s_filtered_%dmin.h5", path_out, label(experiment), minimum_duration/60)
    cas, ets, list_contacts = load_processed_data(experiment, minimum_duration, path_dat);
    do_it(ets, filename, "measurements")

    # use randomized data
    ets = surrogate_randomize_per_train(ets, 1000)
    do_it(ets, filename, "measurements_randomized_per_train")
end

###############################################################################
###############################################################################
### quick default
analyse_all() = analyse_all(Copenhagen())


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

using Markdown
using InteractiveUtils

# ╔═╡ 5bd9e6aa-6820-11eb-2b20-3f17c0045841
begin
	using Test
	using BenchmarkTools
	using PlutoUI
end

# ╔═╡ 19046b20-6820-11eb-3fac-e993b8847dc0
md"""
### Sum of Intervals

Write a function called `sum_intervals()` that accepts an array of intervals, and returns the sum of all the interval lengths. Overlapping intervals should only be counted once.

**Intervals**:

Intervals are represented by a pair of integers in the form of a tuple. The first value of the interval will always be less than the second value. Interval example: `(1, 5)` is an interval from 1 to 5. The length of this interval is 4.

**Overlapping Intervals**:

List containing overlapping intervals:

```julia
[
  (1, 4),
  (7, 10),
  (3, 5)
]
```

The sum of the lengths of these intervals is 7. Since `(1, 4)` and `(3, 5)` overlap, we can treat the interval as `(1, 5)`, which has a length of 4.

**Examples**:

```julia
sum_intervals([
   (1, 2),
   (6, 10),
   (11, 15)
]); # => 9

sum_intervals([
   (1, 4),
   (7, 10),
   (3, 5)
]); # => 7

sum_intervals([
   (1, 5),
   (10, 20),
   (1, 6),
   (16, 19),
   (5, 11)
]); # => 19
```
"""

# ╔═╡ 698441e4-6820-11eb-2ed2-e90fbf53d7f0
md"""
#### Defining a parametric type
"""

# ╔═╡ 7da89592-6820-11eb-0ad7-699807d1d952
struct Interval{T <: Integer}
    lb::T  # lower bound
    ub::T  # upper bound

    function Interval{T}(lb, ub) where {T <: Integer}
		(lb, ub) = lb > ub ? (ub, lb) : (lb, ub)
		@assert 0 ≤ lb ≤ ub
		new(lb, ub)
	end
end

# ╔═╡ 9e07bbc4-6820-11eb-07ca-77f29e8a3030
IntervInt64 = Interval{Int64}

# ╔═╡ b9125d20-6820-11eb-3773-07dd4728988d
interv0 = IntervInt64(3, 1)

# ╔═╡ f6264b68-6820-11eb-1fc4-afcb22381c09
md"""
##### Length of an Interval
"""

# ╔═╡ 1640dbfe-6821-11eb-30bb-37b1f76609a3
function interv_len(interv::Interval{T})::T where T <: Integer
    interv.ub - interv.lb
end

# ╔═╡ 3cdd4d22-6821-11eb-3c20-a7997b70c44d
begin
	interv1 = IntervInt64(3, 1)
	@test interv_len(interv1) == 2

	interv2 = IntervInt64(4, 7)
	@test interv_len(interv2) == 3

	interv3 = IntervInt64(1, 1)
	@test interv_len(interv3) == 0

	interv4 = IntervInt64(5, 1)
	@test interv_len(interv4) == 4
end

# ╔═╡ 46c62dea-6821-11eb-0b18-771fd704ee84
md"""
##### Compute if two intervals overlap
"""

# ╔═╡ 91d63a66-6821-11eb-07a4-97c4a000b482
begin

function interv_inter(i₁::Interval{T}, i₂::Interval{T})::Bool where T <: Integer
	i₁ ⊆ i₂ || i₂ ⊆ i₁ || i₁ ⊆ i₂ || i₁ ∩ i₂ || i₂ ∩ i₁
end

function ⊆(i₁::Interval{T}, i₂::Interval{T})::Bool where T <: Integer
	lb1, ub1 = i₁.lb, i₁.ub
  	lb2, ub2 = i₂.lb, i₂.ub

	lb2 ≤ lb1 ≤ ub1 ≤ ub2
end

# \cap ≡ ∩
function ∩(i₁::Interval{T}, i₂::Interval{T})::Bool where T <: Integer
	lb1, ub1 = i₁.lb, i₁.ub
  	lb2, ub2 = i₂.lb, i₂.ub

	lb1 ≤ lb2 ≤ ub1 ≤ ub2
end

end

# ╔═╡ 14f1261a-6823-11eb-3d78-ebc8784d5ef8
begin
	@test interv_inter(IntervInt64(1, 10), IntervInt64(3, 7))
	@test interv_inter(IntervInt64(3, 7), IntervInt64(1, 8))
	@test interv_inter(IntervInt64(1, 6), IntervInt64(3, 9))
	@test interv_inter(IntervInt64(1, 5), IntervInt64(3, 4))
	@test interv_inter(IntervInt64(1, 5), IntervInt64(5, 6))
	@test interv_inter(IntervInt64(6, 10), IntervInt64(5, 8))
	@test !interv_inter(IntervInt64(1, 5), IntervInt64(6, 8))
	@test !interv_inter(IntervInt64(1, 5), IntervInt64(10, 20))
end

# ╔═╡ aac494f8-6823-11eb-138e-59b1641e416b
md"""
##### Adding two intervals length
"""

# ╔═╡ c2440f3a-6823-11eb-1fac-07057f45e3e3
function Base.:+(i₁::Interval{T}, i₂::Interval{T})::T where T <: Integer
	lb1, ub1 = i₁.lb, i₁.ub
	lb2, ub2 = i₂.lb, i₂.ub

	i₁ ⊆ i₂ && return ub2 - lb2
	i₂ ⊆ i₁ && return ub1 - lb1

   	s_inter = interv_len(i₁) + interv_len(i₂)

	# There might be some overlap: i₁.lb ---- i₂.lb ----- i₁.ub ---- i₂.ub
   	i₁ ∩ i₂ && return s_inter - (ub1 - lb2)
   	i₂ ∩ i₁ && return s_inter - (ub2 - lb1) # symmetry

   	s_inter
end

# ╔═╡ a68980ee-6824-11eb-06e2-b5a9527b24b5
begin
	@test IntervInt64(1, 10) + IntervInt64(3, 7) == 9
	@test IntervInt64(3, 7) + IntervInt64(1, 8) == 7
	@test IntervInt64(1, 6) + IntervInt64(3, 9) == 8
	@test IntervInt64(6, 10) + IntervInt64(5, 8) == 5
	@test IntervInt64(1, 5) + IntervInt64(6, 8) == 6
end

# ╔═╡ dda42a42-6825-11eb-1896-c7c0c2f6577d
md"""
##### Adding interval (length) and a length
"""

# ╔═╡ eb0fe112-6825-11eb-1d1c-d5c2fc96e381
begin

function Base.:+(i::Interval{T}, len::T)::T where T <: Integer
   interv_len(i) + len
end

function Base.:+(len::T, i::Interval{T})::T where T <: Integer
   interv_len(i) + len
end

end

# ╔═╡ 2be6f644-6826-11eb-168f-31248b1f7b57
begin
	@test IntervInt64(1, 10) + 10 == 19
	@test 6 + IntervInt64(2, 10) == 14
end

# ╔═╡ 647479b4-6826-11eb-0e26-8ffb421b8aba
md"""
##### Sorting an array of Intervals

We will sort on ascending order given the lower bound.
"""

# ╔═╡ 73948ede-6826-11eb-3414-e184635bb6b4
begin
	ary = [
    	IntervInt64(2, 4), IntervInt64(8, 10), IntervInt64(1, 5),
    	IntervInt64(1, 3), IntervInt64(4, 8), IntervInt64(3, 7),
    	IntervInt64(5, 11)
	]

	with_terminal() do
		sort!(ary, by = interv -> interv.lb, rev=false)
		for e ∈ ary
			println(e)
		end
	end
end

# ╔═╡ f204da94-6826-11eb-3007-c1a3d607e835
md"""
#### Merging overlapping intervals (in linear time)
"""

# ╔═╡ 30bf00fc-6827-11eb-14ea-8bf3ef119dfe
begin

const VT{T <: Integer} = Vector{Interval{T}}

function merge_interv(il::VT{T})::VT{T} where T <: Integer
	"""
	Merge all the intervals from given vector
	"""

	size(il, 1) == 1 && return il   ## identity

	n = size(il, 1)
	n_il = VT{T}(undef, n)  ## at most the size of given vector
	i₁, jx = il[1], 1

	for ix ∈ 2:n
		i₂ = il[ix]
		if interv_inter(i₁, i₂)
			i₁ = Interval{T}(min(i₁.lb, i₂.lb), max(i₁.ub, i₂.ub))
		else
			n_il[jx] = i₁
			jx += 1
			i₁ = i₂
		end
	end

	n_il[jx] = i₁
	view(n_il, 1:jx)
end

end

# ╔═╡ 344e9cda-6829-11eb-3c5d-61e4908b38a0
begin
	ary1 = [
    	IntervInt64(2, 4), IntervInt64(8, 10), IntervInt64(1, 5),
    	IntervInt64(1, 3), IntervInt64(4, 8), IntervInt64(3, 7),
    	IntervInt64(5, 11)
	]

	sort!(ary1, by = interv -> interv.lb)

	ary2 = [
    	IntervInt64(2, 4), IntervInt64(8, 10), IntervInt64(1, 5),
    	IntervInt64(17, 20), IntervInt64(11, 15), IntervInt64(3, 7),
	]

	sort!(ary2, by = interv -> interv.lb)

	with_terminal() do
		for e ∈ ary1
			println(e)
		end
		println()

		for e ∈ ary2
			println(e)
		end
	end
end

# ╔═╡ 88e56bca-6829-11eb-3da4-d9b2299f2801
@test merge_interv(ary1) == [IntervInt64(1, 11)]

# ╔═╡ fc900436-6829-11eb-1783-b1e82c817e8d
@test merge_interv(ary2) == [IntervInt64(1, 7),
	IntervInt64(8, 10),
    IntervInt64(11, 15),
	IntervInt64(17, 20)
]

# ╔═╡ 158a8cb8-682a-11eb-1d50-d922aa7274a2
begin
	ary3 = [IntervInt64(1, 2), IntervInt64(11, 15)]
	sort!(ary3, by = interv -> interv.lb)

	@test merge_interv(ary3) == ary3 # Identity
end

# ╔═╡ 27a6a1ac-682a-11eb-39f5-9b9159b03088
with_terminal() do
	@btime merge_interv(ary1)
end

# ╔═╡ 33bdfd96-682a-11eb-3f55-93f22f7b1fdf
with_terminal() do
	@btime merge_interv(ary2)
end

# ╔═╡ 8ff62b88-682a-11eb-355e-03ea9b57ab05
# import Base: convert #promote_rule

# overload + already done

Base.convert(::Type{T}, i::Interval{T}) where T <: Integer = i.ub - i.lb
# promote_rule(::Interval{T}, ::Interval{<:Integer}) = T

# ╔═╡ c7b6c8fc-682a-11eb-2cc9-5706878cf4a3
md"""
##### Compute Sum of Interval array

###### Take 1
"""

# ╔═╡ d359aff8-682a-11eb-1061-89a949359934
function sum_intervals(il::VT{T}; merger_fn=merge_interv)::T where T <: Integer
    size(il, 1) == 1 && return interv_len(il[1])

	foldl((s, i) -> s += interv_len(i),
		merger_fn(sort(il, by=interv -> interv.lb));
		init=0)
end

# ╔═╡ 697050a0-682b-11eb-0ec0-650a0a68b11a
begin
	s_ary11 = [
    	IntervInt64(2, 4), IntervInt64(8, 10), IntervInt64(1, 5),
    	IntervInt64(1, 3), IntervInt64(4, 8), IntervInt64(3, 7),
    	IntervInt64(5, 11)
	]

	s_ary12 = [
    	IntervInt64(2, 4), IntervInt64(8, 10), IntervInt64(1, 5),
    	IntervInt64(1, 3), IntervInt64(4, 8), IntervInt64(3, 7),
    	IntervInt64(5, 11)
	]

	@test sum_intervals(s_ary11) == sum_intervals(s_ary12)
end

# ╔═╡ 9cc8c4b4-682b-11eb-046f-b9f2e3e82908
begin
	_ary11 = [IntervInt64(1, 2), IntervInt64(11, 15)]
	@test sum_intervals(_ary11) == 5

	_ary12 = [IntervInt64(1, 2), IntervInt64(11, 15), IntervInt64(6, 10)]
	@test sum_intervals(_ary12) == 9

	_ary13 = [
    	IntervInt64(1, 5), IntervInt64(10, 20), IntervInt64(1, 6),
    	IntervInt64(16, 19), IntervInt64(5, 11)
	]
	@test sum_intervals(_ary13) == 19
end

# ╔═╡ ca5c153e-682b-11eb-390b-ef6bbc821d95
with_terminal() do
	@btime sum_intervals(_ary13)
end

# ╔═╡ 1ae7fe64-682c-11eb-349c-3f3919aafd76
md"""
###### Take 2
"""

# ╔═╡ 2d3db23e-682c-11eb-2baa-91d0f9caf485
function sum_intervals_2(il::VT{T}; merger_fn=merge_interv)::T where T <: Integer
 	size(il, 1) == 1 && return interv_len(il[1])

	map(iv -> interv_len(iv),
		merger_fn(sort(il, by=interv -> interv.lb))) |> sum

end

# ╔═╡ 21c48f2c-682c-11eb-2971-a34d13b77f1c
@test sum_intervals_2(s_ary11) == sum_intervals_2(s_ary12)

# ╔═╡ 9e3b563a-682c-11eb-03b6-292ba8b72db5
begin
	@test sum_intervals_2(_ary11) == 5
	@test sum_intervals_2(_ary12) == 9
	@test sum_intervals_2(_ary13) == 19
end

# ╔═╡ 134299a0-6852-11eb-17ad-f7fa7323406b
with_terminal() do
	@btime sum_intervals_2(_ary13)
end

# ╔═╡ 218593a8-682c-11eb-1eb8-579b44dfd57c
md"""
###### Take 3
"""

# ╔═╡ 04684c04-6852-11eb-3522-bf07f9428110
function sum_intervals_3(il::VT{T}; merger_fn=merge_interv)::T where T <: Integer
 	size(il, 1) == 1 && return interv_len(il[1])

    reduce(+,
		merger_fn(sort(il, by=interv -> interv.lb))
	)
end

# ╔═╡ 40cfa4dc-6853-11eb-2469-ef9c84d3a4f6
@test sum_intervals_3(s_ary11) == sum_intervals_3(s_ary12)

# ╔═╡ 39406316-6853-11eb-2783-35b77b3b62cc
begin
	@test sum_intervals_3(_ary11) == 5
	@test sum_intervals_3(_ary12) == 9
	@test sum_intervals_3(_ary13) == 19
end

# ╔═╡ 47aa8d46-6853-11eb-2505-5727d9c83eb2
with_terminal() do
	@btime sum_intervals_3(_ary13)
end

# ╔═╡ Cell order:
# ╟─19046b20-6820-11eb-3fac-e993b8847dc0
# ╠═5bd9e6aa-6820-11eb-2b20-3f17c0045841
# ╟─698441e4-6820-11eb-2ed2-e90fbf53d7f0
# ╠═7da89592-6820-11eb-0ad7-699807d1d952
# ╠═9e07bbc4-6820-11eb-07ca-77f29e8a3030
# ╠═b9125d20-6820-11eb-3773-07dd4728988d
# ╟─f6264b68-6820-11eb-1fc4-afcb22381c09
# ╠═1640dbfe-6821-11eb-30bb-37b1f76609a3
# ╠═3cdd4d22-6821-11eb-3c20-a7997b70c44d
# ╟─46c62dea-6821-11eb-0b18-771fd704ee84
# ╠═91d63a66-6821-11eb-07a4-97c4a000b482
# ╠═14f1261a-6823-11eb-3d78-ebc8784d5ef8
# ╟─aac494f8-6823-11eb-138e-59b1641e416b
# ╠═c2440f3a-6823-11eb-1fac-07057f45e3e3
# ╠═a68980ee-6824-11eb-06e2-b5a9527b24b5
# ╟─dda42a42-6825-11eb-1896-c7c0c2f6577d
# ╠═eb0fe112-6825-11eb-1d1c-d5c2fc96e381
# ╠═2be6f644-6826-11eb-168f-31248b1f7b57
# ╟─647479b4-6826-11eb-0e26-8ffb421b8aba
# ╠═73948ede-6826-11eb-3414-e184635bb6b4
# ╟─f204da94-6826-11eb-3007-c1a3d607e835
# ╠═30bf00fc-6827-11eb-14ea-8bf3ef119dfe
# ╠═344e9cda-6829-11eb-3c5d-61e4908b38a0
# ╠═88e56bca-6829-11eb-3da4-d9b2299f2801
# ╠═fc900436-6829-11eb-1783-b1e82c817e8d
# ╠═158a8cb8-682a-11eb-1d50-d922aa7274a2
# ╠═27a6a1ac-682a-11eb-39f5-9b9159b03088
# ╠═33bdfd96-682a-11eb-3f55-93f22f7b1fdf
# ╠═8ff62b88-682a-11eb-355e-03ea9b57ab05
# ╟─c7b6c8fc-682a-11eb-2cc9-5706878cf4a3
# ╠═d359aff8-682a-11eb-1061-89a949359934
# ╠═697050a0-682b-11eb-0ec0-650a0a68b11a
# ╠═9cc8c4b4-682b-11eb-046f-b9f2e3e82908
# ╠═ca5c153e-682b-11eb-390b-ef6bbc821d95
# ╟─1ae7fe64-682c-11eb-349c-3f3919aafd76
# ╠═2d3db23e-682c-11eb-2baa-91d0f9caf485
# ╠═21c48f2c-682c-11eb-2971-a34d13b77f1c
# ╠═9e3b563a-682c-11eb-03b6-292ba8b72db5
# ╠═134299a0-6852-11eb-17ad-f7fa7323406b
# ╟─218593a8-682c-11eb-1eb8-579b44dfd57c
# ╠═04684c04-6852-11eb-3522-bf07f9428110
# ╠═40cfa4dc-6853-11eb-2469-ef9c84d3a4f6
# ╠═39406316-6853-11eb-2783-35b77b3b62cc
# ╠═47aa8d46-6853-11eb-2505-5727d9c83eb2

ext: jl

max_stars_repo_path: Julia Kata/06_sum_of_intervals.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Julia Kata/06_sum_of_intervals.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Julia Kata/06_sum_of_intervals.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
@testset "966.vowel-spellchecker.jl" begin
    @test spellchecker(["KiTe", "kite", "hare", "Hare"],
                       ["kite", "Kite", "KiTe", "Hare", "HARE", "Hear", "hear", "keti",
                        "keet", "keto"]) ==
          ["kite", "KiTe", "KiTe", "Hare", "hare", "", "", "KiTe", "", "KiTe"]
end


text:
function monte_carlo_covariance(𝒞::Tuple{AbstractArray, Vararg{AbstractArray}}, 𝒟::Tuple{AbstractArray, Vararg{AbstractArray}},s,ntrial)
    ℳ, ℳʹ = collect(𝒟)
    Λ₁, Λ₂ = collect(𝒞)

    ℱ = zeros(9,ntrial)
    for itrial = 1:ntrial
      # Monte-carlo estimate of the covariance matrix.

      𝒪 = map(ℳ) do 𝐦
          Point2DH(𝑛(𝐦) + Point2DH(s*vcat(randn(2,1),0.0)))
      end
      𝒪ʹ = map(ℳʹ) do 𝐦ʹ
          Point2DH(𝑛(𝐦ʹ) + Point2DH(s*vcat(randn(2,1),0.0)))
      end

      𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(), (𝒪, 𝒪ʹ))
      𝐅 = estimate(FundamentalMatrix(),
                              FundamentalNumericalScheme(reshape(𝐅₀,9,1), 5, 1e-10),
                                                                (Λ₁,Λ₂), (𝒪, 𝒪ʹ))
      𝐟 = reshape(𝐅,9,1)
      𝐟 = 𝐟 / norm(𝐟)
      𝐟 = 𝐟 / sign(𝐟[2])
      ℱ[:,itrial] = 𝐟
      # for correspondence in zip(M, Mʹ)
      #     m , mʹ = correspondence
      #     # Add zero-mean Gaussian noise to the coordinates.
      #     𝐦  = 𝑛(collect(Float64,m.coords)) + s*vcat(rand(2,1),0.0)
      #     𝐦ʹ = 𝑛(collect(Float64,mʹ.coords)) + s*vcat(rand(2,1),0.0)
      #     𝒪[i] =  HomogeneousPoint(tuple(𝐦...))
      #     𝒪ʹ[i] =  HomogeneousPoint(tuple(𝐦ʹ...))
      #
      #     #@show 𝐦
      #     #@show 𝐦ʹ
      #     # Test the Fundamental Numerical Scheme on the Fundamental matrix problem.
      #     𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(), 𝒪, 𝒪ʹ)
      #     # 𝐅 = estimate(FundamentalMatrix(),
      #     #                         FundamentalNumericalScheme(reshape(𝐅₀,9,1), 5, 1e-10),
      #     #                         [eye(4) for i = 1:length(ℳ)],  𝒪, 𝒪ʹ)
      #     # 𝐟 = reshape(𝐅,9,1)
      #     # 𝐟 = 𝐟 / norm(𝐟)
      #     i = i + 1
      #     # 𝐟 = 𝐟 / sign(𝐟[2])
      # end
      #ℱ[:,itrial] = 𝐟
    end
    ℱ
end

#


text:
module TestOrchestraTransformers

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
fcp = FeatureClassification()
nfcp = NumericFeatureClassification()

using FactCheck

importall Orchestra.Transformers.OrchestraTransformers

facts("Orchestra transformers") do
  context("OneHotEncoder transforms nominal features") do
    instances = [
      2 "a" 1 "c";
      1 "b" 2 "d";
    ]
    labels = [
      "y";
      "z";
    ]
    expected_transformed = [
      2 1 0 1 1 0;
      1 0 1 2 0 1;
    ]
    encoder = OneHotEncoder()
    fit!(encoder, instances, labels)
    transformed = transform!(encoder, instances)

    @fact transformed => expected_transformed
  end
  context("OneHotEncoder transforms with options") do
    nominal_columns = [3, 5]
    nominal_column_values_map = {
      3 => ["a", "b", "c", "d", "e"],
      5 => ["a", "b", "c", "d", "e"]
    }
    encoder = OneHotEncoder({
      :nominal_columns => nominal_columns,
      :nominal_column_values_map => nominal_column_values_map
    })
    fit!(encoder, fcp.train_instances, fcp.train_labels)
    transformed = transform!(encoder, fcp.test_instances)

    @fact size(transformed, 2) => 13
  end
  context("OneHotEncoder handles no nominal features") do
    encoder = OneHotEncoder()
    fit!(encoder, nfcp.train_instances, nfcp.train_labels)
    transformed = transform!(encoder, nfcp.test_instances)

    @fact size(transformed, 2) => 2
  end

  context("Imputer replaces NA") do
    instances = [
      1.0      1.0;
      nan(1.0) 1.0;
      0.0      0.0;
      nan(0.0) 0.0;
    ]
    labels = ["x";"x";"y";"y"]
    expected_transformed = [
      1.0 1.0;
      0.5 1.0;
      0.0 0.0;
      0.5 0.0;
    ]
    imputer = Imputer()
    fit!(imputer, instances, labels)
    transformed = transform!(imputer, instances)

    @fact transformed => expected_transformed
  end

  context("Pipeline chains transformers") do
    pipe = Pipeline()
    fit!(pipe, fcp.train_instances, fcp.train_labels)
    transformed = transform!(pipe, fcp.test_instances)

    @fact size(transformed, 2) => 11
    @fact true => !any(map(x -> isnan(x), transformed))
  end

  context("Wrapper delegates to transformer") do
    wrapper = Wrapper({
      :transformer => OneHotEncoder(),
      :transformer_options => Dict()
    })
    fit!(wrapper, fcp.train_instances, fcp.train_labels)
    transformed = transform!(wrapper, fcp.test_instances)

    @fact size(transformed, 2) => 11
  end
end

end # module


text:
# ---------------------------- fit_3c_local ---------------------------------------
# Function to conduct local non-linear least squares curve fitting for each
# pixel in a two-color SAIM image stack
# Images sequences should be saved as a tif stack
# Sequences should be organized as (λ1₁...λ1ₙ λ2₁...λ2ₙ) where n indicates the
# nth and last angle; λ1ᵢ and λ2ᵢ are the ith images acquired with the first and
# and second excitation wavelengths, respectively
# INPUTS:
# name: A string that stores base file name of tif stack
# file_ext: file extension for iumage stack; Only ".tif" has been tested so far
# optic: SAIMOptics structure
# angles: 1D array containing incidence angles in radians for the acquisition
# 	sequence
# init_params: Intial parameters, [A1ₒ, B1ₒ, A2ₒ, B2ₒ, A3ₒ, B3ₒ]
# lower_bounds: Lower bounds for [A1, B1, A2, B2, A3, B3, H]
# upper_bounds: Lower bounds for [A1, B1, A2, B2, A3, B3, H]
# show: When true, a heatmap for the fit heights will be displayed
# OUTPUTS: 1) Heatmap with fitted heights; 2) 16-bit image image with intensities
# proportional to best-fit pixel height -  Pixel intensity = 100*height in nm;
# 3) JLD file with results for each pixel - Fields are "A1," "B1," "A2," "B2," "A3", "B3",
# "H," and "errors"; errors are the standard errors for fit [A1, B1, A2, B2, H]
function fit_3c_local(file_path::String, file_name::String, optic::SAIMOptics,
    angles::AbstractArray, init_params::AbstractArray, lower_bounds::AbstractArray,
	upper_bounds::AbstractArray, show::Bool=false)

    #Add code
end

# ---------------------------- fit_3c_global ----------------------------------
# Function to conduct global non-linear least squares curve fitting for each
# pixel in a three-color SAIM image stack
# Executes a simple grid search, stepping through initial guesses for H; search
# is conducted along H from lower_bounds to upper_bounds
# Images sequences should be saved as a tif stack
# Sequences should be organized as (λ1₁...λ1ₙ λ2₁...λ2ₙ) where n indicates the
# nth and last angle; λ1ᵢ and λ2ᵢ are the ith images acquired with the first and
# and second excitation wavelengths, respectively
# INPUTS:
# name: A string that stores base file name of tif stack
# file_ext: file extension for iumage stack; Only ".tif" has been tested so far
# optic: SAIMOptics structure
# angles: 1D array containing incidence angles in radians for the acquisition
# 	sequence
# init_params: Intial parameters, [A1ₒ, B1ₒ, A2ₒ, B2ₒ, A3ₒ, B3ₒ]
# lower_bounds: Lower bounds for [A1, B1, A2, B2, A3, B3, H]
# upper_bounds: Lower bounds for [A1, B1, A2, B2, A3, B3, H]
# step: Step size for initial guesses for H in grid search
# show: When true, a heatmap for the fit heights will be displayed
# OUTPUTS: 1) Heatmap with fitted heights; 2) 16-bit image image with intensities
# proportional to best-fit pixel height -  Pixel intensity = 100*height in nm;
# 3) JLD file with results for each pixel - Fields are "A1," "B1," "A2," "B2," "A3", "B3",
# "H," and "errors"; errors are the standard errors for fit [A1, B1, A2, B2, H]
function fit_3c_global(file_path::String, file_name::String, optic::SAIMOptics,
    angles::AbstractArray, init_params::AbstractArray, lower_bounds::AbstractArray,
	upper_bounds::AbstractArray, step::Float64=40., show::Bool=false)

    #Add code
end

# ------------------------ Calculate_Constants_3c -----------------------------
# Function that returns constants for calculating Fresnel coefficients (real
# and imag components) and phase constants in two-color SAIM fitting routines
# INPUTS:
# optic: SAIMOptics structure
# angles: 1D array containing incidence angles in radians for the acquisition
# 	sequence
function calculate_constants_3c(optic, angles)

    #Add code
end

#Optical model for three-color experiments: y = A*Pex(H,theta)+B
#p = [A1, B1, A2, B2, A3, B3, H]
#x[i,:] holds the optical model constants for each excitation incidence angle / wavelength combination
function model_3c(x, p)
    #Add code
end

#Generate the analytical Jacobian for model_3c
#p = [A1, B1, A2, B2, A3, B3, H]
#x[i,:] holds the optical model constants for each excitation incidence angle / wavelength combination
function jacobian_model_3c(x,p)
    #Add code
end


text:
module DenseSynchrony

    using LightGraphs, DifferentialEquations

end


text:
# Type file for a Tweet
# Placeholder for now

type Tweet
  #annotations #unused
  #contributors::Nullable{Array{TwContributor}}
  #coordinates::Nullable{TwCoordinate}
  created_at::AbstractString
  current_user_retweet::Dict{AbstractString, Any}
  #entities::Array{TwEntities}
  favorite_count::Nullable{Int}
  favorited::Bool
  filter_level::AbstractString
  #geo # DEPRECEATED
  id::Int64
  id_str::AbstractString
  in_reply_to_screen_name::Nullable{AbstractString}
  in_reply_to_status_id::Nullable{Int64}
  in_reply_to_status_id_str::Nullable{AbstractString}
  in_reply_to_user_id::Nullable{Int64}
  in_reply_to_user_id_str::Nullable{AbstractString}
  lang::AbstractString
  #place::Nullable{TwPlace}
  possibly_sensitive::Nullable{Bool}
  quoted_status_id::Int64
  quoted_status_id_str::AbstractString
  quoted_status::Dict{AbstractString, Any} #Tweet
  #scopes::Dict
  retweet_count::Int
  retweeted::Bool
  retweeted_status::Dict{AbstractString, Any} #Tweet
  source::AbstractString
  text::AbstractString
  truncated::Bool
  user::Dict{AbstractString, Any} #User
  withheld_copyright::Bool
  withheld_in_countries::Array{AbstractString}
  withheld_scope::AbstractString

  Tweet(d::Dict{AbstractString, Any}) = new(
    # annotations #unused
    #get(d,"contributors", Nullable{Array{TwContributor}}()),
    #get(d,"coordinates", Nullable{TwCoordinate}()),
    get(d,"created_at", ""),
    get(d,"current_user_retweet", Dict{AbstractString, Any}()),
    # get(d,"entities", Array{TwEntities}()),
    get(d,"favorite_count", Nullable{Int}()),
    get(d,"favorited", false),
    get(d,"filter_level", ""),
    # get(d,"geo", # DEPRECEATED
    get(d,"id", 0),
    get(d,"id_str", ""),
    get(d,"in_reply_to_screen_name", Nullable{AbstractString}()),
    get(d,"in_reply_to_status_id", Nullable{Int64}()),
    get(d,"in_reply_to_status_id_str", Nullable{AbstractString}()),
    get(d,"in_reply_to_user_id", Nullable{Int64}()),
    get(d,"in_reply_to_user_id_str", Nullable{AbstractString}()),
    get(d,"lang", ""),
    #get(d,"place", Nullable{TwPlace}()),
    get(d,"possibly_sensitive", Nullable{Bool}()),
    get(d,"quoted_status_id", 0),
    get(d,"quoted_status_id_str", ""),
    get(d,"quoted_status",Dict{AbstractString, Any}()),
    #get(d,"scopes", Dict()),
    get(d,"retweet_count", 0),
    get(d,"retweeted", false),
    get(d,"retweeted_status",Dict{AbstractString, Any}()),
    get(d,"source", ""),
    get(d,"text", ""),
    get(d,"truncated", false),
    get(d,"user",Dict{AbstractString, Any}()),
    get(d,"withheld_copyright", false),
    get(d,"withheld_in_countries", Array{AbstractString}()),
    get(d,"withheld_scope","")
  )
end

# Type file for a twitter user

function parse_user(str::AbstractString)
  "parse successful"
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
<!DOCTYPE html>
<declare config=angel_admin_configuration>

            <meta name="viewport" content="width=device-width, initial-scale=1">
            <link rel="stylesheet" href="https://unpkg.com/material-components-web@latest/dist/material-components-web.min.css">
            <link href="https://fonts.googleapis.com/icon?family=Material+Icons" rel="stylesheet">
            <link href="https://fonts.googleapis.com/css?family=Roboto" rel="stylesheet">

            <title>{{ title }} | {{ angel_admin_title }}</title>

            <block name="head"></block>

            <block name="body">
                <block name="app_bar">

                          menu
                          {{ title }}

                        <block name="actions"></block>

                </block>

                            {{ angel_admin_username }}

                            dashboard
                            Dashboard

                                <declare service=config.services[name]>
                                  {{ service.icon }}
                                  {{ service.name }}
                                </declare>

                            exit_to_app
                            Exit Admin Panel

                        <block name="content"></block>

                <block name="scripts"></block>
            </block>

</declare>

max_stars_repo_name: spectra-dart/spectra

max_issues_repo_name: spectra-dart/spectra

max_forks_repo_name: spectra-dart/spectra


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DataGen

makedocs(sitename="DataGen.jl documentation")
deploydocs(
    repo = "github.com/vly/DataGen.jl.git",
    push_preview = true,
    devbranch = "main",
)


text:
using Brillouin, Test
using Brillouin.KPaths: Bravais
using Spglib

@testset "KPath for Spglib cell" begin
    # Trigonal lattice, space group R-3m, 166
    # Example: Bi2Se3 (https://materialsproject.org/materials/mp-541837/)
    a = 1.0
    c = 8.0
    pRs_standard    = [[a*√3/2,  a/2, c/3],
                       [-a*√3/2, a/2, c/3],
                       [0, -a, c/3]]
    pRs_nonstandard = [[a*√3/2, -a/2, c/3],
                       [0, a, c/3],
                       [-a*√3/2, -a/2, c/3]]

    θ = 10 / 180 * pi
    rotation = [[cos(θ) sin(θ) 0]; [-sin(θ) cos(θ) 0]; [0 0 1]]
    pRs_rotated = Ref(rotation) .* pRs_standard

    cell_standard = Spglib.Cell(pRs_standard, [[0, 0, 0]], [0])
    cell_nonstandard = Spglib.Cell(pRs_nonstandard, [[0, 0, 0]], [0])
    cell_rotated = Spglib.Cell(pRs_rotated, [[0, 0, 0]], [0])

    kp_standard = irrfbz_path(cell_standard)
    kp_nonstandard = irrfbz_path(cell_nonstandard)
    kp_rotated = irrfbz_path(cell_rotated)

    @test Bravais.reciprocalbasis(basis(kp_standard)) ≈ pRs_standard
    @test Bravais.reciprocalbasis(basis(kp_nonstandard)) ≈ pRs_nonstandard
    @test Bravais.reciprocalbasis(basis(kp_rotated)) ≈ pRs_rotated

    @test paths(kp_nonstandard) == paths(kp_standard)
    @test paths(kp_rotated) == paths(kp_standard)

    for (lab, kv) in points(kp_standard)
        @test points(kp_nonstandard)[lab] ≈ kv
        @test points(kp_rotated)[lab] ≈ kv
    end

    # The case where the input cell is a supercell of a smaller primitive cell
    pRs_super = [pRs_standard[1], pRs_standard[2], 2 .* pRs_standard[3]]
    cell_super = Spglib.Cell(pRs_super, [[0, 0, 0], [0, 0, 0.5]], [0, 0])
    kp_super = irrfbz_path(cell_super)

    @test Bravais.reciprocalbasis(basis(kp_super)) ≈ pRs_super
    @test paths(kp_super) == paths(kp_standard)

    # Check that the Cartesian coordinates of the high-symmetry k points are identical
    kp_standard_cart = cartesianize(kp_standard)
    kp_super_cart = cartesianize(kp_super)
    for (lab, kv) in points(kp_standard_cart)
        @test points(kp_super_cart)[lab] ≈ kv
    end
end


text:
module GadflyDiff

using ComposeDiff
using Reactive
using Compose
using Gadfly

import Base: writemime

function writemime(io::IO, m::MIME"text/html", ctx::Signal{Gadfly.Plot})
    writemime(io, m, map(c -> draw(
        Patchable(
            Compose.default_graphic_width,
            Compose.default_graphic_height
        ), c), ctx))
end

using Requires

@require IJulia begin
    IJulia.metadata(::Reactive.Signal{Plot}) = Dict()
end

function __init__()
    for f in Requires.__inits__
        f()
    end
end

end # module


text:
push!(LOAD_PATH,"../src/") # load module from local directory

using TheDancer
using Base.Test

init!()

SCFthenUnitary(0.2, 10, 500) # This isn't really a test - it just runs the full GUI!

println("Tests finished succesfully.")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libjlnode"
version = v"14.17.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sunoru/jlnode.git", "d89e011f0ea6e5575c2309e8c53fc93abddffd2c"),
    ArchiveSource("https://github.com/nodejs/node-addon-api/archive/refs/tags/4.0.0.tar.gz", "a61019de219cfbb4943b109fd1c56466c48dedbfcce10567f8e7826992be9c0d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

cd node-addon-api-*
export NAPI_INC=`pwd`

cd ../jlnode
mkdir build
cd build
cmake .. -G"Unix Makefiles" \
    -DCMAKE_BUILD_TYPE="Release" \
    -DCMAKE_LIBRARY_OUTPUT_DIRECTORY="`pwd`/Release" \
    -DCMAKE_JS_INC="${prefix}/include/node" \
    -DNAPI_INC="$NAPI_INC" \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
cmake --build .
cmake --install .
install_license ../LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),

    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),

    # Platform("x86_64", "macos"),
    # Platform("aarch64", "macos"),

    # Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libjlnode", :libjlnode),
    FileProduct("lib/jlnode_addon.node", :jlnode_addon)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("libnode_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS, name, version, sources, script, platforms, products, dependencies;
    julia_compat="1.6",
    preferred_gcc_version=v"6"
)

max_stars_repo_name: remi-garcia/Yggdrasil

max_issues_repo_name: remi-garcia/Yggdrasil

max_forks_repo_name: remi-garcia/Yggdrasil


text:
using Documenter, Perfusion

makedocs(;
    modules = [Perfusion],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Introduction" => "index.md",
        "Pre-processing" => Any[
            "Loading DICOM"=>"preprocessing/loading.md",
            "T1 mapping"=>"preprocessing/T1mapping.md",
            "Signal to Concentration"=>"preprocessing/signal2concentration.md",
        ],
        "AIF-based Models" => Any[
            "AIF"=>"models/aif.md",
            "Tofts"=>"models/tofts.md",
            "Extended Tofts"=>"models/extendedtofts.md",
            "Exchange"=>"models/exchange.md",
            "Uptake"=>"models/uptake.md",
            "Filtration"=>"models/filtration.md",
        ],
        "Functions" => "functions.md",
    ],
    repo = "https://github.com/notZaki/Perfusion.jl/blob/{commit}{path}#L{line}",
    sitename = "Perfusion.jl",
)

# Apply citations
include("apply_citeproc.jl")
apply_citeproc("./build")

deploydocs(; repo = "github.com/notZaki/Perfusion.jl")


text:
# Parameter iterators
postagv(m)=m[1];
deprelv(m)=m[2];
lcountv(m)=m[3];
rcountv(m)=m[4];
distancev(m)=m[5];
bufferm(m)=m[6];
parserv(m)=m[7];

# To initialize parser model from scratch
function initmodel1(o, s; ftype=Float32)
    model = Any[]

    # vector-level initialization
    dpostag, ddeprel, dcount = o[:embed]
    for (k,n,d) in ((:postag,17,dpostag),(:deprel,37,ddeprel),(:lcount,10,dcount),(:rcount,10,dcount),(:distance,10,dcount))
        push!(model, [ initr(d) for i=1:n ])
    end

    # buffer-lstm initialization
    r_b, wr_b = initbuf(o[:bufembed], o[:bufhidden]) # you may change lstm - gru with kwargs
    push!(model, [r_b, wr_b])

    # decision module initialization
    p = o[:arctype](s)
    f = features([p], o[:feats], model)
    mlpdims = (length(f)+r_b.hiddenSize, o[:hidden]..., p.nmove)
    info("mlpdims: $mlpdims")
    decider = Any[]
    for i=2:length(mlpdims)
        push!(decider, initx(mlpdims[i], mlpdims[i-1]))
        push!(decider, initx(mlpdims[i], 1))
    end
    push!(model, decider)
    optims = optimizers(model, o[:optimization])
    return model, optims
end

# To seperate feature vectors and all others
function splitmodel(pmodel)
    # feat-model
    featmodel = Array{Any}(5)
    for k in 1:5 # (:postag,:deprel,:lcount,:rcount,:distance)
        featmodel[k] = Any[]
        pmodel_k = pmodel[k]
        for i in 1:length(pmodel_k)
            push!(featmodel[k], pmodel_k[i])
        end
    end

    # buffer-lstm
    bufmodel = Any[]; buftemp = bufferm(pmodel);
    for i in 1:length(buftemp); push!(bufmodel, buftemp[i]);end;

    # mlp-model
    mlpmodel = Any[]; mlptemp = parserv(pmodel);
    for i in 1:length(mlptemp); push!(mlpmodel, mlptemp[i]); end;
    return (featmodel, bufmodel, mlpmodel)
end

function mlp(w,x; pdrop=(0,0))
    x = dropout(x,pdrop[1])
    for i=1:2:length(w)-2
        x = relu.(w[i]*x .+ w[i+1])
        x = dropout(x,pdrop[2])
    end
    return w[end-1]*x .+ w[end]
end

function scan_buffer(bufmodel, featmodel, parser)
    rbuf, wbuf, pvecs = bufmodel[1], bufmodel[2], postagv(featmodel)
    sentence = parser.sentence
    seqe = length(sentence)
    seqs = parser.wptr
    range = seqe:-1:seqs
    #instate = zeros(rbuf.dataType, rbuf.hiddenSize, 1)
    instate = zeros(Float32, 256, 1)
    ybuf = hout = cout = (gpu() >=0 ? KnetArray(instate) : instate)
    for i in range
        input = vcat(sentence.cavec[i], pvecs[sentence.postag[i]])
        input = (gpu() >=0 ? KnetArray(input) : input)
        ybuf, hout, cout, _ = rnnforw(rbuf, wbuf, input, hout, cout, hy=true, cy=true)
    end
    return ybuf
end

function scan_bufbatch(bufmodel, featmodel, parsers)
    # Todo: we may cache some parsers in future
    yalls = Any[]
    for p in parsers
        ybuf = scan_buffer(bufmodel, featmodel, p)
        push!(yalls, ybuf)
    end
    yalls = cat1d(yalls...)
    ncols = length(parsers)
    nrows = div(length(yalls), ncols)
    return reshape(yalls, nrows, ncols)
end

function oracleloss(allmodel, sentences, arctype, feats; losses=nothing, pdrop=(0.0, 0.0))
    parsers =  map(arctype, sentences)
    mcosts  =  Array{Cost}(parsers[1].nmove)
    parserdone = falses(length(parsers))
    featmodel, bufmodel, mlpmodel = splitmodel(allmodel)
    totalloss = 0.0

    while !all(parserdone)
        fmatrix = features(parsers, feats, featmodel)
        if gpu() >= 0
            fmatrix = KnetArray(fmatrix)
        end
        ybufs = scan_bufbatch(bufmodel, featmodel, parsers)
        input = vcat(fmatrix, ybufs)
        scores = mlp(mlpmodel, input, pdrop=pdrop)
        logprobs = logp(scores, 1)
        for (i, p) in enumerate(parsers)
            if parserdone[i]
                continue
            end
            movecosts(p, p.sentence.head, p.sentence.deprel, mcosts)
            goldmove = indmin(mcosts)

            if mcosts[goldmove] == typemax(Cost)
                parserdone[i] = true
                p.sentence.parse = p
            else
                totalloss -= logprobs[goldmove, i]
                move!(p, goldmove)
                if losses != nothing
                    loss1 = -getval(logprobs)[goldmove,i]
                    losses[1] += loss1
                    losses[2] += 1
                    if losses[2] < 1000
                        losses[3] = losses[1]/losses[2]
                    else
                        losses[3] = 0.999 * losses[3] + 0.001 * loss1
                    end
                end
            end
        end
    end
    return totalloss / length(sentences)
end

oraclegrad = grad(oracleloss)

# Implement test related

# endofparse to use in test
endofparse(p)=(p.sptr == 1 && p.wptr > p.nword)

# TODO: You need to modify score calculations
function oracletest(pmodel, corpus, arctype, feats, batchsize; pdrop=(0.0, 0.0))
    sentbatches = minibatch(corpus, batchsize)
    featmodel, bufmodel, mlpmodel = splitmodel(pmodel)

    for sentences in sentbatches
        parsers = map(arctype, sentences)
        mcosts = Array{Cost}(parsers[1].nmove) #Array(Cost, parsers[1].nmove)
        parserdone = falses(length(parsers))
        while !all(parserdone)
            fmatrix = features(parsers, feats, featmodel)
            if gpu()>=0
                fmatrix = KnetArray(fmatrix)
            end
            ybufs = scan_bufbatch(bufmodel, featmodel, parsers)
            input = vcat(fmatrix, ybufs)
            scores = mlp(mlpmodel, input, pdrop=pdrop)
            logprob = Array(logp(scores, 1))

            for (i,p) in enumerate(parsers)
                if parserdone[i]
                    continue
                end

                isorted = sortperm(logprob[:, i], rev=true) # ith column for ith instance
                for m in isorted # take the best&valid action
                    if moveok(p, m)
                        move!(p, m)
                        break
                    end
                end
                if endofparse(p)
                    parserdone[i] = true
                    p.sentence.parse = p
                end
            end
        end
    end
    return las(corpus)
end

# labeled-attachment score calculator
function las(corpus)
    nword = ncorr = 0
    for s in corpus
        p = s.parse
        nword += length(s)
        ncorr += sum((s.head .== p.head) .& (s.deprel .== p.deprel))
    end
    ncorr / nword
end

function empty_parses!(corpus)
    for s in corpus
        s.parse = nothing
    end
end


text:
# to be a drop in replacement we need to support iterables
struct TakeWhile{Func, Iter}
  f::Func
  iter::Iter
end

function Base.iterate(tw::TakeWhile)
  (nextval, nextstate) = @ifsomething Base.iterate(tw.iter)
  tw.f(nextval) ? (nextval, nextstate) : nothing
end
function Base.iterate(tw::TakeWhile, state)
  (nextval, nextstate) = @ifsomething Base.iterate(tw.iter, state)
  tw.f(nextval) ? (nextval, nextstate) : nothing
end

# to be a drop in replacement we need to support iterables
struct DropWhile{Func, Iter}
  f::Func
  iter::Iter
end
function Base.iterate(dw::DropWhile)
  (nextval, nextstate) = @ifsomething Base.iterate(dw.iter)
  while dw.f(nextval)
    (nextval, nextstate) = @ifsomething Base.iterate(dw.iter, nextstate)
  end
  (nextval, nextstate)
end
Base.iterate(dw::DropWhile, state) = Base.iterate(dw.iter, state)

# IteratorEltype and IteratorSize can be defined for both simultanuously

const DropWhile_or_TakeWhile{F, Iter} = Union{DropWhile{F, Iter}, TakeWhile{F, Iter}}
const TypeDropWhile_or_TypeTakeWhile{F, Iter} = Union{Type{DropWhile{F, Iter}}, Type{TakeWhile{F, Iter}}}

Base.IteratorEltype(::TypeDropWhile_or_TypeTakeWhile{F, Iter}) where {F, Iter} = Base.IteratorEltype(Iter)
Base.eltype(::TypeDropWhile_or_TypeTakeWhile{F, Iter}) where {F, Iter} = Base.eltype(Iter)

# defaulting to Base.SizeUnknown
function Base.IteratorSize(::TypeDropWhile_or_TypeTakeWhile)
    Base.SizeUnknown()
end

function Base.length(tw::DropWhile_or_TakeWhile)
  s::Int = 0
  for x in tw
    s += 1
  end
  s
end

function Base.collect(tw::DropWhile_or_TakeWhile)
  everything = []
  for x in tw
    push!(everything, x)
  end
  everything
end


text:
include("../examples/2D\ Cauchy\ Distribution.jl")
include("../examples/Airy\ equation.jl")
if isdir(Pkg.dir("RandomMatrices")) && isdir(Pkg.dir("PyPlot"))
    include("../examples/GUE\ Sampling.jl")
end
include("../examples/Lanczos.jl")
include("../examples/Lee\ &\ Greengard\ equation.jl")
include("../examples/Newton iteration for Nonlinear BVPs.jl")


text:
"""
    struct DefaultContext <: AbstractContext end

The `DefaultContext` is used by default to compute log the joint probability of the data
and parameters when running the model.
"""
struct DefaultContext <: AbstractContext end

"""
    struct PriorContext{Tvars} <: AbstractContext
        vars::Tvars
    end

The `PriorContext` enables the computation of the log prior of the parameters `vars` when
running the model.
"""
struct PriorContext{Tvars} <: AbstractContext
    vars::Tvars
end
PriorContext() = PriorContext(nothing)

"""
    struct LikelihoodContext{Tvars} <: AbstractContext
        vars::Tvars
    end

The `LikelihoodContext` enables the computation of the log likelihood of the parameters when
running the model. `vars` can be used to evaluate the log likelihood for specific values
of the model's parameters. If `vars` is `nothing`, the parameter values inside the `VarInfo` will be used by default.
"""
struct LikelihoodContext{Tvars} <: AbstractContext
    vars::Tvars
end
LikelihoodContext() = LikelihoodContext(nothing)

"""
    struct MiniBatchContext{Tctx, T} <: AbstractContext
        ctx::Tctx
        loglike_scalar::T
    end

The `MiniBatchContext` enables the computation of
`log(prior) + s * log(likelihood of a batch)` when running the model, where `s` is the
`loglike_scalar` field, typically equal to `the number of data points / batch size`.
This is useful in batch-based stochastic gradient descent algorithms to be optimizing
`log(prior) + log(likelihood of all the data points)` in the expectation.
"""
struct MiniBatchContext{Tctx, T} <: AbstractContext
    ctx::Tctx
    loglike_scalar::T
end
function MiniBatchContext(ctx = DefaultContext(); batch_size, npoints)
    return MiniBatchContext(ctx, npoints/batch_size)
end


text:
module SLOBExec

using SequentialLOB

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    parsed_args = parse_commandline()
    slob_model = SLOB(
        1,
        parsed_args["T"],
        parsed_args["p₀"],
        parsed_args["M"],
        parsed_args["L"],
        parsed_args["D"],
        parsed_args["σ"],
        parsed_args["nu"],
        parsed_args["α"],
        SourceTerm(parsed_args["λ"], parsed_args["μ"]),
    )

    mid_price_paths = slob_model(parsed_args["SEED"])
    print(mid_price_paths[:, 1])
    return 0
end

end


text:
module Layers

using Knet
using Photon: Tensor, Shape

include("core.jl")
include("conv.jl")
include("recurrent.jl")
include("container.jl")

export Dense, Sequential, Flatten, Activation, BatchNorm, Dropout, Concurrent, Softmax, LogSoftmax,
    Conv1D, Conv2D, Conv3D,  Conv2DTranspose, MaxPool2D, AvgPool2D, MaxPool2D, AdaptiveAvgPool, AdaptiveMaxPool,
    RNN, LSTM, GRU, add, output_size,
    relu, elu, selu, sigm # some Knet activation functions

end


text:
"Random Variables are functions from the sample space to some other space"
abstract type RandVar{T} end

"All RandVars have an id which is used to make them independent"
Id = Integer

"A symbolic *canonical* representation of a random variable"
abstract type SymbolicRandVar{T} <: RandVar{T} end

"Can be excuted as a normal julia function"
type JuliaRandVar{T} <: RandVar{T}
  func::Function
  dims::Set{Int}
end

# """An array of random variables (and also a random variable itself)
#   `T` is the range type of elements (e.g for multivariate normal, T = Float64)
#   `N` is the dimensionality of array"""
# type RandArray{T,N} <: DenseArray{T,N}
#   array::Array{RandVar{T},N}
# end
#
# "A matrix of compiled random variables"
# type ExecutableRandArray{T, N} <: DenseArray{T,N}
#   array::Array{JuliaRandVar{T},N}
# end

"The type of the range of a random variable"
rangetype{T}(::RandVar{T}) = T

"Number of dimensions of a random variable"
ndims(X::RandVar) = length(dims(X))

"Set Precision - Generic does nothing"
function set_precision!(::RandVar, precision::Float64) end

## Aliases
Lift{T} = Union{T,SymbolicRandVar{T}}
include("randvar/randarray.jl")
AllRandVars = Union{RandVar, RandArray}

include("randvar/symbolic.jl")
include("randvar/elementary.jl")
include("randvar/executable.jl")

## Specific solver-based randvars
include("randvar/dreal.jl")

# Set the default rrand var
global global_randvar = DRealRandVar
default_randvar() = (global global_randvar; global_randvar)
function set_default_randvar!{T<:RandVar}(R::Type{T})
  println("Setting default solver to $R")
  global global_randvar
  global_randvar = T
end

# Call An Arbitrary Simple Composite type with an ω
# This calling is quite naive, it just calls all the fields of the type
# with ω.  It does not do any recursive lookup
# And the type must be constructable just by using the parameters (default constructor)
# Use call_type instead of call to avoid mass ambiguity
function call_type{T}(X::T, ω::Omega)
  T.abstract && error("Cannot use abstract type as RandVar")
  properties = Any[]
  for fieldname in fieldnames(X)
    field = getfield(X, fieldname)
    if isa(field, AllRandVars)
      push!(properties, call(field, ω))
    else
      push!(properties, field)
    end
  end
  T(properties...)
end


text:
function getthing(mod::Module, name::Vector{Symbol}, default = nothing)
  thing = mod
  for sym in name
    if isdefined(thing, sym)
      thing = getfield(thing, sym)
    else
      return default
    end
  end
  return thing
end

getthing(mod::Module, name::AbstractString, default = nothing) =
  name == "" ?
    default :
    getthing(mod, map(Symbol, split(name, ".", keep=false)), default)


text:
function _empty!(
    tree   :: LeafOrNodeCausalH{S}
    ) where {S}
    if is_leaf(tree)
        tree.inds_build = Int[]
        tree.inds_pred = Int[]
    else
        _empty!(tree.left)
        _empty!(tree.right)
    end
end

function _empty!(
    tree   :: LeafOrNodeCausalNH{S}
    ) where {S}
    if is_leaf(tree)
        tree.inds_build = Int[]
    else
        _empty!(tree.left)
        _empty!(tree.right)
    end
end

function _fill_proj!(
    tree       :: LeafOrNodeCausalH{S},
    inds_build :: Vector{Int},
    inds_pred  :: Vector{Int},
    X          :: AbstractMatrix{S},
    set_u      :: Set{Int}
    ) where {S}

    if is_leaf(tree)
        tree.inds_build = inds_build
        tree.inds_pred = inds_pred
    elseif tree.featid in set_u
        _fill_proj!(tree.left, inds_build, inds_pred, X, set_u)
        _fill_proj!(tree.right, inds_build, inds_pred, X, set_u)
    else
        n_build = length(inds_build)
        n_pred = length(inds_pred)
        build_left = Int[]
        build_right = Int[]
        pred_left = Int[]
        pred_right = Int[]
        for i in 1:n_build
            if X[inds_build[i], tree.featid] < tree.featval
                push!(build_left, inds_build[i])
            else
                push!(build_right, inds_build[i])
            end
        end
        for i in 1:n_pred
            if X[inds_pred[i], tree.featid] < tree.featval
                push!(pred_left, inds_pred[i])
            else
                push!(pred_right, inds_pred[i])
            end
        end
        _fill_proj!(tree.left, build_left, pred_left, X, set_u)
        _fill_proj!(tree.right, build_right, pred_right, X, set_u)
    end
end

function _fill_proj!(
    tree       :: LeafOrNodeCausalNH{S},
    inds_build :: Vector{Int},
    X          :: AbstractMatrix{S},
    set_u      :: Set{Int}
    ) where {S}

    if is_leaf(tree)
        tree.inds_build = inds_build
    elseif tree.featid in set_u
        _fill_proj!(tree.left, inds_build, X, set_u)
        _fill_proj!(tree.right, inds_build, X, set_u)
    else
        n_build = length(inds_build)
        build_left = Int[]
        build_right = Int[]
        for i in 1:n_build
            if X[inds_build[i], tree.featid] < tree.featval
                push!(build_left, inds_build[i])
            else
                push!(build_right, inds_build[i])
            end
        end
        _fill_proj!(tree.left, build_left, X, set_u)
        _fill_proj!(tree.right, build_right, X, set_u)
    end
end

function build_projected_tree(
    tree   :: TreeCausalH{S},
    set_u  :: Set{Int},
    X      :: AbstractMatrix{S}
    ) where {S}

    inds_build = tree.inds_build
    inds_pred = tree.inds_pred
    root = deecopy(tree.tree) # TODO verif si pas deepcopy plutot si juste copy il faut la definir a la main  Base.copy(s::Nodeetc)

    _empty!(root)
    _fill_proj!(root, inds_build, inds_pred, X, set_u)

    return TreeCausalH{S}(root, inds_build, inds_pred, tree.oob)

end

function build_projected_tree(
    tree   :: TreeCausalNH{S},
    set_u  :: Set{Int},
    X      :: AbstractMatrix{S}
    ) where {S}

    inds = tree.inds
    root = deepcopy(tree.tree) # TODO verif si pas deepcopy plutot

    _empty!(root)
    _fill_proj!(root, inds, X, set_u)

    return TreeCausalNH{S}(root, inds, tree.oob)

end

function apply_forest_oob(
    forest :: EnsembleCausal{S}
    ) where {S}

    centering = forest.centering
    n_samples = length(forest.Y)
    n_trees = length(forest)
    predictions = Array{Float64}(undef, n_samples)
    if centering
        Y_vec =  forest.Y_center
        T_vec =  forest.T_center
    else
        Y_vec = forest.Y
        T_vec = forest.T
    end

    for j in 1:n_samples

        alpha = zeros(n_samples)
        nb_oob_samples = 0

        for b in 1:n_trees
            if j in forest.trees[b].oob
                nb_oob_samples += 1
                N = neighbours(forest, b, forest.X[j, :]) #TODO attention x avant mais mort  --> forest.X[j,:] ??
                l = length(N)
                for e in N
                    alpha[e] += 1/l
                end
            end
        end
        alpha = alpha/nb_oob_samples
        T_alpha = 0
        Y_alpha = 0
        for i in 1:n_samples
            T_alpha += alpha[i]*T_vec[i]
            Y_alpha += alpha[i]*Y_vec[i]
        end
        num = 0
        denom = 0
        for i in 1:n_samples
            diffT = T_vec[i]-T_alpha
            num += alpha[i]*diffT*(Y_vec[i]-Y_alpha)
            denom += alpha[i]*diffT^2
        end
        predictions[j] = num/denom
    end
    return predictions
end

function prf_causal_effect( # causal effect on oob
    forest   :: EnsembleCausal{S},
    set_u    :: Set{Int}
    ) where {S}

    n_trees = length(forest.trees)

    if forest.trees isa Vector{TreeCausalH{S}}
        trees = Vector{TreeCausalH{S}}(undef, n_trees)
    else
        trees = Vector{TreeCausalNH{S}}(undef, n_trees)
    end

    for i in 1:n_trees
        trees[i] = build_projected_tree(forest.trees[i], set_u, forest.X)
    end

    prf = EnsembleCausal{S}(trees, forest.centering, forest.bootstrap, forest.honest,
        forest.X, forest.Y, forest.T, forest.model_Y, forest.model_T, forest.Y_center,
        forest.T_center)

    return apply_forest_oob(prf)
end