text:
immutable EigenprecodingState
    V::Array{Matrix{Complex128},1}
end

function Eigenprecoding(channel::SinglecarrierChannel, network)
    assignment = get_assignment(network)

    I = get_num_BSs(network)
    Ps = get_transmit_powers(network)
    sigma2s = get_receiver_noise_powers(network)
    alphas = get_user_priorities(network)
    aux_params = get_aux_precoding_params(network)

    state = EigenprecodingState(Array(Matrix{Complex128}, channel.K))

    for i = 1:I
        served = served_MS_ids(i, assignment)
        Kc = length(served)

        for k in served
            _, sing_vals, r_sing_vecs = svd(channel.H[k,i], thin=true)
            sort!(sing_vals, rev=true) # just in case
            max_d = length(sing_vals)

            # Due to sorting, the strongest channel is first
            noise_over_channel_power = sigma2s[k]./sing_vals

            # Activate all channels
            active_channels = collect(1:max_d)

            Palloc = zeros(Float64, max_d)
            while true
                # Current water level and power allocation
                mu = (1/length(active_channels))*(Ps[i]/Kc + sum(noise_over_channel_power[active_channels]))
                Palloc = mu - noise_over_channel_power[active_channels]

                # If water level high enough, done. Otherwise, turn off weakest channel
                Palloc[end] > 0 ? break : pop!(active_channels)
            end

            # Final precoder for this user
            state.V[k] = r_sing_vecs*diagm(vcat(sqrt(Palloc), zeros(Float64, max_d - length(Palloc))))
        end
    end

    return calculate_logdet_rates(state, channel, sigma2s, alphas, assignment, aux_params)
end

function calculate_logdet_rates(state::EigenprecodingState,
    channel::SinglecarrierChannel, sigma2s, alphas, assignment, aux_params)

    max_d = min(maximum(channel.Ns), maximum(channel.Ms)) # might not be tight..

    if aux_params["output_protocol"] == :all_iterations
        intercell_tdma_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        intracell_tdma_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        uncoord_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])

        intercell_tdma_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        intracell_tdma_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        uncoord_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])

        intercell_tdma_weighted_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        intracell_tdma_weighted_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        uncoord_weighted_logdet_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])

        intercell_tdma_weighted_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        intracell_tdma_weighted_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
        uncoord_weighted_MMSE_rates = Array(Float64, channel.K, max_d, aux_params["max_iters"])
    elseif aux_params["output_protocol"] == :final_iteration
        intercell_tdma_logdet_rates = Array(Float64, channel.K, max_d)
        intracell_tdma_logdet_rates = Array(Float64, channel.K, max_d)
        uncoord_logdet_rates = Array(Float64, channel.K, max_d)

        intercell_tdma_MMSE_rates = Array(Float64, channel.K, max_d)
        intracell_tdma_MMSE_rates = Array(Float64, channel.K, max_d)
        uncoord_MMSE_rates = Array(Float64, channel.K, max_d)

        intercell_tdma_weighted_logdet_rates = Array(Float64, channel.K, max_d)
        intracell_tdma_weighted_logdet_rates = Array(Float64, channel.K, max_d)
        uncoord_weighted_logdet_rates = Array(Float64, channel.K, max_d)

        intercell_tdma_weighted_MMSE_rates = Array(Float64, channel.K, max_d)
        intracell_tdma_weighted_MMSE_rates = Array(Float64, channel.K, max_d)
        uncoord_weighted_MMSE_rates = Array(Float64, channel.K, max_d)
    end

    for i = 1:channel.I
        served = served_MS_ids(i, assignment)
        Kc = length(served)

        for k in served
            Phi_intracell = complex(sigma2s[k]*eye(channel.Ns[k]))
            Phi_uncoord = complex(sigma2s[k]*eye(channel.Ns[k]))
            for j in delete!(IntSet(1:channel.I), i)
                for l in served_MS_ids(j, assignment)
                    F = channel.H[k,j]*state.V[l]
                    FFh = F*F'
                    Phi_intracell += FFh
                    Phi_uncoord += FFh
                end
            end
            for l in served_MS_ids_except_me(k, i, assignment)
                F = channel.H[k,i]*state.V[l]
                Phi_uncoord += F*F'
            end

            d = size(state.V[k], 2)
            W_intercell = UniformScaling(1.) + (channel.K/Kc)*state.V[k]'*channel.H[k,i]'*(1/sigma2s[k])*channel.H[k,i]*state.V[k]
            W_intracell = UniformScaling(1.) + state.V[k]'*channel.H[k,i]'*(Phi_intracell\channel.H[k,i])*state.V[k]
            W_uncoord = UniformScaling(1.) + state.V[k]'*channel.H[k,i]'*(Phi_uncoord\channel.H[k,i])*state.V[k]

            r_intercell_logdet = (1/channel.K)*log2(max(1, abs(eigvals(W_intercell))))
            r_intracell_logdet = (1/Kc)*log2(max(1, abs(eigvals(W_intracell))))
            r_uncoord_logdet = log2(max(1, abs(eigvals(W_uncoord))))

            r_intercell_MMSE = (1/channel.K)*log2(max(1, abs(1./diag(inv(W_intercell)))))
            r_intracell_MMSE = (1/Kc)*log2(max(1, abs(1./diag(inv(W_intracell)))))
            r_uncoord_MMSE = log2(max(1, abs(1./diag(inv(W_uncoord)))))

            r_intercell_weighted_logdet = alphas[k]*r_intercell_logdet
            r_intracell_weighted_logdet = alphas[k]*(1/Kc)*r_intracell_logdet
            r_uncoord_weighted_logdet = alphas[k]*r_uncoord_logdet

            r_intercell_weighted_MMSE = alphas[k]*r_intercell_MMSE
            r_intracell_weighted_MMSE = alphas[k]*r_intracell_MMSE
            r_uncoord_weighted_MMSE = alphas[k]*r_uncoord_MMSE

            if aux_params["output_protocol"] == :all_iterations
                for iter = 1:aux_params["max_iters"]
                    intercell_tdma_logdet_rates[k,:,iter] = cat(1, r_intercell_logdet, zeros(Float64, max_d - d))
                    intracell_tdma_logdet_rates[k,:,iter] = cat(1, r_intracell_logdet, zeros(Float64, max_d - d))
                    uncoord_logdet_rates[k,:,iter] = cat(1, r_uncoord_logdet, zeros(Float64, max_d - d))

                    intercell_tdma_MMSE_rates[k,:,iter] = cat(1, r_intercell_MMSE, zeros(Float64, max_d - d))
                    intracell_tdma_MMSE_rates[k,:,iter] = cat(1, r_intracell_MMSE, zeros(Float64, max_d - d))
                    uncoord_MMSE_rates[k,:,iter] = cat(1, r_uncoord_MMSE, zeros(Float64, max_d - d))

                    intercell_tdma_weighted_logdet_rates[k,:,iter] = cat(1, r_intercell_weighted_logdet, zeros(Float64, max_d - d))
                    intracell_tdma_weighted_logdet_rates[k,:,iter] = cat(1, r_intracell_weighted_logdet, zeros(Float64, max_d - d))
                    uncoord_weighted_logdet_rates[k,:,iter] = cat(1, r_uncoord_weighted_logdet, zeros(Float64, max_d - d))

                    intercell_tdma_weighted_MMSE_rates[k,:,iter] = cat(1, r_intercell_weighted_MMSE, zeros(Float64, max_d - d))
                    intracell_tdma_weighted_MMSE_rates[k,:,iter] = cat(1, r_intracell_weighted_MMSE, zeros(Float64, max_d - d))
                    uncoord_weighted_MMSE_rates[k,:,iter] = cat(1, r_uncoord_weighted_MMSE, zeros(Float64, max_d - d))
                end
            elseif aux_params["output_protocol"] == :final_iteration
                intercell_tdma_logdet_rates[k,:] = cat(1, r_intercell_logdet, zeros(Float64, max_d - d))
                intracell_tdma_logdet_rates[k,:] = cat(1, r_intracell_logdet, zeros(Float64, max_d - d))
                uncoord_logdet_rates[k,:] = cat(1, r_uncoord_logdet, zeros(Float64, max_d - d))

                intercell_tdma_MMSE_rates[k,:] = cat(1, r_intercell_MMSE, zeros(Float64, max_d - d))
                intracell_tdma_MMSE_rates[k,:] = cat(1, r_intracell_MMSE, zeros(Float64, max_d - d))
                uncoord_MMSE_rates[k,:] = cat(1, r_uncoord_MMSE, zeros(Float64, max_d - d))

                intercell_tdma_weighted_logdet_rates[k,:] = cat(1, r_intercell_weighted_logdet, zeros(Float64, max_d - d))
                intracell_tdma_weighted_logdet_rates[k,:] = cat(1, r_intracell_weighted_logdet, zeros(Float64, max_d - d))
                uncoord_weighted_logdet_rates[k,:] = cat(1, r_uncoord_weighted_logdet, zeros(Float64, max_d - d))

                intercell_tdma_weighted_MMSE_rates[k,:] = cat(1, r_intercell_weighted_MMSE, zeros(Float64, max_d - d))
                intracell_tdma_weighted_MMSE_rates[k,:] = cat(1, r_intracell_weighted_MMSE, zeros(Float64, max_d - d))
                uncoord_weighted_MMSE_rates[k,:] = cat(1, r_uncoord_weighted_MMSE, zeros(Float64, max_d - d))
            end
        end
    end

    results = PrecodingResults()
    results["intercell_tdma_logdet_rates"] = intercell_tdma_logdet_rates
    results["intracell_tdma_logdet_rates"] = intracell_tdma_logdet_rates
    results["uncoord_logdet_rates"] = uncoord_logdet_rates
    results["intercell_tdma_MMSE_rates"] = intercell_tdma_MMSE_rates
    results["intracell_tdma_MMSE_rates"] = intracell_tdma_MMSE_rates
    results["uncoord_MMSE_rates"] = uncoord_MMSE_rates
    results["intercell_tdma_weighted_logdet_rates"] = intercell_tdma_weighted_logdet_rates
    results["intracell_tdma_weighted_logdet_rates"] = intracell_tdma_weighted_logdet_rates
    results["uncoord_weighted_logdet_rates"] = uncoord_weighted_logdet_rates
    results["intercell_tdma_weighted_MMSE_rates"] = intercell_tdma_weighted_MMSE_rates
    results["intracell_tdma_weighted_MMSE_rates"] = intracell_tdma_weighted_MMSE_rates
    results["uncoord_weighted_MMSE_rates"] = uncoord_weighted_MMSE_rates
    return results
end


text:
function get_result(rdir::AbstractString, var::Union{String, Symbol}; filter_by::Dict=Dict())

    df = joinpath(rdir,String(var)*".feather") |> Feather.read

    for n in names(df)
        if typeof(df[!,n]) <: Arrow.BitPrimitive
            df[!,n] = df[!,n] |> collect
        end
    end

    for (k,v) in filter_by
        if typeof(v) <: AbstractArray
            df = filter(x-> x[k] in v, df)
        else
            df = filter(x-> x[k] == v, df)
        end
    end

    return df
end

function c_gradient(g::Array{RGB{Float64}, 1}, n::Int)
    if n > 1
        return RGB[get(ColorScheme(g), z) for z=range(1, stop=0, length=n)]
    elseif n==1
        return RGB[get(ColorScheme(g), 1.0)]
    end
end

function c_gradient(g::ColorScheme, n::Int)
    if n > 1
        return RGB[get(g, z) for z=range(1, stop=0, length=n)]
    elseif n==1
        return RGB[get(g, 1.0)]
    end
end

function create_filter_dict(sector::Symbol)
    filter_dict = Dict(:ev => 0, :heat => 0, :h2 => 0)
    delete!(filter_dict, sector)
    return filter_dict
end

function plot_generation_investments(rdir::AbstractString,
    color_dict::Dict;
    techs=["Solar PV", "Wind Onshore", "Wind Offshore", "Biomasse"],
    sector::Symbol,
    diff=false,
    kwargs...)

    capacity = get_result(rdir, :N, filter_by=create_filter_dict(sector))

    renewable = filter(row-> row[:fuel] in techs, capacity)
    sort!(renewable, [:min_res, sector])

    renewable = by(renewable, [:min_res, sector, :fuel], :Value => sum)
    # Note: The above line creates a column-name "Value_sum" in the results.

    plt = plot(legend=:topleft)
    fuels = unique(renewable[!,:fuel])
    len_sector = unique(capacity[!,sector]) |> length

    for f in fuels
        df = renewable[renewable.fuel .== f, :]
        c = c_gradient(color_dict[f], len_sector) |> permutedims
        x = df[!,:min_res]
        y =  df[!,:Value_sum]
        g = df[!,sector]
        l = vcat(f, fill("", len_sector-1)) |> permutedims
        t = "Generation Capacity"
        if sum(df[!,:Value_sum]) > 0
            plot!(x, y; group=g, c=c, label=l, width=2, title=t, grid=false, kwargs...)
        end
    end

    return plt
end

function plot_storage_investments(rdir::AbstractString,
    color_dict::Dict;
    sector::Symbol,
    kwargs...)

    fd = create_filter_dict(sector)
    storage_power = get_result(rdir, :N_STO_P, filter_by=fd)
    storage_energy = get_result(rdir, :N_STO_E, filter_by=fd)
    p2g = get_result(rdir, :N_P2G, filter_by=fd)
    g2p = get_result(rdir, :N_G2P, filter_by=fd)
    gs = get_result(rdir, :N_GS, filter_by=fd)

    for df in [storage_power, storage_energy, p2g, g2p, gs]
        if !isempty(df)
            sort!(df, [:min_res, sector])
        end
    end

    plt = plot(legend=:topleft)
    len_sector = unique(storage_power[!,sector]) |> length

    storage_power = by(storage_power, [:min_res, sector], :Value => sum)
    x = storage_power[!,:min_res]
    y = storage_power[!,:Value_sum]
    g = storage_power[!,sector]
    c = c_gradient(color_dict["Storages"], len_sector)
    plt = plot(x, y; group=g, color=permutedims(c), width=2,
        title="Storage Power", legend=false, grid=false, kwargs...)

    if !isempty(p2g)
        p2g = by(p2g, [:min_res, sector], :Value => sum)
        x = p2g[!,:min_res]
        y = p2g[!,:Value_sum]
        g = p2g[!,sector]
        c = c_gradient(color_dict["Hydrogen"], len_sector)
        plot!(plt, x, y; group=g, color=permutedims(c), width=2,
            label="Power to gas", legend=false, grid=false, linestyle=:dash)
    end

    if !isempty(g2p)
        g2p = by(g2p, [:min_res, sector], :Value => sum)
        x = g2p[!,:min_res]
        y = g2p[!,:Value_sum]
        g = g2p[!,sector]
        c = c_gradient(color_dict["Hydrogen"], len_sector)
        plot!(plt, x,y; group=g, color=permutedims(c), width=2,
            label="Gas to Power", legend=false, grid=false, linestyle=:dot)
    end

    storage_energy = by(storage_energy, [:min_res, sector], :Value => sum)
    x = storage_energy[!,:min_res]
    y = storage_energy[!,:Value_sum]
    g = storage_energy[!,sector]
    c = c_gradient(color_dict["Storages"], len_sector)
    plt2 = plot(x, y; group=g, color=permutedims(c), width=2, legend=:topleft,
        title="Storage Energy", grid=false, kwargs...)

    # if !isempty(gs)
    #     gs = by(gs, [:min_res, sector], :Value => sum)
    #     x = gs[:min_res]
    #     y = gs[:Value_sum]
    #     g = gs[sector]
    #     c = c_gradient(color_dict["Hydrogen"], len_sector)
    #     plot!(plt2, x,y; group=g, color=permutedims(c), width=2,
    #         label="h2 Storage", legend=false, grid=false, linestyle=:dash)
    # end

    return plt, plt2
end

function plot_curtailment(rdir::AbstractString,
    color_dict::Dict;
    sector::Symbol,
    kwargs...)

    fd = create_filter_dict(sector)
    cu = get_result(rdir, "summary_generation", filter_by=fd)
    sort!(cu, [:min_res, sector])

    plt = plot(legend=false)
    len_sector = unique(cu[!,sector]) |> length

    x = cu[!,:min_res]
    y = cu[!,:relative_curtailment]
    g = cu[!,sector]
    c = c_gradient(color_dict["Curtailment"], len_sector)
    plt = plot(x, y; group=g, color=permutedims(c), width=2, legend=:topleft,
        title="Curtailment [%]", grid=false, kwargs...)

    return plt
end

function plot_all(rdir::AbstractString,
    c::Dict;
    sector::Symbol,
    kwargs...)

    cap = plot_generation_investments(rdir, c, sector=sector)
    pow, ener = plot_storage_investments(rdir, c, sector=sector)
    cur = plot_curtailment(rdir, c, sector=sector)

    plt = plot(cap, cur, pow, ener; grid=(2,2), kwargs...)
    return plt
end

ext: jl


text:
# Tested in Julia 1.5.1
using DelimitedFiles, Unicode
data = readdlm("./data/adventInput04.txt", ',', skipblanks= false)
expectedFields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid", "cid"]

function passportsDict(data)
    allpassports = []
    for row in data
        push!(allpassports, split(row, ' ')...)
    end
    passports = Dict()
    passKeys = Dict{String,String}()
    c = 1
    lenpass = length(allpassports)
    for i in 1:lenpass
        keyval = allpassports[i]
        if keyval == "" || i == lenpass
            if i == lenpass
               key, val = split(keyval, ':')
                passKeys[key] = val
            end
            passports["Passport_$(c)"] = passKeys
            passKeys = Dict{String,String}()
            c += 1
        else
            key, val = split(keyval, ':')
            passKeys[key] = val
        end
    end
    passports
end
function validpassport(passport; m = "cid" )
    for k in expectedFields
        keyin = haskey(passport, k)
       if keyin == false && k == m
           return true
        elseif keyin == false
            return false
        elseif keyin == true && k == m
            return true
        end
    end
end
function validCount(passports)
    c = 0
    for (k, passport) in passports
        c += validpassport(passport)
    end
    c
end

passports = passportsDict(data)
validCount(passports)
# More rules !
eyeColors = ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
paramsDigits = Dict()
paramsDigits["byr"] = [1920,2002]
paramsDigits["iyr"] = [2010,2020]
paramsDigits["eyr"] = [2020,2030]

function validDigits(passport; key = "byr", min = 1920, max = 2002)
    is4digits = tryparse(Int, passport[key])
    if is4digits == nothing
        return false
    elseif min <= is4digits <= max
        return true
    else
        return false
    end
end
function validYears(passport)
    cdigtis = 0
    for (k, v)  in paramsDigits
        c = validDigits(passport; key = k, min = v[1], max = v[2])
        if c == 1
            cdigtis += 1
        else
            break
        end
    end
    cdigtis == 3 ? true : false
end
function validUnits(passport)
    unitsRange = Dict()
    unitsRange["in"] = [59,76]
    unitsRange["cm"] = [150,193]

    digit = []
    for s in passport["hgt"]
        if isdigit(s)
            push!(digit, "$(s)")
        end
    end
    isin = occursin(prod(digit), passport["hgt"])
    digit = parse(Int, prod(digit))

    units = [match.match for match in eachmatch(r"(in|cm)", passport["hgt"])]
    if isin == true && length(units) > 0
        cm_in = units[1]
        min, max = unitsRange[cm_in][1], unitsRange[cm_in][2]
        if cm_in == "in" &&  min <= digit <= max
            return true
        elseif cm_in == "cm" &&  min <= digit <= max
            return true
        else
            return false
        end
    elseif isin == false || length(units) == 0
        return false
    end
end

function validColor(passport)
    color = passport["hcl"]
    if color[1] == '#'
        return checkChars(color)
    else
        return false
    end
end
function checkChars(color; lenChars = 6, minchar = 'a', maxchar = 'f')
    vc = 0
    for c in color[2:end]
        if isdigit(c) == false && isuppercase(c) == false && minchar <= c <= maxchar
            vc += 1
        elseif isdigit(c) == true && 0 <= parse(Int,c) <= 9
            vc += 1
        else
            break
        end
    end
    vc == lenChars ? true : false
end
validEyeColor(passport; eyeColors = eyeColors) = passport["ecl"] in eyeColors
function validPid(passport; numdigit = 9)
   lendigit = 0
    for d in passport["pid"]
        lendigit += isdigit(d)
    end
    lendigit == numdigit ? true : false
end
function validCountRules(passports)
    c = 0
    for (k, p) in passports
        if validpassport(p)
            isvalid = validYears(p) + validUnits(p) + validColor(p) + validEyeColor(p) + validPid(p) == 5
            c += isvalid
        end
    end
    c
end
validCountRules(passports)


text:
function combine_spouses_data(spouse_value_matrix)
        n_ins = length(spouse_value_matrix[:,1])

        combined_data = Array(Any,n_ins)
        for i = 1 : n_ins
                combined_data[i] = tuple(spouse_value_matrix[i,:]...)
        end

        return combined_data
end

function parent_class_combine(data_matrix,parent_set)
        n_parent = length(parent_set)
        N = length(data_matrix[:,1])

        parent_data_set = Array(Int64,N,n_parent)

        for i = 1 : n_parent
                parent_data_set[:,i] = data_matrix[:,parent_set[i]]
        end

        parent_data_set = combine_spouses_data(parent_data_set)

        class_combine = length(unique(parent_data_set))
        return class_combine
end

function log_prob_single_edge_last_term(class)
        n = length(class)
        uniq_cl = unique(class)
        val_cl = length(uniq_cl)

        distr_table = Array(Array,n)

        # Survey distribution for each data point

        for index_data = 1 : n
                class_cl = findfirst(uniq_cl,class[index_data])
                z = zeros(Int64,val_cl)
                z_modify = zeros(Int64,val_cl)
                z_modify[class_cl] = 1
                distr_table[index_data] = z_modify
        end

        # Survey class distribution between point i and point j

        distr_intval_table = Array(Array,n,n)

        for ind_init = 1 : n
                for ind_end = ind_init : n
                        if ind_init == ind_end
                                current_distr = distr_table[ind_init]
                        else
                                current_distr = distr_intval_table[ind_init,ind_end-1] +
                                                distr_table[ind_end]
                        end
                        distr_intval_table[ind_init,ind_end] = current_distr
                end
        end
        # println(distr_intval_table[1,n])
        # Get 1/P(D|M) for single edge case and count only last term in objective func

        inv_p = Array(Float64,n,n)

        for ind_init = 1 : n
                for ind_end = 1 : n
                        if ind_end < ind_init
                                inv_p[ind_init,ind_end] = Inf
                        else
                                distr_in_intval = distr_intval_table[ind_init,ind_end]
                                current = lfact(ind_end-ind_init+1)
                                for ind_cl = 1 : val_cl
                                        current -= lfact(distr_in_intval[ind_cl])
                                end

                                inv_p[ind_init,ind_end] = current
                        end
                end
        end

        return inv_p
end

function log_prob_spouse_child_data(child,spouse)
        n = length(child)
        uniq_sp = unique(spouse)
        uniq_ch = unique(child)
        val_sp = length(uniq_sp)
        val_ch = length(uniq_ch)

        distr_table = Array(Array,n,val_sp)

        # Survey distribution for each data point

        for index_data = 1 : n
                class_ch = findfirst(uniq_ch,child[index_data])
                class_sp = findfirst(uniq_sp,spouse[index_data])
                z = zeros(Int64,val_ch)
                z_modify = zeros(Int64,val_ch)
                z_modify[class_ch] = 1

                for class_sp_index = 1 : val_sp
                        if class_sp_index == class_sp
                                distr_table[index_data,class_sp_index] = z_modify
                        else
                                distr_table[index_data,class_sp_index] = z
                        end
                end
        end

        # Survey class distribution in between point i and point j

        distr_intval_table = Array(Array,n,n,val_sp)

        for ind_val_sp = 1 : val_sp
                for ind_init = 1 : n
                for ind_end = ind_init : n
                        if ind_init == ind_end
                                current_distr = distr_table[ind_end,ind_val_sp]
                        else
                                current_distr = distr_intval_table[ind_init,ind_end-1,ind_val_sp] +
                                                distr_table[ind_end,ind_val_sp]
                        end
                        #println(current_distr)
                        distr_intval_table[ind_init,ind_end,ind_val_sp] = current_distr
                end
                end
        end

        # Get the 1/P(D|M) for this child-spouse set

        inv_p = Array(Float64,n,n)

        for ind_init = 1 : n
        for ind_end = 1 : n
                if ind_end < ind_init
                        inv_p[ind_init,ind_end] = Inf
                else
                        current_val = 0.0
                        for ind_val_sp = 1 : val_sp
                                distr_in_intval = distr_intval_table[ind_init,ind_end,ind_val_sp]
                                total_num = sum(distr_in_intval)

                                # The forth term in objective function
                                current = lfact(total_num)
                                for ind_val_ch = 1 : val_ch
                                        current -= lfact(distr_in_intval[ind_val_ch])
                                end

                                # The third term in objective function
                                current += lfact(total_num + val_ch -1) - lfact(val_ch-1) - lfact(total_num)

                                # The 3rd and 4th terms for this value of spouse
                                current_val += current
                        end
                        inv_p[ind_init,ind_end] = current_val
                end
        end
        end

        return inv_p
end

function prior_of_intval(continuous,lambda)
        N = length(continuous)
        d_1_N = continuous[N] - continuous[1]
        prior = Array(Float64,N)
        for i = 1 : N-1
                d_i = continuous[i+1] - continuous[i]
                prior[i] = 1 - exp(-lambda * d_i / d_1_N)
        end
        prior[N] = 1

        return prior
end

function largest_class_value(data_matrix)
        n = length(data_matrix[1,:])
        largest = 0
        for i = 1 : n
                class_value = length(unique(data_matrix[:,i]))
                if class_value > largest
                        largest = class_value
                end
        end

        return largest
end

function BN_discretizer_p_data_model(data_matrix,parent_set,child_spouse_set,approx = true)

        N = length(data_matrix[:,1])
        n = length(data_matrix[1,:])
        n_p = length(parent_set)
        n_c = length(child_spouse_set)
        n_r = n_p + n_c

        parent_class_number = 1
        if length(parent_set) > 0
                for i = 1 : n_p
                        multi = length(unique(data_matrix[:,parent_set[i]]))
                        parent_class_number = parent_class_number * multi
                end

                parent_class_number_2 = parent_class_combine(data_matrix,parent_set)
        end

        # -log(P(D|M)) part:

        log_P_data_model = zeros(Float64,N,N)
        nearest_var_set = [parent_set;child_spouse_set]
        if approx == true
        for ind_parent = 1 : n_p
                table = 0
                single_variable_data = data_matrix[:,parent_set[ind_parent]]
                table = log_prob_single_edge_last_term(single_variable_data)
                log_P_data_model += table
        end

        else
                if n_p > 0
                table = 0
                parent_data = Array(Int64,N,n_p)
                for p = 1 : n_p
                        parent_data[:,p] = data_matrix[:,parent_set[p]]
                end
                combined_parent = combine_spouses_data(parent_data)
                table = log_prob_single_edge_last_term(combined_parent)
                log_P_data_model += table
                end
        end

        for ind_child = 1 : n_c
                table = 0
                if length(child_spouse_set[ind_child]) == 1
                        child_data = data_matrix[:,child_spouse_set[ind_child]]
                        spouse_data = zeros(Int64,N)
                        table = log_prob_spouse_child_data(child_data,spouse_data)

                elseif length(child_spouse_set[ind_child]) == 2
                        (child,spouse) = child_spouse_set[ind_child]
                        child_data = data_matrix[:,child]
                        spouse_data = data_matrix[:,spouse]
                        table = log_prob_spouse_child_data(child_data,spouse_data)
                else
                        # A child has more than 2 parents
                        child = child_spouse_set[ind_child][1]
                        child_data = data_matrix[:,child]
                        spouse_set = child_spouse_set[ind_child][2:end]
                        spouse_matrix = Array(Int64,N,length(spouse_set))
                        for spouse_index = 1 : length(spouse_set)
                                spouse_matrix[:,spouse_index] = data_matrix[:,spouse_set[spouse_index]]
                        end
                        spouse_data = combine_spouses_data(spouse_matrix)
                        table = log_prob_spouse_child_data(child_data,spouse_data)
                end

                log_P_data_model += table
        end

        # Add -log(p_class_distribution_for_parent)

        for i = 1 : N
                for j = i : N
                        log_P_data_model[i,j] += lfact(j - i + parent_class_number) -
                                                 lfact(j-i+1) - lfact(parent_class_number-1)
                end
        end

        return log_P_data_model

end

function BN_discretizer_free_number_rep(continuous,data_matrix,parent_set,child_spouse_set,approx = true)
        p_data_model = BN_discretizer_p_data_model(data_matrix,parent_set,child_spouse_set,approx)
        #lambda = div(largest_class_value(data_matrix),2)
        lambda = largest_class_value(data_matrix)
        split_on_intval = prior_of_intval(continuous,lambda)
        N = length(continuous)
        not_split_on_intval = ones(Float64,N) - split_on_intval

        conti_norep = unique(continuous)
        conti_head = Array(Int64,length(conti_norep))
        conti_tail = Array(Int64,length(conti_norep))
        N_norep = length(conti_norep)

        conti_head[1] = 1
        index_head = 1
        for i = 2 : N
                if continuous[i] != continuous[i-1]
                        index_head += 1
                        conti_head[index_head] = i
                end
        end

        conti_tail[end] = N
        index_tail = length(conti_norep)
        for i = N-1 : -1 : 1
                if continuous[i] != continuous[i+1]
                        index_tail -= 1
                        conti_tail[index_tail] = i
                end
        end

        #println(conti_tail==conti_head)
        smallest_value = Array(Float64,N_norep)
        optimal_disc = Array(Array,N_norep)
        length_conti_data = continuous[N] - continuous[1]

        for a = 1 : N_norep
                if a == 1
                        smallest_value[1] = -log(split_on_intval[conti_tail[1]]) +
                                                        p_data_model[conti_head[1],conti_tail[1]]
                        optimal_disc[1] = [conti_tail[1]]
                else
                        current_value = Inf
                        current_disc = [0]
                        for b = 1 : a
                                if b == a
                                        value =
                                        ( ( (continuous[conti_tail[a]] - continuous[1])/
                                                length_conti_data)*lambda) -
                                        log(split_on_intval[conti_tail[a]]) +
                                        p_data_model[1,conti_tail[a]]
                                else

                                        value = smallest_value[b] +
                                        ( ( (continuous[conti_tail[a]] - continuous[conti_head[b+1]])/
                                        length_conti_data)*lambda) -
                                        log(split_on_intval[conti_tail[a]]) +
                                        p_data_model[conti_head[b+1],conti_tail[a]]
                                end

                                if value < current_value
                                        current_value = value
                                        if b == a
                                                current_disc = [conti_tail[a]]
                                        else
                                                current_disc = [optimal_disc[b];conti_tail[a]]
                                        end
                                end
                        end
                        smallest_value[a] = current_value
                        optimal_disc[a] = current_disc
                end
        end

        #println(optimal_disc[end])
        optimal_disc_value = Array(Float64,length(optimal_disc[end])+1)

        for i = 0 : length(optimal_disc[end])
                if i == 0
                        optimal_disc_value[i+1] = continuous[1]
                elseif i == length(optimal_disc[end])
                        optimal_disc_value[i+1] = continuous[end]
                else
                        optimal_disc_value[i+1] = 0.5 * (continuous[optimal_disc[end][i]]
                                                      +continuous[optimal_disc[end][i]+1])
                end
        end
        #println(smallest_value[end])
        return optimal_disc_value

end

# Equal width discretization
function equal_width_disc(continuous,m)

        N = length(continuous)
        min_value = minimum(continuous)
        max_value = maximum(continuous)
        span = (max_value - min_value)/m
        class_list = Array(Int64,N)

        for i = 1 : N
                class = div((continuous[i] - min_value),span)
                if class >= m
                        class -= 1
                end
                class_list[i] = class + 1
        end
        return class_list
end

function equal_width_edge(continuous,m)
        min_value = continuous[1]
        max_value = continuous[end]
        span = max_value - min_value
        edge = [continuous[1]]
        for i = 1 : m
                edge = [edge, min_value + (span/m)*i]
        end
        return edge
end

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

function cartesian_product(product_set)
        N = length(product_set)
        M = 1
        for i = 1 : N
                M *= length(product_set[i])
        end

        product_returned = [product_set[1]...]

        for class_index = 2 : N
                class_number = length(product_set[class_index])
                current_length = length(product_returned[:,1])

                enlarged_matrix = Array(Int64,current_length*class_number,class_index)

                if class_number == 1
                        enlarged_matrix[:,1:class_index-1] = product_returned
                        for i = 1 : current_length
                                enlarged_matrix[i,class_index] = product_set[class_index][1]
                        end
                else

                        for enlarge_times = 1 : class_number
                                enlarged_matrix[(enlarge_times-1)*current_length+1:enlarge_times*current_length,1:class_index-1] =
                                product_returned
                        end

                        for i = 1 : class_number * current_length
                                item_index = div(i-1,current_length) + 1
                                enlarged_matrix[i,class_index] = product_set[class_index][item_index]
                        end
                end
                product_returned = enlarged_matrix

        end

        return product_returned
end

function find_intval(x,disc)
    if x < disc[2]
        return 1
    end
    if x >= disc[end-1]
        return length(disc)-1
    end

    for i = 2 : length(disc)-2
           if (x >= disc[i])&(x < disc[i+1])
                   return i
           end
    end
end

function continuous_to_discrete(data,bin_edge_1)
        #bin_edge_extend
        bin_edge = copy(bin_edge_1)
        bin_edge[1] = -Inf
        bin_edge[end] = Inf
        data_discrete = Array(Int64,length(data))
        for i = 1 : length(data)
                index = 0
                for j = 2 : length(bin_edge)
                        if (data[i] > bin_edge[j-1])&(data[i] <= bin_edge[j])
                                index = j-1
                        end
                end
                data_discrete[i] = index
        end
        return data_discrete
end

function sort_disc_by_vorder(continuous_order,disc_edge)
      reorder_disc_edge = Array(Any,length(continuous_order))
      for i = 1 : length(continuous_order)
            num_less = 1
            for j = 1 : length(continuous_order)
                  if continuous_order[j] < continuous_order[i]
                        num_less += 1
                  end
            end
            reorder_disc_edge[num_less] = disc_edge[i]
    end
    return reorder_disc_edge
end

function rand_seq(N)

        seq = Array(Int64,N)

        i = 1
        seq[1] = 1 + round(Int64, div(N * rand(),1) )

        while i < N
                number = 1 + round(Int64,div( N*rand(),1))
                if ~(number in seq[1:i])
                        i += 1
                        seq[i] = number
                end
        end
        return seq
end

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

function disc_intval_seq(table)
        N = length(table)
        seq = Array(Int64,N)
        for i = 1 : N
                seq[i] = length(table[i])
        end
        return seq
end

function graph_to_reverse_order(graph)
        order = Array(Int64,length(graph))
        for i = 1 : length(graph)
                order[length(graph) - i + 1] = graph[i][end]
        end
        return order
end

function graph_to_reverse_conti_order(graph,continuous_index)
        order = graph_to_reverse_order(graph)
        Order = Array(Int64,length(continuous_index))
        ind = 0
        for i = 1 : length(order)
                if order[i] in continuous_index
                        ind += 1
                        Order[ind] = order[i]
                end
        end
        return Order
end

function graph_to_markov(graph,target)
        parent_set = []; child_spouse_set = [];
        for i = 1 : length(graph)

                condi_graph = graph[i]
                if length(condi_graph) == 1
                        if target == condi_graph
                                parent_set = []
                        end

                elseif target in condi_graph

                        index = findfirst(condi_graph,target)
                        if index == length(condi_graph)
                                parent_set = [condi_graph[1:end-1]...]
                        else

                                child = condi_graph[end]
                                spouse = []

                                for j = 1 : length(condi_graph)-1
                                        if j != index
                                                spouse = [spouse,condi_graph[j]]
                                        end
                                end
                                child_spouse_set = [child_spouse_set,tuple([child,spouse]...)] ###### change to ; ######

                        end
                else

                end
        end

        for i = 1 : length(child_spouse_set)
                if length(child_spouse_set[i]) == 1
                        child_spouse_set = [child_spouse_set[1:i-1],child_spouse_set[i][1],child_spouse_set[i+1:end]]
                end
        end

        return (parent_set,child_spouse_set)
end

function one_iteration(data,data_integer,graph,discrete_index,continuous_index,lcard,approx = true)

        # Save discretization edge
        disc_edge_collect = Array(Any,length(continuous_index))

        for i = 1 : length(continuous_index)
                target = continuous_index[i]

                increase_order = sortperm(data[:,target])
                conti = data[:,target][increase_order]
                data_integer_sort = Array(Int64,size(data))

                # sort data_integer properly
                for j = 1 : length(data_integer[1,:])
                        data_integer_sort[:,j] = data_integer[:,j][increase_order]
                end
                sets = graph_to_markov(graph,target)
                parent_set = sets[1]; child_spouse_set = sets[2];

                if (length(sets[1]) + length(sets[2])) == 0
                        disc_edge = equal_width_edge(conti,lcard)
                        disc_edge_collect[i] = disc_edge
                else
                        disc_edge = BN_discretizer_free_number_rep(conti,data_integer_sort,parent_set,child_spouse_set,approx)
                        disc_edge_collect[i] = disc_edge
                end

                # Update by current discretization
                data_integer[:,target] = continuous_to_discrete(data[:,target],disc_edge)
        end

        return (data_integer,disc_edge_collect)
end

function BN_discretizer_iteration_converge(data,graph,discrete_index,continuous_index,cut_time,approx = true)
        # intital the first data_integer
        l_card = 0
        for i = 1 : length(discrete_index)
                card = length(unique(data[:,discrete_index[i]]))
                if card > l_card
                        l_card = card
                end
        end

        # pre-equal-width-discretize on continuous variables
        data_integer = zeros(Int64,size(data))
        for i = 1 : length(discrete_index)
                index = discrete_index[i]
                data_integer[:,index] = data[:,index]
        end

        for i = 1 : length(continuous_index)
                index = continuous_index[i]
                data_integer[:,index] = equal_width_disc(data[:,index],l_card)
        end

        disc_edge_collect = Array(Any,length(continuous_index))

        disc_edge_previous = Array(Any,length(continuous_index))

        for i = 1 : length(continuous_index)
                disc_edge_previous[i] = []
                disc_edge_collect[i] = [0]
        end

        times = 0
        while (disc_edge_previous != disc_edge_collect)&(times<cut_time)

            times += 1
            #println(("iteration times = ",times))
            disc_edge_previous = disc_edge_collect

            X = one_iteration(data,data_integer,graph,discrete_index,continuous_index,l_card,approx)
            data_integer = X[1]
            disc_edge_collect = X[2]

        end

        return (data_integer,disc_edge_collect)
end

function K2_f(x,parent)
        N = length(x)
        n = length(parent[1,:])

        # build carti product
        carti_set = [tuple(unique(parent[:,1])...)]

        for i = 2 : n
                carti_set = [carti_set;tuple(unique(parent[:,i])...)]
        end

        uni_x = unique(x); r = length(uni_x)
        carti_set = [carti_set;tuple(uni_x...)]

        CP = cartesian_product(carti_set)

        r_p = length(CP[:,1]) / r
        r_p = round(Int64,r_p)

        # Survey distribution
        count = zeros(Int64,length(CP[:,1]))

        for i = 1 : N
                reform_data = [parent[i,:] x[i]]
                for j = 1 : length(CP[:,1])
                        if reform_data == CP[j,:]
                                count[j] += 1
                        end
                end
        end

        # Sum over same parent value
        count_p = zeros(Int64,r_p)
        for i = 1 : r_p
                count_pp = 0
                for j = 1 : r
                        index = i + r_p * (j-1)
                        count_pp += count[index]
                end
                count_p[i] = count_pp
        end

        # Evaluate
        f = 0
        for i = 1 : r_p
                f += lfact(count_p[i]+r-1) - lfact(r-1)
        end
        for i = 1 : length(CP[:,1])
                f -= lfact(count[i])
        end

        return f
end

function parent_table_to_graph(parent_table)
        m = length(parent_table)
        # Parent_table -> graph
        graph = []
        for i = 1 : m
                child = i
                if length(parent_table[i]) == 0
                        graph = [graph;child]
                else
                parent = []
                for j = 1 : length(parent_table[i])
                        parent = [parent;parent_table[i][j]]
                end

                add = [parent,child]
                graph = [graph;tuple(add...)]
                end
        end
        return graph
end

function K2_one_iteration_discretization(order,u,data_matrix,continuous_index,cut_time,approx = true)
        N = length(data_matrix[:,1])
        m = length(data_matrix[1,:]) # number of variables
        mc = length(continuous_index)

        # continuous/discrete index in the new order
        conti_index = Array(Int64,mc)
        disc_index = Array(Int64,m-mc)
        for i = 1 : mc
                conti_index[i] = findfirst(order,continuous_index[i])
        end
        conti_index = sort(conti_index, rev=true)

        index_disc = 0
        for i = 1 : m
                if ~(i in continuous_index)
                        index_disc += 1
                        disc_index[index_disc] = findfirst(order,i)
                end
        end

        # data that is sorted by order
        data = Array(Any,size(data_matrix))
        for i = 1 : m
                data[:,i] = data_matrix[:,order[i]]
                #println(i)
        end

        # Initialize graph structure
        initial_lcard= 0
        data_discretized = Array(Int64,size(data_matrix))
        if (m - length(conti_index)) == 0
            initial_lcard = round(Int64, log(m)+1)
        else
            for i = 1 : m
                    if ~(i in conti_index)
                            data_discretized[:,i] = data[:,i]
                            current_initial_lcard = length(unique(data[:,i]))

                            if current_initial_lcard > initial_lcard
                                    initial_lcard = current_initial_lcard
                            end
                    end
            end
        end
        for i = 1 : m
                if i in conti_index
                        data_discretized[:,i] = equal_width_disc(data[:,i],initial_lcard)
                end
        end

        # Initialize graph structure
        parent_table = Array(Any,m)
        for i = 1 : m
                parent_table[i] = [] # No parent in the begining
        end
        # store discretization edges
        disc_edge = 0

        # Iterate for all nodes
        for i = 1 : m
                p_old = K2_f(data_discretized[:,i],zeros(Int64,N))
                OKToProceed = true
                current_parent = []

                while OKToProceed & (length(parent_table[i]) < u)
                        # find node z for iteration
                        iteration_list = []
                        for j = 1 : i-1
                                if ~(j in current_parent)
                                        iteration_list = [iteration_list;j]
                                end
                        end

                        # iteration to find most probable
                        current_value = Inf
                        current_best_parent = []

                        for j = 1 : length(iteration_list)

                                iteration_parent = [current_parent;iteration_list[j]]
                                iteration_parent_data = data_discretized[:,iteration_parent[1]]
                                #if i ==4 ; println(iteration_parent); end;

                                for k = 2 : length(iteration_parent)
                                        iteration_parent_data = [iteration_parent_data data_discretized[:,iteration_parent[k]]]
                                end
                                this_iteration_value = K2_f(data_discretized[:,i],iteration_parent_data)

                                if this_iteration_value < current_value
                                        current_value = this_iteration_value
                                        current_best_parent = iteration_parent
                                end
                        end

                        if current_value < p_old
                                p_old = current_value
                                current_parent = current_best_parent
                                parent_table[i] = current_parent
                                ###### Rediscretize ######
                                #println((i,"th varaible"))
                                current_graph = parent_table_to_graph(parent_table)
                                disc_result = BN_discretizer_iteration_converge(data,current_graph,
                                                                               disc_index,conti_index,cut_time,approx)

                                ##### Replace data_discretized #####

                                data_discretized = disc_result[1]
                                disc_edge = disc_result[2]

                        else

                                OKToProceed = false

                        end

                end
                #println(("parent of",i,"is",parent_table[i]))
        end

        score = 0
        # Calculate score
        for i = 1 : m
                if length(parent_table[i]) == 0
                        score += K2_f(data_discretized[:,i],zeros(Int64,N))
                else
                        parent_data_final = Array(Int64,N,length(parent_table[i]))
                        for j = 1 : length(parent_table[i])
                                parent_data_final[:,j] = data_discretized[:,parent_table[i][j]]
                        end
                        score += K2_f(data_discretized[:,i],parent_data_final)
                end
        end

        # Parent_table -> graph
        graph = []
        for i = 1 : m
                child = order[i]
                if length(parent_table[i]) == 0
                        graph = [graph,child]
                else
                parent = []
                for j = 1 : length(parent_table[i])
                        parent = [parent,order[parent_table[i][j]]]
                end

                add = [parent,child]
                graph = [graph,tuple(add...)]
                end
        end

        # Disc_edge is modified to correct order
        conti_index_correspond = Array(Any,mc,3)
        for i = 1 : mc
              conti_index_correspond[i,2] = continuous_index[i]
              conti_index_correspond[i,1] = findfirst(order,continuous_index[i])
        end

        zz = sortperm(conti_index_correspond[:,1], rev = true)
        conti_index_correspond[:,1] = conti_index_correspond[:,1][zz]
        conti_index_correspond[:,2] = conti_index_correspond[:,2][zz]

        for i = 1 : mc
              conti_index_correspond[i,3] = disc_edge[i]
        end

        zzz = sortperm(conti_index_correspond[:,2])
        disc_edge_result = conti_index_correspond[:,3][zzz]

        return (score,graph,disc_edge_result)
end

function K2_w_discretization(data_matrix,u,continuous_index,times,cut_time,approx = true)

        score = Inf
        graph = 0
        disc_edge = 0

        for time = 1 : times
                # Produce random sequence of indexes
                println(("Times of K2 along with discretization method:",time))
                order = rand_seq(length(data_matrix[1,:]))
                iteration_result = K2_one_iteration_discretization(order,u,data_matrix,
                                                          continuous_index,cut_time,approx)

                if iteration_result[1] < score
                        score = iteration_result[1]
                        graph = iteration_result[2]
                        disc_edge = iteration_result[3]
                end
        end

        return (score,graph,disc_edge)
end

function Discretize_All(data_matrix,graph,continuous_index,cut_time)
    discrete_index = Array(Int64,length(data_matrix[1,:]) - length(continuous_index))
    num = 0
    for i = 1 : length(data_matrix[1,:])
        if ~(i in continuous_index)
            num += 1
            discrete_index[num] = i
        end
    end
    sort_continuous = graph_to_reverse_conti_order(graph,continuous_index)
    X = BN_discretizer_iteration_converge(data_matrix,graph,discrete_index,sort_continuous,cut_time,false)
    edge = X[2]
    reorder_edge = sort_disc_by_vorder(sort_continuous,edge)

    return reorder_edge
end

function Learn_DVBN(data_matrix,continuous_index,order,u,cut_time)
    X = K2_one_iteration_discretization(order,u,data_matrix,continuous_index,cut_time,false)
    return (X[2],X[3])
end

function Learn_Discrete_Bayesian_Net(data_matrix,continuous_index,u,cut_time,times)
    X = K2_w_discretization(data_matrix,u,continuous_index,times,cut_time,false)
    return (X[2],X[3],X[1])
end


text:
export iauPvppv
"""
This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a        double[2][3]      first pv-vector
   b        double[2][3]      second pv-vector

Returned:
   apb      double[2][3]      a + b

Note:
   It is permissible to re-use the same array for any of the
   arguments.

Called:
   iauPpp       p-vector plus p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvppv(double a[2][3], double b[2][3], double apb[2][3])
function iauPvppv(a::AbstractMatrix{<:Real}, b::AbstractMatrix{<:Real})
   # Preallocate return values
   amb = zeros(Float64, 3, 2)

   ccall((:iauPvppv, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Ref{Cdouble}),
          convert(Array{Float64, 2}, a'),
          convert(Array{Float64, 2}, b'),
          amb)

   return SMatrix{2,3}(amb')
end


text:
function compute()
	a=BigInt(1)
	n=a<<1000
	str=string(n)
	ans=0;
	for i=1:length(str)
		ans+=str[i]-'0'
	end
	return ans

end

println(compute())


text:
# Dispatches on BVPSystem
function BVPSystem(fun, bc, p, x, M::Integer, order)
    T = eltype(x)
    N = size(x,1)
    y = vector_alloc(T, M, N)
    BVPSystem(order, M, N, fun, bc, p, x, y, vector_alloc(T, M, N), vector_alloc(T, M, N), eltype(y)(undef, M))
end

# If user offers an intial guess
function BVPSystem(fun, bc, p, x, y, order)
    T, U = eltype(x), eltype(y)
    M, N = size(y)
    BVPSystem{T,U}(order, M, N, fun, bc, p, x, y, vector_alloc(T, M, N), vector_alloc(T, M, N), eltype(y)(M))
end

# Auxiliary functions for evaluation
@inline function eval_fun!(S::BVPSystem)
    for i in 1:S.N
        S.fun!(S.f[i], S.y[i], S.p, S.x[i])
    end
end

@inline general_eval_bc_residual!(S::BVPSystem) = S.bc!(S.residual[end], S.y, S.p, S.x)
@inline eval_bc_residual!(S::BVPSystem) = S.bc!(S.residual[end], (S.y[1], S.y[end]), S.p, (S.x[1], S.x[end]))

#=
@inline function banded_update_K!(S::BVPSystem, cache::AbstractMIRKCache, TU::MIRKTableau, i, h)
    M, N, residual, x, y, fun!, order, y_new = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order, S.tmp
    c, v, b, X = TU.c, TU.v, TU.b, TU.x
    ## K, LJ, RJ, Jacobian = cache.K, cache.LJ, cache.RJ, cache.Jacobian
    K = cache.K

    ## index = (M*(i-1)+2):(M*i+1)
    ## Lindex = (M*(i-1)+1):(M*i)
    ## Rindex = (M*(i-1)+1+M):(M*i+M)

    function Kᵣ!(Kr, y, y₁, r)
        x_new = x[i] + c[r]*h
        y_new = (1-v[r])*y + v[r]*y₁
        if r > 1
          y_new += h * sum(j->X[r, j]*K[j], 1:r-1)
        end
        fun!(x_new, y_new, Kr)
    end

    # L is the left strip, and R is the right strip
    # Lᵢ = -I - hᵢ*Σᵣbᵣ*(∂Kᵣ/∂yᵢ)
    # Rᵢ = I - hᵢ*Σᵣbᵣ*(∂Kᵣ/∂y_{i+1})
    # From the paper "A Runge-Kutta Type Boundary Value ODE Solver with Defect Control"
    # by W.H. Enright and Paul Muir
    for r in 1:order
        Kᵣ!(K[r], y[i], y[i+1], r)
        # ∂Kᵣ/∂yᵢ
        ### ForwardDiff.jacobian!(LJ[r], (Kr, y₀)->Kᵣ!(Kr, y₀, y[i+1], r), K[r], y[i])
        # ∂Kᵣ/∂y_{i+1}
        ## ForwardDiff.jacobian!(RJ[r], (Kr, y₁)->Kᵣ!(Kr, y[i], y₁, r), K[r], y[i+1])
        # h*bᵣ*(∂Kᵣ/∂yᵢ)
        ## scale!(-b[r]*h, LJ[r])
        # hᵢ*Σᵣbᵣ*(∂Kᵣ/∂y_{i+1})
        ## scale!(-b[r]*h, RJ[r])
        # sum them up
        ## Jacobian[index,Lindex] += LJ[r]
        ## Jacobian[index,Rindex] += RJ[r]
        # fun_jac!(LJ[r], fun!, x_new, y_new, K[r])
        # fun_jac!(RJ[r], fun!, x_new, y_new, K[r])
    end
    # Lᵢ = -I - ...
    # Rᵢ = I - ...
    ## Jacobian[index,Lindex] -= I
    ## Jacobian[index,Rindex] += I
end

function banded_Φ!(S::BVPSystem, TU::MIRKTableau, cache::AbstractMIRKCache)
    order, residual, N, y, x = S.order, S.residual, S.N, S.y, S.x
    K, b = cache.K, TU.b
    for i in 1:N-1
        h = x[i+1] - x[i]
        banded_update_K!(S, cache, TU, i, h)
        # Update residual
        residual[i] = y[i+1] - y[i] - h * sum(j->b[j]*K[j], 1:order)
    end
    eval_bc_residual!(S)
    #ForwardDiff.jacobian!(@view(cache.Jacobian[1:S.M, 1:S.M]),                     (x,y)->S.bc!(x,y,S.y[1]),   residual[1],   S.y[1])
    #ForwardDiff.jacobian!(@view(cache.Jacobian[(end-S.M+1):end, (end-S.M+1):end]), (x,y)->S.bc!(x,S.y[end],y), residual[end], S.y[end])
    #display(cache.Jacobian)
end
=#

@inline function update_K!(S::BVPSystem, cache::AbstractMIRKCache, TU::MIRKTableau, i, h)
    M, N, residual, x, y, fun!, order = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order
    K, b = cache.K, TU.b
    c, v, X = TU.c, TU.v, TU.x

    function Kᵣ!(Kr, y, y₁, r)
        x_new = x[i] + c[r]*h
        y_new = (1-v[r])*y + v[r]*y₁
        if r > 1
          y_new += h * sum(j->X[r, j]*K[j], 1:r-1)
        end
        fun!(Kr, y_new, S.p, x_new)
    end

    for r in 1:order
        Kᵣ!(K[r], y[i], y[i+1], r)
    end
end

function Φ!(S::BVPSystem{T}, TU::MIRKTableau, cache::AbstractMIRKCache) where T
    M, N, residual, x, y, fun!, order = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order
    K, b = cache.K, TU.b
    c, v, X = TU.c, TU.v, TU.x
    for i in 1:N-1
        h = x[i+1] - x[i]
        # Update K
        update_K!(S, cache, TU, i, h)
        # Update residual
        residual[i] = y[i+1] - y[i] - h * sum(j->b[j]*K[j], 1:order)
    end
end


text:
module TestDaemonTask

using Base.Test
using SimpleTasks.Types
using SimpleTasksTests.Utils.MockServices

import SimpleTasks.Tasks.DaemonTask
import SimpleTasks.Tasks.DaemonTask

type NewTask <: DaemonTaskDetails end
function test_execute_undefined_task()
    task = NewTask()
    datasource = MockDatasourceService()
    @test_throws ErrorException DaemonTask.execute(task, datasource)
end

function __init__()
    test_execute_undefined_task()
end

end # module TestDaemonTask

max_forks_repo_name: seung-lab/SimpleTasks


text:
include("utils.jl")

using Manifolds: default_metric_dispatch

@testset "Symmetric Positive Definite Matrices" begin
    M1 = SymmetricPositiveDefinite(3)
    @test repr(M1) == "SymmetricPositiveDefinite(3)"
    M2 = MetricManifold(SymmetricPositiveDefinite(3), Manifolds.LinearAffineMetric())
    M3 = MetricManifold(SymmetricPositiveDefinite(3), Manifolds.LogCholeskyMetric())
    M4 = MetricManifold(SymmetricPositiveDefinite(3), Manifolds.LogEuclideanMetric())

    @test (@inferred default_metric_dispatch(M2)) === Val(true)
    @test (@inferred default_metric_dispatch(M1, Manifolds.LinearAffineMetric())) ===
          Val(true)
    @test (@inferred default_metric_dispatch(M1, Manifolds.LogCholeskyMetric())) ===
          Val(false)
    @test (@inferred default_metric_dispatch(M3)) === Val(false)
    @test is_default_metric(M2)
    @test is_default_metric(M1, Manifolds.LinearAffineMetric())
    @test !is_default_metric(M1, Manifolds.LogCholeskyMetric())
    @test !is_default_metric(M3)

    @test injectivity_radius(M1) == Inf
    @test injectivity_radius(M1, one(zeros(3, 3))) == Inf
    @test injectivity_radius(M1, ExponentialRetraction()) == Inf
    @test injectivity_radius(M1, one(zeros(3, 3)), ExponentialRetraction()) == Inf
    @test zero_vector(M1, one(zeros(3, 3))) == zero_vector(M2, one(zeros(3, 3)))
    @test zero_vector(M1, one(zeros(3, 3))) == zero_vector(M3, one(zeros(3, 3)))
    metrics = [M1, M2, M3]
    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{3,3,Float64,9})

    for M in metrics
        basis_types = if M == M3
            ()
        else
            (DefaultOrthonormalBasis(),)
        end
        @testset "$(typeof(M))" begin
            @test representation_size(M) == (3, 3)
            for T in types
                exp_log_atol_multiplier = 8.0
                if T <: MMatrix{3,3,Float64}
                    # eigendecomposition of 3x3 SPD matrices from StaticArrays is not very accurate
                    exp_log_atol_multiplier = 5.0e7
                end
                if M == M3 && T <: MMatrix
                    # Cholesky or something does not work in vector_transport yet for MMatrix
                    continue
                end
                A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)]
                ptsF = [#
                    [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1],
                    [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1],
                    A(π / 6) * [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] * transpose(A(π / 6)),
                ]
                pts = [convert(T, a) for a in ptsF]
                test_manifold(
                    M,
                    pts;
                    test_default_vector_transport=true,
                    vector_transport_methods=typeof(M) == SymmetricPositiveDefinite{3} ?
                                             [ParallelTransport()] : [],
                    test_forward_diff=false,
                    test_reverse_diff=false,
                    test_vee_hat=M === M2,
                    exp_log_atol_multiplier=exp_log_atol_multiplier,
                    basis_types_vecs=basis_types,
                    basis_types_to_from=basis_types,
                    is_tangent_atol_multiplier=1,
                )
            end
            @testset "Test Error cases in is_point and is_vector" begin
                pt1f = zeros(2, 3) # wrong size
                pt2f = [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 1.0] # not positive Definite
                pt3f = [2.0 0.0 1.0; 0.0 1.0 0.0; 0.0 0.0 4.0] # not symmetric
                pt4 = [2.0 1.0 0.0; 1.0 2.0 0.0; 0.0 0.0 4.0]
                @test !is_point(M, pt1f)
                @test !is_point(M, pt2f)
                @test !is_point(M, pt3f)
                @test is_point(M, pt4)
                @test !is_vector(M, pt4, pt1f)
                @test is_vector(M, pt4, pt2f)
                @test !is_vector(M, pt4, pt3f)
            end
        end
    end
    x = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
    y = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
    @testset "Convert SPD to Cholesky" begin
        v = log(M1, x, y)
        (l, w) = Manifolds.spd_to_cholesky(x, v)
        (xs, vs) = Manifolds.cholesky_to_spd(l, w)
        @test isapprox(xs, x)
        @test isapprox(vs, v)
    end
    @testset "Preliminary tests for LogEuclidean" begin
        @test representation_size(M4) == (3, 3)
        @test isapprox(distance(M4, x, y), sqrt(2) * log(2))
        @test manifold_dimension(M4) == manifold_dimension(M1)
    end
    @testset "Test for tangent ONB on LinearAffineMetric" begin
        v = log(M2, x, y)
        donb = get_basis(base_manifold(M2), x, DiagonalizingOrthonormalBasis(v))
        X = get_vectors(base_manifold(M2), x, donb)
        k = donb.data.eigenvalues
        @test isapprox(0.0, first(k))
        for i in 1:length(X)
            @test isapprox(1.0, norm(M2, x, X[i]))
            for j in (i + 1):length(X)
                @test isapprox(0.0, inner(M2, x, X[i], X[j]))
            end
        end
        d2onb = get_basis(M2, x, DiagonalizingOrthonormalBasis(v))
        @test donb.data.eigenvalues == d2onb.data.eigenvalues
        @test get_vectors(base_manifold(M2), x, donb) == get_vectors(M2, x, d2onb)
    end
    @testset "Vector transport and transport along with Schild and Pole ladder" begin
        A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)]
        M = SymmetricPositiveDefinite(3)
        p1 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
        p2 = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
        p3 = A(π / 6) * [1.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] * transpose(A(π / 6))
        X1 = log(M, p1, p3)
        Y1 = vector_transport_to(M, p1, X1, p2)
        @test is_vector(M, p2, Y1)
        Y2 = vector_transport_to(M, p1, X1, p2, PoleLadderTransport())
        @test is_vector(M, p2, Y2, atol=10^-16)
        Y3 = vector_transport_to(M, p1, X1, p2, SchildsLadderTransport())
        @test is_vector(M, p2, Y3)
        @test isapprox(M, p2, Y1, Y2) # pole is exact on SPDs, i.e. identical to parallel transport
        @test norm(M, p1, X1) ≈ norm(M, p2, Y1) # parallel transport is length preserving
        # test isometry
        X2 = log(M, p1, p2)
        Y4 = vector_transport_to(M, p1, X2, p2)
        @test norm(M, p1, X2) ≈ norm(M, p2, Y4)
        @test is_vector(M, p2, Y4)
        Y5 = vector_transport_to(M, p1, X2, p2, PoleLadderTransport())
        @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y1, Y4) # parallel transport isometric
        @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y2, Y5) # pole ladder transport isometric
    end
end


text:
## Example for using LSTNet

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg;
Pkg.activate(".");
Pkg.instantiate()

@info "Loading packages"
using FluxArchitectures
using Plots

# Load some sample data
@info "Loading data"
poollength = 10
horizon = 15
datalength = 1000
input, target = get_data(:exchange_rate, poollength, datalength, horizon) |> gpu

# Define the network architecture
@info "Creating model and loss"
inputsize = size(input, 1)
convlayersize = 2
recurlayersize = 3
skiplength = 120

# Define the neural net
model = LSTnet(inputsize, convlayersize, recurlayersize, poollength, skiplength,
    init = Flux.zeros32, initW = Flux.zeros32) |> gpu

# MSE loss
function loss(x, y)
    Flux.reset!(model)
    return Flux.mse(model(x), permutedims(y))
end

# Callback for plotting the training
cb = function ()
    Flux.reset!(model)
    pred = model(input) |> permutedims |> cpu
    Flux.reset!(model)
    p1 = plot(pred, label = "Predict")
    p1 = plot!(cpu(target), label = "Data", title = "Loss $(loss(input, target))")
    display(plot(p1))
end

# Training loop
@info "Start loss" loss = loss(input, target)
@info "Starting training"
Flux.train!(loss, Flux.params(model), Iterators.repeated((input, target), 20),
    ADAM(0.01), cb = cb)

@info "Finished"
@info "Final loss" loss = loss(input, target)


text:
using BoundingBoxes
using Test

b_ref = BoundingBox([1.0, 1.0], [2.0, 2.0])
b_inside = BoundingBox([1.2, 1.2], [1.8, 1.8])
b_overlap1 = BoundingBox([1.8, 1.8], [2.2, 2.2])
b_overlap2 = BoundingBox([0.8, 0.8], [1.2, 1.2])
b_overlap3 = BoundingBox([1.8, 0.8], [2.2, 1.2])
b_overlap4 = BoundingBox([0.8, 1.8], [1.2, 2.2])
b_outside_all1 = BoundingBox([3.0, 3.0], [4.0, 4.0])
b_outside_all2 = BoundingBox([-4.0, -4.0], [-3.0, -3.0])
b_outside_all3 = BoundingBox([3.0, -4.0], [4.0, -3.0])
b_outside_all4 = BoundingBox([-4.0, 3.0], [-3.0, 4.0])
b_outside_partial1 = BoundingBox([1.1, 3.0], [2.1, 4.0])
b_outside_partial2 = BoundingBox([3.0, 1.1], [4.0, 2.1])
b_outside_partial3 = BoundingBox([1.2, 3.0], [1.8, 4.0])
b_outside_partial4 = BoundingBox([3.0, 1.2], [4.0, 1.8])
b_outside_inf = BoundingBox([3.0, 3.0], [Inf, Inf])
b_overlap_inf = BoundingBox([1.8, 1.8], [Inf, Inf])
b_overlap_inf2 = BoundingBox([-Inf, -Inf], [Inf, Inf])
b_empty = BoundingBox([NaN, NaN], [NaN, NaN])

b_true = [b_inside, b_overlap1, b_overlap2, b_overlap3, b_overlap4, b_overlap_inf, b_overlap_inf2]
b_false = [b_outside_all1, b_outside_all2, b_outside_all3, b_outside_all4, b_outside_partial1,
  b_outside_partial2, b_outside_partial3, b_outside_partial4, b_outside_inf, b_empty]

function overlap_sym(a, b, answer)
  @test overlaps(a, b)==answer
  @test overlaps(b, a)==answer
end

for b in b_true
  overlap_sym(b_ref, b, true)
end

for b in b_false
  overlap_sym(b_ref, b, false)
end

for b in vcat(b_true, b_false)
  overlap_sym(b_empty, b, false)
end

max_stars_repo_name: goretkin/BoundingBoxes

max_issues_repo_name: goretkin/BoundingBoxes

max_forks_repo_name: goretkin/BoundingBoxes


text:
export Model, SyntheticData

"""
        Model

A structure to perform the twin-experiment

## Fields
- `Nx` : Dimension of the state variable
- `Ny` : Dimension of the observation variable
- `Δtdyn` : Time-step for the dynamical model
- `Δtobs` : Time step between two observations of the state
- `ϵx` : Process noise
- `ϵy` : Observation noise
- `π0` : Multivariate distribution for the initial condition
- `Tburn` : Number of steps to burn from the end of the spin up to compute the metrics
- `Tstep` : Number of steps for which the filter to be tested is applied
- `Tspinup` : Number of steps of spin-up phase,
   i.e. number of steps to generate the initial ensemble for the filtering algorithms
- `F` : State-Space Model
"""
struct Model

    "Dimension of the state variable"
    Nx::Int64

	"Dimension of the observation variable"
    Ny::Int64

	"Time-step for the dynamical model"
    Δtdyn::Float64

	"Time step between two observations of the state"
    Δtobs::Float64

	"Process noise"
    ϵx::InflationType

	"Observation noise"
    ϵy::AdditiveInflation

    "Multivariate distribution for the initial condition"
    π0::ContinuousMultivariateDistribution

    "Number of steps to burn from the end of the spin up to compute the metrics"
    Tburn::Int64
    Tstep::Int64
    Tspinup::Int64

    " State-Space Model"
    F::StateSpace
end

"""
    SyntheticData

A structure to store the synthetic data in a twin-experiment

## Fields
- `tt` : time history
- `Δt` : time step
- `x0` : the initial condition
- `xt` : history of the state
- `yt` : history of the observations
"""
struct SyntheticData
	tt::Array{Float64,1}
	Δt::Float64
	x0::Array{Float64,1}
	xt::Array{Float64,2}
	yt::Array{Float64,2}
end


text:
"""
    se_bootstrap(semfit::SemFit; n_boot = 3000, data = nothing, kwargs...)

Return boorstrap standard errors.
Only works for single models.

# Arguments
- `n_boot`: number of boostrap samples
- `data`: data to sample from. Only needed if different than the data from `sem_fit`
- `kwargs...`: passed down to `swap_observed`
"""
function se_bootstrap(semfit::SemFit; n_boot = 3000, data = nothing, specification = nothing, kwargs...)

    if isnothing(data)
        data = get_data(observed(model(semfit)))
    end

    data = prepare_data_bootstrap(data)

    start = solution(semfit)

    new_solution = zero(start)
    sum = zero(start)
    squared_sum = zero(start)

    n_failed = 0.0

    converged = true

    for _ = 1:n_boot

        sample_data = bootstrap_sample(data)
        new_model = swap_observed(
            model(semfit);
            data = sample_data,
            specification = specification,
            kwargs...
        )

        new_solution .= 0.0

        try
            new_solution = solution(sem_fit(new_model; start_val = start))
        catch
            n_failed += 1
        end

        @. sum += new_solution
        @. squared_sum += new_solution^2

        converged = true
    end

    n_conv = n_boot - n_failed
    sd = sqrt.(squared_sum/n_conv - (sum/n_conv).^2)
    print("Number of nonconverged models: ", n_failed, "\n")
    return sd

end

function prepare_data_bootstrap(data)
    return Matrix(data)
end

function bootstrap_sample(data)
    nobs = size(data, 1)
    index_new = rand(1:nobs, nobs)
    data_new = data[index_new, :]
    return data_new
end


text:
module WindAttackTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Snowballs" => Ahorn.EntityPlacement(
        Maple.WindAttackTrigger,
        "rectangle",
    ),
)

end


text:
using Uhlarm
using Base.Test


text:
using LinearAlgebra, OrdinaryDiffEq, Test, PreallocationTools, ForwardDiff, GalacticOptim, Optim

randmat = rand(5, 3)
sto = similar(randmat)
function claytonsample!(sto, τ, α; randmat=randmat)
    sto = get_tmp(sto, τ)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 1] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)*α
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

#= taking the second derivative of claytonsample! with respect to τ with manual chunk_sizes.
In setting up the dualcache, we are setting chunk_size to [1, 1], because we differentiate
only with respect to τ. This initializes the cache with the minimum memory needed. =#
stod = dualcache(sto, [1, 1])
df3 = ForwardDiff.derivative(τ -> ForwardDiff.derivative(ξ -> claytonsample!(stod, ξ, 0.0), τ), 0.3)

#= taking the second derivative of claytonsample! with respect to τ with auto-detected chunk-size.
For the given size of sto, ForwardDiff's heuristic chooses chunk_size = 8. Since this is greater
than what's needed (1+1), the auto-allocated cache is big enough to handle the nested dual numbers, even
if we don't specify the keyword argument levels = 2. This should in general not be relied on to work,
especially if more levels of nesting occur (see optimization example below). =#
stod = dualcache(sto)
df4 = ForwardDiff.derivative(τ -> ForwardDiff.derivative(ξ -> claytonsample!(stod, ξ, 0.0), τ), 0.3)

@test df3 ≈ df4

#= taking the second derivative of claytonsample! with respect to τ with auto-detected chunk-size.
For the given size of sto, ForwardDiff's heuristic chooses chunk_size = 8 and with keyword arg levels = 2,
the created cache size is larger than what's needed (even more so than the last example). =#
stod = dualcache(sto, levels = 2)
df5 = ForwardDiff.derivative(τ -> ForwardDiff.derivative(ξ -> claytonsample!(stod, ξ, 0.0), τ), 0.3)

@test df3 ≈ df5

#= Checking nested dual numbers using optimization problem involving Optim.jl's Newton() (involving Hessians);
so, we will need one level of AD for the ODE solver (TRBDF2) and two more to calculate the Hessian =#
function foo(du, u, p, t)
    tmp = p[2]
    A = reshape(p[1], size(tmp.du))
    tmp = get_tmp(tmp, u)
    mul!(tmp, A, u)
    @. du = u + tmp
    nothing
end

ps = 2 #use to specify problem size; don't go crazy on this, because of the compilation time...
coeffs = -collect(0.1:0.1:(ps^2/10))
cache = dualcache(zeros(ps,ps), levels = 3)
prob = ODEProblem(foo, ones(ps, ps), (0., 1.0), (coeffs, cache))
realsol = solve(prob, TRBDF2(), saveat = 0.0:0.1:10.0, reltol = 1e-8)

function objfun(x, prob, realsol, cache)
    prob = remake(prob, u0 = eltype(x).(prob.u0), p = (x, cache))
    sol = solve(prob, TRBDF2(), saveat = 0.0:0.1:10.0, reltol = 1e-8)

    ofv = 0.0
    if any((s.retcode != :Success for s in sol))
        ofv = 1e12
    else
        ofv = sum((sol.-realsol).^2)
    end
    return ofv
end
fn(x,p) = objfun(x, p[1], p[2], p[3])
optfun = OptimizationFunction(fn, GalacticOptim.AutoForwardDiff())
optprob = OptimizationProblem(optfun, zeros(length(coeffs)), (prob, realsol, cache))
newtonsol = solve(optprob, Newton())

@test all(abs.(coeffs .- newtonsol.u) .< 1e-3)

#an example where chunk_sizes are not the same on all differentiation levels:
cache = dualcache(zeros(ps,ps), [4, 4, 2])
prob = ODEProblem(foo, ones(ps, ps), (0., 1.0), (coeffs, cache))
realsol = solve(prob, TRBDF2(chunk_size = 2), saveat = 0.0:0.1:10.0, reltol = 1e-8)

function objfun(x, prob, realsol, cache)
    prob = remake(prob, u0 = eltype(x).(prob.u0), p = (x, cache))
    sol = solve(prob, TRBDF2(chunk_size = 2), saveat = 0.0:0.1:10.0, reltol = 1e-8)

    ofv = 0.0
    if any((s.retcode != :Success for s in sol))
        ofv = 1e12
    else
        ofv = sum((sol.-realsol).^2)
    end
    return ofv
end

fn(x,p) = objfun(x, p[1], p[2], p[3])

optfun = OptimizationFunction(fn, GalacticOptim.AutoForwardDiff())
optprob = OptimizationProblem(optfun, zeros(length(coeffs)), (prob, realsol, cache))
newtonsol2 = solve(optprob, Newton())

@test all(abs.(coeffs .- newtonsol2.u) .< 1e-3)


text:
module UniDerivative
"""
Compute derivative of a univariate function.
"""

abstract type DerivativeMethod end

struct CentralDiff <: DerivativeMethod
    h
    CentralDiff(h = cbrt(eps())) = new(h)
end

function df(param::CentralDiff, f, x::Real)::Real
    (f(x + param.h / 2) - f(x - param.h / 2)) / param.h
end

struct ComplexDiff <: DerivativeMethod
    h
    ComplexDiff(h = 1e-20) = new(h)
end

function df(param::ComplexDiff, f, x::Real)::Real
    imag(f(x + param.h * im)) / param.h
end

end # module


text:
import Base: writemime

# HTML utils

using .DOM

fade(x) = span(".fade", x)

# Text

function displayinline(t::Text)
  lines = split(string(t), "\n")
  Collapsible(span(".text", lines[1]),
              span(".text", join(lines[2:end], "\n")))
end

# Tables

type Table{T}
  class::ASCIIString
  data::AbstractMatrix{T}
end

Table(data::AbstractMatrix) = Table("", data)

const MAX_CELLS = 500

function getsize(h, w, maxcells)
  (h == 0 || w == 0) && return 0, 0
  swap = false
  w > h && ((w, h, swap) = (h, w, true))
  h = min(maxcells ÷ w, h)
  h ≥ w ?
    (swap ? (w, h) : (h, w)) :
    (ifloor(sqrt(maxcells)), ifloor(sqrt(maxcells)))
end

function Base.writemime(io::IO, m::MIME"text/html", table::Table)
  println(io, """
""")
  h, w = size(table.data)
  (h == 0 || w == 0) && @goto none
  h′, w′ = getsize(h, w, MAX_CELLS)
  for i = (h′ == h ? (1:h′) : [1:(h′÷2), h-(h′÷2)+1:h])
    println(io, """
""")
    for j = (w′ == w ? (1:w′) : [1:(w′÷2), w-(w′÷2)+1:w])
      println(io, """""")
      if isdefined(table.data, i, j)
        item = applydisplayinline(table.data[i, j])
        writemime(io, bestmime(item), item)
      else
        print(io, "#undef")
      end
      println(io, """	""")

      w > w′ && j == (w′÷2) && println(io, """⋯	""")
    end
    println(io, """
""")

    h > h′ && i == (h′÷2) && println(io, "
","⋮	"^(w≤w′?w:w′+1),"
")
  end
  @label none
  println(io, """
""")
end

# Void

displayinline(::Void) = Text("✓")

# Floats

function round3(n)
  n = string(n)
  n = replace(n, r"\.0$", ".")
  ismatch(r"[0-9]\.0*999", n) && return n
  zero = ismatch(r"^[^\d]*0\.0+", n) # Special case for e.g. 0.0001
  r = zero ? r"0[1-9][0-9]{3,}[^\.]*$" : r"\.[0-9]{4,}"
  n = replace(n, r, s->string(s[1], @sprintf("%03d", round(Integer, parse(Int, s[2:5])/10 ))), 1)
end

function writemime(io::IO, m::MIME"text/html", x::AbstractFloat)
  print(io, """""")
  print(io, round3(x))
  print(io, """""")
end

displayinline!(x::AbstractFloat, opts) =
  showresult(stringmime("text/html", x), opts, html=true)

# Functions

name(f::Function) =
  isgeneric(f) ? string(f.env.name) :
  isdefined(f, :env) && isa(f.env,Symbol) ? string(f.env) :
  "λ"

displayinline(f::Function) =
  isgeneric(f) ?
    Collapsible(strong(name(f)), methods(f)) :
    Text(name(f))

# Arrays

sizestr(a::AbstractArray) = join(size(a), "×")

displayinline(a::Matrix) =
  Collapsible(span(strong("Matrix "), fade("$(eltype(a)), $(sizestr(a))")),
              Table("array", a))

function copytranspose(xs::Vector)
  result = similar(xs, length(xs), 1)
  for i = 1:length(xs)
    isdefined(xs, i) && (result[i] = xs[i])
  end
  return result
end

displayinline(a::Vector, t = "Vector") =
  Collapsible(span(strong(t), fade(" $(eltype(a)), $(length(a))")),
              Table("array", copytranspose(a)))

displayinline(s::Set) = displayinline(collect(s), "Set")

@cond if VERSION < v"0.4-"
  keytype(d) = eltype(d)[1]
  valtype(d) = eltype(d)[2]
else
  keytype(d) = eltype(d).parameters[1]
  valtype(d) = eltype(d).parameters[2]
end

displayinline(d::Dict) =
  Collapsible(span(strong("Dictionary "), fade("$(keytype(d)) → $(valtype(d)), $(length(d))")),
              HTML() do io
                println(io, """
""")
                kv = collect(d)
                for i = 1:(min(length(kv), MAX_CELLS÷2))
                  print(io, "
")
                  item = displayinline(kv[i][1])
                  writemime(io, bestmime(item), item)
                  print(io, "	")
                  item = displayinline(kv[i][2])
                  writemime(io, bestmime(item), item)
                  print(io, "
")
                end
                length(kv) ≥ MAX_CELLS÷2 && println(io, """⋮	⋮	""")
                println(io, """
""")
              end)

# Others

import Jewel: @require

# Data Frames

@require DataFrames begin
  displayinline(f::DataFrames.DataFrame) =
    isempty(f) ? Collapsible(span(strong("DataFrame "), fade("Empty"))) :
      Collapsible(span(strong("DataFrame "), fade("($(join(names(f), ", "))), $(size(f,1))")),
                  Table("data-frame", vcat(map(s->HTML(string(s)), names(f))',
                                           DataFrames.DataArray(f))))
end

# Colors

@require Color begin
  displayinline(c::Color.ColourValue) =
    Collapsible(span(strong(@d(:style => "color: #$(Color.hex(c))"),
                            "#$(Color.hex(c)) "),
                     fade(string(c))),
                tohtml(MIME"image/svg+xml"(), c))

  displayinline{C<:Color.ColourValue}(cs::VecOrMat{C}) = tohtml(MIME"image/svg+xml"(), cs)
end

# Gadfly

@require Gadfly begin
  displayinline(p::Gadfly.Plot) = DOM.div(p, style = "background: white")
end

# PyPlot

@require PyPlot begin
  try
    PyPlot.pygui(true)
  catch e
    warn("PyPlot is set to display in the console")
  end
end

# Images

@require Images begin
  displayinline{T,N,A}(img::Images.Image{T,N,A}) =
    Collapsible(HTML("""$(strong("Image")) $(sizestr(img)), $T"""),
                HTML(applydisplayinline(img.properties),tohtml(MIME"image/png"(), img)))
end

# Expressions

fixsyms(x) = x
fixsyms(x::Symbol) = @> x string replace(r"#", "_") symbol
fixsyms(ex::Expr) = Expr(ex.head, map(fixsyms, ex.args)...)

function displayinline(x::Expr)
  rep = stringmime(MIME"text/plain"(), x |> fixsyms)
  lines = split(rep, "\n")
  html = span(".code.text", @d("data-lang" => "julia2"), rep)
  length(lines) == 1 && length(lines[1]) ≤ 50 ?
    Collapsible(html) :
    Collapsible(strong("Julia Code"), html)
end

# Profile tree

function toabspath(file)
  isabspath(file) && file
  path = basepath(file)
  return path == nothing ? file : path
end

@require Jewel.ProfileView begin
  function displayinline!(tree::Jewel.ProfileView.ProfileTree, opts)
    raise(opts[:editor], "julia.profile-result",
          @d("value" => stringmime("text/html", tree),
             "bounds" => opts[:bounds],
             "lines" => [@d(:file => toabspath(li.file),
                            :line => li.line,
                            :percent => p) for (li, p) in Jewel.ProfileView.fetch() |> Jewel.ProfileView.flatlines]))
  end
end


text:
# First order explicit upwind scheme for the non-linear Shallow-Water Equations
import Pkg
ENV["PLOTS_USE_ATOM_PLOTPANE"] = "false"
Pkg.add("Plots")
using Plots

function thomas_alg(a,b,c,d)
    #
    # Thomas algorithm
    # Name of the function: thomas_alg
    # Name of the file:     thomas_alg.m
    # Arguments:
    #   a = vector of lower diagonal elements
    #   b = vector of diagonal elements
    #   c = vector of upper diagonal elements
    #   d = known right hand side vector
    # Results / output:
    #   x = solution of linear system
    #

    c[1]=c[1]/b[1];
    d[1]=d[1]/b[1];
    #
    N = length(d);      # number of elements in the vector
    #
    x = zeros(Float64,N);
    #
    # Part I: Forward elimination
    for i = 2:N
        K    = 1/(b[i]-c[i-1]*a[i]);
        c[i] = c[i]*K;
        d[i] = (d[i]-a[i]*d[i-1])*K;
    end
    #
    # Part II: Backward elimination
    x[N] = d[N];
    #
    for i = N-1:-1:1
        x[i] = d[i]-c[i]*x[i+1];
    end
    return x
end

function IntTrajectory(xstart, x, dx, dt, q, BCL, BCR)
    # Name of the fuction is IntTrajectory
    # Name of the file (mandatory) is IntTrajectory.m
    # (Can also be a vector function for all points in the domain together, using as vectors xstart and q)
    # Input / arguments:
    #   xstart = starting point of integration trajectory (scalar)
    #   x      = vector of spatial coordinates
    #   dx     = mesh spacing (for the equidistant grid)
    #   dt     = final time of the ODE integrator = physical time
    #   q      = vector of state variables
    #   BCL    = boundary condition on the left
    #   BCR    = boundary condition on the right
    #
    CFL  = 0.1;          # local Courant-Friedrichs-Levy number
    xL   = x[1];         # location of left boundary
    xR   = x[end];       # location of the right boundary
    tau  = 0;            # initialize local time tau
    x0   = xstart;       # starting point of the characteristic
    NMAX = 10000;        # max. number of integration steps
    #
    q0 = 0.0
    #
    for n = 1:NMAX
        if (x0 <= xL)
            q0 = BCL;
            break
        end
        if (x0 >= xR)
            q0 = BCR;
            break
        end
        j    = ceil(Int64,(x0-xL)/dx);                    # current nearest index j associated to point x0
        p    = (x0-x[j])/dx;                          # remainder (real number 0 <= p < 1)
        q0   = (1-p)*q[j] + p*q[j+1];               # linear interpolation of q at x0
        qmax = max( abs(q[j]),abs(q[j+1]) );        # max. local velocity
        dtau = CFL*dx/qmax;                         # local time step
        if (tau + dtau > dt)
           dtau = dt - tau;                         # match final time dt EXACTLY
        end
        if (tau >= dt)
            break
        end
        #x1   = x0 - dtau*q0;                                     # update x0 with first order explicit Euler method
        x1   = x0 - dtau*q0 + 0.5*dtau^2*(q[j+1]-q[j])/dx*q0;    # update x0 with second order explicit Euler method
        x0   = x1;
        tau  = tau + dtau;                          # update the local time
    end
    return q0
end

function timeloop(NMAX, times, time0, dx, grid0, grid, state, bc, physics, checking)
    # Loop in time

    # Unpack the time constants
    dt0 = time0["dt0"]
    tend = time0["tend"]

    # Unpack the grid constants
    imax = grid0["imax"]
    IMAX = grid0["IMAX"]

    # Unpack the grid arrays
    x = grid["x"]
    xb = grid["xb"]

    # Unpack the state
    h = state["h"]
    u = state["u"]
    hu = state["hu"]
    hi = state["hi"]
    eta = state["eta"]
    etaj = state["etaj"]

    # Unpack the boundary conditions
    BCR = bc["BCR"]
    BCL = bc["BCL"]

    # Unpack the physics
    g = physics["grav"]

    # Unpack checking variables
    masstot = checking["masstot"]
    momtot  = checking["momtot"]

    b = Array{Float64,1}(undef, imax)
    d = Array{Float64,1}(undef, imax)
    ld = Array{Float64,1}(undef, imax)
    ud = Array{Float64,1}(undef, imax)

    eta1 = Array{Float64,1}(undef, imax)
    u1 = Array{Float64,1}(undef, IMAX)
    Fu = Array{Float64,1}(undef, IMAX)

    for n = 1:NMAX
    #    amax = max(abs(q))          # Maximum wave propagation speed
    #    dt = C*dx/amax              # Time step
        dt = min(dt0, tend-times[n])
        hu[1] = max(0, h[1]+eta[1])
        #H(1) = 0
        for i = 2:IMAX-1
            hu[i] = max(0, h[i]+eta[i-1],h[i]+eta[i])
        end
        hu[IMAX] = max(0, h[imax]+eta[imax])
        #H(IMAX) = 0
        #Htot(n) = sum(H)
        if (times[n] >= tend)
        #    Hdiff = Htot(n) - Htot(1)
            break
        end

        for i = 1:IMAX
            if (hu[i] <= 0.0)
                u[i] = 0.0                                  # empty (dry)interfaces do not have a velocity
            end
        end

        # Loop in space over all velocity points
        for i = 1:IMAX
            Fu[i] = IntTrajectory(x[i], x, dx, dt, u, BCL, BCR)
        end

        for i = 1:imax
            d[i]  = eta[i] - dt/dx*( hu[i+1]*Fu[i+1] - hu[i]*Fu[i] )
            ld[i] = -g*dt^2/dx^2*hu[i]
            ud[i] = -g*dt^2/dx^2*hu[i+1]
        end
        ld[1]    = 0
        ud[imax] = 0
        b = ones(Float64, imax) - ld - ud
        eta1 = thomas_alg(ld', b', ud', d')'
        u1[1] = BCL
        for i = 2:IMAX-1
            if (hu[i] > 0.01)
                u1[i]   = Fu[i] - g*dt/dx*(eta1[i]-eta1[i-1])
            else
                u1[i] = 0
            end
        end
        u1[IMAX] = BCR
        #
        times[n+1] = times[n] + dt
        for i = 1:imax
            eta[i] = max(eta1[i], -hi[i])
        end   # overwrite old solution with new one
        u  = u1
        #masstot[n+1]   = sum(max(eta1, -hi) + hi)
        masstot[n+1]   = sum(eta + hi)
        etaj[1] = eta[1]
        for i = 2:IMAX-1
            if (i > 1 && i < IMAX)
                etaj[i] = 0.5*(eta[i-1]+eta[i])
            end
        end
        etaj[IMAX]  = eta[imax]
        for i = 1:IMAX
            hu[i] = max(0.0,h[i] + etaj[i])
        end
        momtot[n+1] = sum(hu.*u)
        # Plot
        #subplot(4,1,1)
        plot(xb,[eta, -hi])
        filnam = string("wl\_t=",string(times[n+1]),".png")
        png(filnam)
        #title(sprintf('Current time t = #f',times[n]))
        #title(sprintf('Water level'))
        #axis([xL,xR,-2,2])
        #title(sprintf('Current time t = #f',times[n]))
        #subplot(4,1,2)
        plot(x,u)
        filnam = string("u\_t=",string(times[n+1]),".png")
        png(filnam)
        #axis([xL,xR,-1,1])
        #subplot(4,1,3)
        plot(times,masstot)
        filnam = string("mass\_t=",string(times[n+1]),".png")
        png(filnam)

        #axis([0,tend,40,50])
        #subplot(4,1,4)
        plot(times,momtot)
        filnam = string("momentum\_t=",string(times[n+1]),".png")
        png(filnam)
        #axis([0,tend,-50,50])
        show = true
        #sleep(5)
        #gui()
    end

end

function swe1d_nonlin()

    # Definition of the computational domain
    xL = 0.0                          # left border of spatial domain
    xR = 100.0                          # right border of spatial domain
    t0 = 0.0                          # initial time
    tend = 10.0                      # final time
    # Definition of the advective speed
    a = 1                           # advective speed
      # Discretization parameters
    IMAX = 100                      # number of velocity points in x-direction
    imax = IMAX-1                    # number of pressure points in x-direction
    dx = (xR-xL)/IMAX               # mesh size / mesh spacing
    dt0 = 0.1                       # time step
    #C = 0.9                         # Courant number
    # Definition of the gravity constant
    g = 9.81
    # Define the meshes
    #x  = range(xL, xR, length=IMAX)                # Distribute IMAX equidistant points between xL and xR for the velocity points
    x = Array{Float64,1}(undef, IMAX)
    #x = [xL:dx:xR]
    #xb = range(xL+dx/2, xR-dx/2, length=imax)      # Distribute IMAX equidistant points between xL and xR for the pressure points
    xb = Array{Float64,1}(undef, imax)
    eta = Array{Float64,1}(undef, imax)
    for i = 1:imax
        x[i]  =  xL + dx*(i-1)
        xb[i] =  xL + dx*(i-1/2)
        if x[i] < (xR-xL)/2.0
            eta[i] = 5.0
        else
            eta[i] = 1.0
        end
    end
    # Define the boundary conditions
    # For the velocities
    BCL = 0.0
    BCR = 0.0
    # Define the inital condition
    # First the free surface
    # Now initialize the velocities and the topography (at velocity points)
    u = Array{Float64,1}(undef, IMAX)
    h = Array{Float64,1}(undef, IMAX)
    etaj = Array{Float64,1}(undef, IMAX)
    hu = Array{Float64,1}(undef, IMAX)
    for i = 1:IMAX
        u[i] = 0                # Initial velocity
    #    h(i) = 0                # Bottom profile
        h[i] = 0.0
        if (i > 1 && i < IMAX)
            etaj[i] = 0.5*(eta[i-1]+eta[i])
            hu[i]   = max(0,h[i] + etaj[i])
        end
    end
    hu[1]    = max(0,h[1] + etaj[2])
    hu[IMAX] = max(0,h[IMAX] + etaj[imax])
    hi = Array{Float64,1}(undef, imax)
    for i = 1:imax
        hi[i]      = 0.5*(h[i]+h[i+1])
    end
    eta = max(eta,-hi)

    NMAX = 10000
    masstot = Array{Float64,1}(undef, NMAX)
    momtot  = Array{Float64,1}(undef, NMAX)
    masstot[1] = sum(eta + hi)
    momtot[1] = sum(hu.*u)
    #subplot(2,1,1)
    #plot(xb,eta,'o-')
    #subplot(2,1,2)
    #plot(x,u,'o-')
    # Initialize the current time with the initial time t0
    times = Array{Float64,1}(undef, NMAX)
    times[1] = t0
    i = 1
    # Define boundary conditions
    # Boundary condition on the left: q(xL,t) = q(xL,0)

    time0 = Dict{String,Float64}("tend"=>tend,"dt0"=>dt0)
    grid0 = Dict{String,Int64}("imax"=>imax,"IMAX"=>IMAX)
    grid  = Dict{String,Array{Float64,1}}("x"=>x,"xb"=>xb)
    state = Dict{String,Array{Float64,1}}("h"=>h,"u"=>u, "hi"=>hi, "hu"=>hu, "eta"=>eta, "etaj"=>etaj)
    bc    = Dict{String,Float64}("BCL"=>BCL,"BCR"=>BCR)
    physics = Dict{String,Float64}("grav"=>g)
    checking = Dict{String,Array{Float64,1}}("masstot"=>masstot,"momtot"=>momtot)

    # Loop in time
    timeloop(NMAX, times, time0, dx, grid0, grid, state, bc, physics, checking)
end

swe1d_nonlin()


text:
"""
$(SIGNATURES)

Apply geometry and time calibration to given hits.
"""
function calibrate(calibration::Calibration, hits)
    calibrated_hits = Vector{CalibratedHit}()
    has_trigger_mask = :trigger_mask ∈ fieldnames(eltype(hits))
    for hit in hits
        dom_id = hit.dom_id
        channel_id = hit.channel_id
        tot = hit.tot
        pos = calibration.pos[dom_id][channel_id+1]
        dir = calibration.dir[dom_id][channel_id+1]
        t0 = calibration.t0[dom_id][channel_id+1]
        t = hit.t + t0
        du = calibration.du[dom_id]
        floor = calibration.floor[dom_id]
        trigger_mask = has_trigger_mask ? hit.trigger_mask : 0
        c_hit = CalibratedHit(channel_id, dom_id, du, floor, t, tot,
                              pos, dir, t0, trigger_mask, Multiplicity(0,0))
        push!(calibrated_hits, c_hit)
    end
    calibrated_hits
end

"""
$(SIGNATURES)

Apply geometry and time calibration to given mc_hits.
"""
function calibrate(mc_hits::Vector{T},
                   calibration::Calibration,
                   event_info::Union{Nothing,MCEventInfo,DAQEventInfo} = nothing) where {T<:MCHit}
    calibrated_hits = Vector{CalibratedHit}()
    if event_info != nothing
        mctime = make_mc_time_converter(event_info)
    else
        mctime = x->x
    end
    for hit in mc_hits
        dom_id, channel_id = calibration.pmts[hit.pmt_id]
        tot = hit.a
        pos = calibration.pos[dom_id][channel_id+1]
        dir = calibration.dir[dom_id][channel_id+1]
        t0 = calibration.t0[dom_id][channel_id+1]
        t = mctime(hit.t)
        du = calibration.du[dom_id]
        floor = calibration.floor[dom_id]
        c_hit = CalibratedHit(channel_id, dom_id, du, floor, t, tot,
                              pos, dir, t0, hit.triggered, Multiplicity(0,0))
        push!(calibrated_hits, c_hit)
    end
    calibrated_hits
end

"""
    function floordist(calib::Calibration)

Calculates the average floor distance between neighboured modules.
"""
function floordist(calib::Calibration)
    floors = unique(sort([floor for floor ∈ values(calib.floor)]))
    floordistances = Float64[]
    for du ∈ values(calib.du)
        # collecting the z position for the same PMT channel on each module
        zs = [calib.pos[omkey2domid(calib, du, floor)][1][3] for floor ∈ floors]
        push!(floordistances, mean(diff(zs)))
    end
    mean(floordistances)
end

const SLEWS = [
    8.01, 7.52, 7.05, 6.59, 6.15, 5.74, 5.33, 4.95, 4.58, 4.22, 3.89, 3.56,
    3.25, 2.95, 2.66, 2.39, 2.12, 1.87, 1.63, 1.40, 1.19, 0.98, 0.78, 0.60,
    0.41, 0.24, 0.07, -0.10, -0.27, -0.43, -0.59, -0.75, -0.91, -1.08,
    -1.24, -1.41, -1.56, -1.71, -1.85, -1.98, -2.11, -2.23, -2.35, -2.47,
    -2.58, -2.69, -2.79, -2.89, -2.99, -3.09, -3.19, -3.28, -3.37, -3.46,
    -3.55, -3.64, -3.72, -3.80, -3.88, -3.96, -4.04, -4.12, -4.20, -4.27,
    -4.35, -4.42, -4.49, -4.56, -4.63, -4.70, -4.77, -4.84, -4.90, -4.97,
    -5.03, -5.10, -5.16, -5.22, -5.28, -5.34, -5.40, -5.46, -5.52, -5.58,
    -5.63, -5.69, -5.74, -5.80, -5.85, -5.91, -5.96, -6.01, -6.06, -6.11,
    -6.16, -6.21, -6.26, -6.31, -6.36, -6.41, -6.45, -6.50, -6.55, -6.59,
    -6.64, -6.68, -6.72, -6.77, -6.81, -6.85, -6.89, -6.93, -6.98, -7.02,
    -7.06, -7.09, -7.13, -7.17, -7.21, -7.25, -7.28, -7.32, -7.36, -7.39,
    -7.43, -7.46, -7.50, -7.53, -7.57, -7.60, -7.63, -7.66, -7.70, -7.73,
    -7.76, -7.79, -7.82, -7.85, -7.88, -7.91, -7.94, -7.97, -7.99, -8.02,
    -8.05, -8.07, -8.10, -8.13, -8.15, -8.18, -8.20, -8.23, -8.25, -8.28,
    -8.30, -8.32, -8.34, -8.37, -8.39, -8.41, -8.43, -8.45, -8.47, -8.49,
    -8.51, -8.53, -8.55, -8.57, -8.59, -8.61, -8.62, -8.64, -8.66, -8.67,
    -8.69, -8.70, -8.72, -8.74, -8.75, -8.76, -8.78, -8.79, -8.81, -8.82,
    -8.83, -8.84, -8.86, -8.87, -8.88, -8.89, -8.90, -8.92, -8.93, -8.94,
    -8.95, -8.96, -8.97, -8.98, -9.00, -9.01, -9.02, -9.04, -9.04, -9.04,
    -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04,
    -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04,
    -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04,
    -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04,
    -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04, -9.04,
    -9.04, -9.04
]

"""
$(SIGNATURES)

Return the time slewing for a ToT.
"""
slew(tot) = @inbounds ifelse(tot < 256, SLEWS[tot + 1], SLEWS[end])
slew(hit::AbstractHit) = slew(hit.tot)

omkey2domid(calib::Calibration, du, floor) = omkey2domid(calib, OMKey(du, floor))
omkey2domid(calib::Calibration, omkey::OMKey) = calib.omkeys[omkey]