text:
#Generic benchmark driver
for (testfunc, testname, longtestname, problem_sizes) in testdata
    for (n, t, size) in problem_sizes
        @timeit apply(testfunc, n, t) string(testname,"_",size) string(uppercase(size[1]),size[2:end]," ",longtestname," test")
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""
This file contains all the scripts used to produce the plots for the QuNet paper.
(Except for the Satellite plot which is in SatPlotDemo.jl)
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using DelimitedFiles

"""
Plot the performance statistics of greedy-multi-path vs the number of end-user pairs
"""
function plot_with_userpairs(max_pairs::Int64,
    num_trials::Int64)

    # The average routing costs between end-users sampled over num_trials for different numbers of end-users
    perf_data = []
    # The associated errors of the costs sampled over num_trials
    perf_err = []
    # Average numbers of paths used, sampled over num_trials for different numbers of end-users
    # e.g. [3,4,5]: 3 end-users found no path on average, 4 end-users found 1 path on average etc.
    path_data = []
    # Associated errors of path_data
    path_err = []

    grid_size = 10

    for i in 1:max_pairs
        println("Collecting for pairsize: $i")

        net = GridNetwork(grid_size, grid_size)

        # Collect performance statistics
        p, p_e, pat, pat_e = net_performance(net, num_trials, i, max_paths=4)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)

    end

    # Collect data for conditional probability of purification: (N2+N3)/∑N_i
    cpp = []
    cpp_err = []
    for i in 1:max_pairs
        tableau = path_data[i]
        errors = path_err[i]
        off = 1

        # Numerator and denominator
        num = tableau[2 + off] + tableau[3 + off] + tableau[4 + off]
        denom = tableau[1 + off] + tableau[2 + off] + tableau[3 + off] + tableau[4 + off]

        # Percentage errors (Ignore factor 100. Not needed)
        num_perr = (errors[2 + off] + errors[3 + off] + errors[4 + off])/num
        denom_perr = (errors[1 + off] + errors[2 + off] + errors[3 + off] + errors[4 + off])/denom

        data = num/denom
        data_err = (num_perr + denom_perr) * data
        push!(cpp, data)
        push!(cpp_err, data_err)
        i += 1
    end

    # Collect data for average number of paths used
    avepath = []
    for i in 1:max_pairs
        tableau = convert(Vector{Float64}, path_data[i])
        data = QuNet.ave_paths_used(tableau)

        # TODO: Include errors
        push!(avepath, data)
        i += 1
    end

    # Get values for x axis
    x = collect(1:max_pairs)

    # Extract data from performance
    loss = collect(map(x->x["loss"], perf_data))
    z = collect(map(x->x["Z"], perf_data))
    loss_err = collect(map(x->x["loss"], perf_err))
    z_err = collect(map(x->x["Z"], perf_err))

    # Extract data from path: PX is the rate of using X paths
    P0 = [path_data[i][1]/i for i in 1:max_pairs]
    P1 = [path_data[i][2]/i for i in 1:max_pairs]
    P2 = [path_data[i][3]/i for i in 1:max_pairs]
    P3 = [path_data[i][4]/i for i in 1:max_pairs]
    P4 = [path_data[i][5]/i for i in 1:max_pairs]

    P0e = [path_err[i][1]/i for i in 1:max_pairs]
    P1e = [path_err[i][2]/i for i in 1:max_pairs]
    P2e = [path_err[i][3]/i for i in 1:max_pairs]
    P3e = [path_err[i][4]/i for i in 1:max_pairs]
    P4e = [path_err[i][5]/i for i in 1:max_pairs]

    # Save data to txt
    open("data/userpairs.txt", "w") do io
        writedlm(io, ["max_pairs = $max_pairs",
        "num_trials = $num_trials",
        "grid_size = $grid_size",
        "perf_data", perf_data,
        "perf_err", perf_err,
        "path_data", path_data,
        "path_err", path_err,
        "cpp", cpp,
        "cpp_err", cpp_err,
        "avepath", avepath])
    end

    # Plot
    plot(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
    legend=:bottomright)
    plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$")
    xaxis!(L"$\textrm{Number of End User Pairs}$")
    savefig("plots/cost_userpair.png")
    savefig("plots/cost_userpair.pdf")

    plot(x, P0, ylims=(0,1), linewidth=2, yerr = P0e, label=L"$P_0$", legend= :right)
    plot!(x, P1, linewidth=2, yerr = P1e, label=L"$P_1$")
    plot!(x, P2, linewidth=2, yerr = P2e, label=L"$P_2$")
    plot!(x, P3, linewidth=2, yerr = P3e, label=L"$P_3$")
    plot!(x, P4, linewidth=2, yerr = P4e, label=L"$P_4$")
    plot!(x, cpp, linewidth=2, yerr = cpp_err, label=L"$P_{P}$")
    xaxis!(L"$\textrm{Number of End User Pairs}$")
    savefig("plots/path_userpair.png")
    savefig("plots/path_userpair.pdf")

    plot(x, avepath, linewidth=2, legend = false)
    xaxis!(L"$\textrm{Number of End User Pairs}$")
    yaxis!(L"$\textrm{Average Number of Paths Used}$")
    savefig("plots/avepath_userpair.png")
    savefig("plots/avepath_userpair.pdf")
end

"""
Plot the performance statistics of greedy-multi-path with respect to edge percolation
rate (The probability that a given edge is removed)
"""
function plot_with_percolations(perc_range::Tuple{Float64, Float64, Float64}, num_trials::Int64)
    # number of end-user pairs
    num_pairs = 1

    # Network to be percolated.
    grid_size = 10
    net = GridNetwork(grid_size, grid_size)

    perf_data = []
    perf_err = []
    path_data = []
    path_err = []

    for p in perc_range[1]:perc_range[2]:perc_range[3]
        println("Collecting for percolation rate: $p")

        # Percolate the network
        # perc_net = QuNet.percolate_edges(net, p)
        # refresh_graph!(perc_net)

        # Collect performance data with error, percolating the network edges
        p, p_e, pat, pat_e = net_performance(net, num_trials, num_pairs, max_paths=4,
        edge_perc_rate = p)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
    end

    # Get values for x axis
    x = collect(perc_range[1]:perc_range[2]:perc_range[3])

    # Extract performance data
    loss = collect(map(x->x["loss"], perf_data))
    z = collect(map(x->x["Z"], perf_data))
    loss_err = collect(map(x->x["loss"], perf_err))
    z_err = collect(map(x->x["Z"], perf_err))

    # Save data to txt
    open("data/percolation.txt", "w") do io
        writedlm(io, ["perc_range = $perc_range",
        "num_trials = $num_trials",
        "num_pairs = $num_pairs",
        "grid_size = $grid_size",
        "perf_data:", perf_data,
        "perf_err:", perf_err,
        "path_data:", path_data,
        "path_err:", path_err])
    end

    # Extract data from path: PX is the rate of using X paths
    P0 = [path_data[i][1]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P1 = [path_data[i][2]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P2 = [path_data[i][3]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P3 = [path_data[i][4]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P4 = [path_data[i][5]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]

    # Extract errors from path:
    P0e = [path_err[i][1]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P1e = [path_err[i][2]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P2e = [path_err[i][3]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P3e = [path_err[i][4]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
    P4e = [path_err[i][5]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]

    # Plot
    plot(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
    legend=:bottomright)
    plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$")
    xaxis!(L"$\textrm{Probability of Edge Removal}$")
    savefig("plots/cost_percolation.pdf")
    savefig("plots/cost_percolation.png")

    plot(x, P0, ylims=(0,1), linewidth=2, yerr = P0e, label=L"$P_0$", legend= :right)
    plot!(x, P1, linewidth=2, yerr = P1e, label=L"$P_1$")
    plot!(x, P2, linewidth=2, yerr = P2e, label=L"$P_2$")
    plot!(x, P3, linewidth=2, yerr = P3e, label=L"$P_3$")
    plot!(x, P4, linewidth=2, yerr = P3e, label=L"$P_4$")
    xaxis!(L"$\textrm{Probability of Edge Removal}$")
    savefig("plots/path_percolation.pdf")
    savefig("plots/path_percolation.png")
end

"""
Plot the performance statistics of greedy-multi-path with respect to the timedepth
of the graph
"""
function plot_with_timedepth(num_trials::Int64, max_depth::Int64)

    """
    For a given grid size, this function runs the greedy_multi_path routing algorithm
    on randomly placed end-user pairs.
    """
    function asymptotic_costs(gridsize::Int64)
        N = 10000
        G = GridNetwork(gridsize, gridsize)
        T = QuNet.TemporalGraph(G, 5, memory_costs = unit_costvector())
        p, dum1, dum2, dum3 = net_performance(T, N, 1, max_paths=4)
        return p["loss"], p["Z"]
    end

    # BEGIN MAIN

    num_pairs = 40
    grid_size = 10

    perf_data = []
    perf_err = []
    path_data = []
    path_err = []

    for i in 1:max_depth
        println("Collecting for time depth $i")
        G = GridNetwork(grid_size, grid_size)
        # Create a Temporal Graph from G with timedepth i
        T = QuNet.TemporalGraph(G, i, memory_costs = unit_costvector())
        # Get random pairs of asynchronus nodes
        user_pairs = make_user_pairs(T, num_pairs)
        # Get data
        p, p_e, pat, pat_e = net_performance(T, num_trials, num_pairs, max_paths=4)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
    end

    # Collect data for horizontal lines:
    println("Collecting data for asymptote")
    as = asymptotic_costs(grid_size)
    e_as = ones(length(1:max_depth)) * as[1]
    f_as = ones(length(1:max_depth)) * as[2]

    # Get values for x axis
    x = collect(1:max_depth)

    # Extract from performance data
    loss = collect(map(x->x["loss"], perf_data))
    z = collect(map(x->x["Z"], perf_data))
    loss_err = collect(map(x->x["loss"], perf_err))
    z_err = collect(map(x->x["Z"], perf_err))

    # Save data to txt
    open("data/temporal.txt", "w") do io
        writedlm(io, ["max_depth = $max_depth",
        "num_trials = $num_trials",
        "num_pairs = $num_pairs",
        "grid_size = $grid_size",
        "perf_data:", perf_data,
        "perf_err:", perf_err,
        "path_data:", path_data,
        "path_err:", path_err,
        "e_as", e_as,
        "f_as", f_as])
    end

    # Extract from path data
    P0 = [path_data[i][1]/num_pairs for i in 1:max_depth]
    P1 = [path_data[i][2]/num_pairs for i in 1:max_depth]
    P2 = [path_data[i][3]/num_pairs for i in 1:max_depth]
    P3 = [path_data[i][4]/num_pairs for i in 1:max_depth]
    P4 = [path_data[i][5]/num_pairs for i in 1:max_depth]

    P0e = [path_err[i][1]/num_pairs for i in 1:max_depth]
    P1e = [path_err[i][2]/num_pairs for i in 1:max_depth]
    P2e = [path_err[i][3]/num_pairs for i in 1:max_depth]
    P3e = [path_err[i][4]/num_pairs for i in 1:max_depth]
    P4e = [path_err[i][5]/num_pairs for i in 1:max_depth]

    # Plot
    # after seriestype: marker = (5)
    plot(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
    legend=:right)
    plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$")

    # Plot asymptote
    plot!(x, e_as, linestyle=:dot, color=:red, linewidth=2, label=L"$\textrm{Asymptotic } \eta$")
    plot!(x, f_as, linestyle=:dot, color=:green, linewidth=2, label=L"$\textrm{Asymptotic } F$")
    xaxis!(L"$\textrm{Time Depth of Tempral Meta-Graph}$")

    savefig("plots/cost_temporal.png")
    savefig("plots/cost_temporal.pdf")

    plot(x, P0, ylims=(0,1), linewidth=2, yerr = P0e, label=L"$P_0$", legend= :right)
    plot!(x, P1, linewidth=2, yerr = P1e, label=L"$P_1$")
    plot!(x, P2, linewidth=2, yerr = P2e, label=L"$P_2$")
    plot!(x, P3, linewidth=2, yerr = P3e, label=L"$P_3$")
    plot!(x, P4, linewidth=2, yerr = P4e, label=L"$P_4$")
    xaxis!(L"$\textrm{Time Depth of Temporal Meta-Graph}$")

    savefig("plots/path_temporal.png")
    savefig("plots/path_temporal.pdf")
end

"""
Plot performance statistics of greedy-multi-path with respect to grid size of the network
"""
function plot_with_gridsize(num_trials::Int64, num_pairs::Int64, min_size::Int64, max_size::Int64)
    @assert min_size < max_size
    @assert num_pairs*2 <= min_size^2 "Graph size too small for num_pairs"

    perf_data = []
    perf_err = []
    path_data = []
    path_err = []

    size_list = collect(min_size:1:max_size)
    for i in size_list
        println("Collecting for gridsize: $i")
        # Generate ixi graph:
        net = GridNetwork(i, i)

        # Collect performance statistics
        p, p_e, pat, pat_e = net_performance(net, num_trials, num_pairs, max_paths=4)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
    end

    # Get values for x axis
    x = collect(min_size:1:max_size)

    # Extract data from performance
    loss = collect(map(x->x["loss"], perf_data))
    z = collect(map(x->x["Z"], perf_data))
    loss_err = collect(map(x->x["loss"], perf_err))
    z_err = collect(map(x->x["Z"], perf_err))

    # Save data to txt
    open("data/gridsize.txt", "w") do io
        writedlm(io, ["min_size = $min_size",
        "max_size = $max_size",
        "num_trials = $num_trials",
        "num_pairs = $num_pairs",
        "perf_data:", perf_data,
        "perf_err:", perf_err,
        "path_data:", path_data,
        "path_err:", path_err])
    end

    # Extract from path data
    P0 = [path_data[i][1]/num_pairs for i in 1:(max_size-min_size)+1]
    P1 = [path_data[i][2]/num_pairs for i in 1:(max_size-min_size)+1]
    P2 = [path_data[i][3]/num_pairs for i in 1:(max_size-min_size)+1]
    P3 = [path_data[i][4]/num_pairs for i in 1:(max_size-min_size)+1]
    P4 = [path_data[i][5]/num_pairs for i in 1:(max_size-min_size)+1]

    P0e = [path_err[i][1]/num_pairs for i in 1:(max_size-min_size)+1]
    P1e = [path_err[i][2]/num_pairs for i in 1:(max_size-min_size)+1]
    P2e = [path_err[i][3]/num_pairs for i in 1:(max_size-min_size)+1]
    P3e = [path_err[i][4]/num_pairs for i in 1:(max_size-min_size)+1]
    P4e = [path_err[i][5]/num_pairs for i in 1:(max_size-min_size)+1]

    # Plot
    plot(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
    legend=:bottomright)
    plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$")
    xaxis!(L"$\textrm{Grid Size}$")
    savefig("plots/cost_gridsize.png")
    savefig("plots/cost_gridsize.pdf")

    plot(x, P0, ylims=(0,1), linewidth=2, yerr = P0e, label=L"$P_0$", legend= :right)
    plot!(x, P1, linewidth=2, yerr = P1e, label=L"$P_1$")
    plot!(x, P2, linewidth=2, yerr = P2e, label=L"$P_2$")
    plot!(x, P3, linewidth=2, yerr = P3e, label=L"$P_3$")
    plot!(x, P4, linewidth=2, yerr = P4e, label=L"$P_4$")
    xaxis!(L"$\textrm{Grid Size}$")
    savefig("plots/path_gridsize.png")
    savefig("plots/path_gridsize.pdf")
end

"""
Plot the performance statistics of greedy_multi_path for 1 random end-user pair
over a range of graph sizes and for different numbers of max_path (1,2,3) (ie. the
maximum number of paths that can be purified by a given end-user pair.)
"""
function plot_maxpaths_with_gridsize(num_trials::Int64, min_size::Int64, max_size::Int64)
    @assert min_size < max_size

    perf_data1 = []
    perf_data2 = []
    perf_data3 = []
    perf_data4 = []

    size_list = collect(min_size:1:max_size)
    for (j, data_array) in enumerate([perf_data1, perf_data2, perf_data3, perf_data4])
        println("Collecting for max_paths: $j")
        for i in size_list
            println("Collecting for gridsize: $i")
            # Generate ixi graph:
            net = GridNetwork(i, i)

            # Collect performance statistics
            performance, dummy, dummy, dummy = net_performance(net, num_trials, 1, max_paths=j)
            push!(data_array, performance)
        end
    end

    # Get values for x axis
    x = collect(min_size:1:max_size)

    # Extract data from performance data
    loss_arr1 = collect(map(x->x["loss"], perf_data1))
    z_arr1 = collect(map(x->x["Z"], perf_data1))
    loss_arr2 = collect(map(x->x["loss"], perf_data2))
    z_arr2 = collect(map(x->x["Z"], perf_data2))
    loss_arr3 = collect(map(x->x["loss"], perf_data3))
    z_arr3 = collect(map(x->x["Z"], perf_data3))
    loss_arr4 = collect(map(x->x["loss"], perf_data4))
    z_arr4 = collect(map(x->x["Z"], perf_data4))

    # Save data to txt
    open("data/maxpaths.txt", "w") do io
        writedlm(io, ["min_size = $min_size",
        "max_size = $max_size",
        "num_trials = $num_trials",
        "num_pairs = 1",
        "perf_data1:", perf_data1,
        "perf_data2:", perf_data2,
        "perf_data3", perf_data3,
        "perf_data4", perf_data4])
    end

    # Plot
    plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), color = :red, legend=:left)
    plot!(x, z_arr1, linewidth=2, label=L"$F_1$", linestyle=:dash, color =:red)
    plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", color =:blue)
    plot!(x, z_arr2, linewidth=2, label=L"$F_2$", linestyle=:dash, color =:blue)
    plot!(x, loss_arr3, linewidth=2, label=L"$\eta_3$", color =:green)
    plot!(x, z_arr3, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:green)
    plot!(x, loss_arr4, linewidth=2, label=L"$\eta_3$", color =:purple)
    plot!(x, z_arr4, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:purple)
    xaxis!(L"$\textrm{Grid Size}$")

    savefig("plots/cost_maxpaths.png")
    savefig("plots/cost_maxpaths.pdf")
end

"""
Draw a network with timedepth 1 and the greedy-paths chosen between 3 end user pairs.
"""
function draw_network_routing()
    timedepth = 3
    grid_size = 10

    # Index offset for asynchronus nodes
    off = grid_size^2 * timedepth
    # Choose asynchronus endusers
    userpairs = [(1 + off, 100 + off), (50 + off, 81 + off), (87 + off, 22 + off)]

    net = GridNetwork(grid_size, grid_size)
    T = QuNet.TemporalGraph(net, timedepth, memory_prob=1.0)
    QuNet.add_async_nodes!(T, userpairs)
    T_copy = deepcopy(T)
    user_paths, dum1, dum2 = QuNet.greedy_multi_path!(T_copy, QuNet.purify, userpairs)
    QuNet.plot_network(T.graph["Z"], user_paths, T.locs_x, T.locs_y)
end

"""
Two different temporal plots, one with memory and one without. While varying the
number of end-users, we compare the ratio of the depths of the graphs
"""
function plot_bandwidth_ratio_with_userpairs(num_trials::Int64, max_pairs::Int64)

    grid_size = 10
    time_depth = 50
    asynchronus_weight = 100

    # Generate ixi graph and extend it in time
    G = GridNetwork(grid_size, grid_size)
    # Extend in time with memory links:
    T_mem = QuNet.TemporalGraph(G, time_depth, memory_costs = unit_costvector())
    # Extend in time without memory
    T = QuNet.TemporalGraph(G, time_depth, memory_prob=0.0)

    plot_data = []
    error_data = []
    for i in 1:max_pairs
        println("Collecting for pairs : $i")
        raw_data = []
        for j in 1:num_trials
            # Get i random userpairs. Ensure src nodes are fixed on T=1, dst nodes are asynchronus.
            mem_user_pairs = make_user_pairs(T, i, src_layer=1, dst_layer=-1)
            user_pairs = make_user_pairs(T, i, src_layer=-1, dst_layer=-1)

            # Make copies of the network
            T_mem_copy = deepcopy(T_mem)
            T_copy = deepcopy(T)

            # Add async nodes
            QuNet.add_async_nodes!(T_mem_copy, mem_user_pairs, ϵ=asynchronus_weight)
            QuNet.add_async_nodes!(T_copy, user_pairs, ϵ=asynchronus_weight)

            # Get pathset data
            pathset_mem, dum1, dum2 = QuNet.greedy_multi_path!(T_mem_copy, QuNet.purify, mem_user_pairs, 4)
            pathset, dum1, dum2 = QuNet.greedy_multi_path!(T_copy, QuNet.purify, user_pairs, 4)
            # Pathset is an array of vectors containing edges describing paths between end-user pairs
            # Objective: find the largest timedepth used in the pathsets

            max_depth_mem = QuNet.max_timedepth(pathset_mem, T)
            max_depth = QuNet.max_timedepth(pathset, T)

            # Get the ratio of these two quantities. Add it to data array
            push!(raw_data, max_depth / max_depth_mem )
        end
        # Average the raw data, add it to plot data:
        push!(plot_data, mean(raw_data))
        # Get standard error
        push!(error_data, std(raw_data)/sqrt(num_trials - 1))
    end

    # Save data to txt
    open("data/bandwidth_with_userpairs.txt", "w") do io
        writedlm(io, ["grid_size = $grid_size",
        "num_trials = $num_trials",
        "time_depth = $time_depth",
        "max_pairs = $max_pairs",
        "asynchronus_weight = $asynchronus_weight",
        "bandwidth_ratios = $plot_data",
        "bandwidth_error = $error_data"])
    end

    # Plot
    x = collect(1:max_pairs)
    plot(x, plot_data, yerr = error_data, legend = false)
    xaxis!(L"$\textrm{Number of End User Pairs}$")

    savefig("plots/bandwidth_with_userpairs.png")
    savefig("plots/bandwidth_with_userpairs.pdf")
end

# For this one, let's just consider bandwidth. No ratio needed!
function plot_bandwidth_ratio_with_memory_rate(num_trials::Int64, perc_range::Tuple{Float64, Float64, Float64})

    grid_size = 10
    time_depth = 8
    num_pairs = 50
    asynchronus_weight = 100

    # Generate ixi graph and extend it in time
    G = GridNetwork(grid_size, grid_size)

    # Extend graph in time without memory
    T = QuNet.TemporalGraph(G, time_depth, memory_prob=0.0)

    plot_data = []
    error_data = []
    for i in perc_range[1]:perc_range[2]:perc_range[3]
        println("Collecting for memory percolation rate : $i")
        raw_data = []
        for j in 1:num_trials

            # Make a network with some probability of memory
            T_mem = QuNet.TemporalGraph(G, time_depth, memory_prob=i, memory_costs = unit_costvector())
            # Make a copy of the network without memory
            T_copy = deepcopy(T)

            # Get i random userpairs with asynchronus src and dst nodes.
            mem_user_pairs = make_user_pairs(T_mem, num_pairs, src_layer=-1, dst_layer=-1)
            user_pairs = make_user_pairs(T, num_pairs, src_layer=-1, dst_layer=-1)

            # Add async nodes
            QuNet.add_async_nodes!(T_mem, mem_user_pairs, ϵ=asynchronus_weight)
            QuNet.add_async_nodes!(T_copy, user_pairs, ϵ=asynchronus_weight)

            # Get pathset data
            pathset_mem, dum1, dum2 = QuNet.greedy_multi_path!(T_mem, QuNet.purify, mem_user_pairs, 1)
            pathset, dum1, dum2 = QuNet.greedy_multi_path!(T_copy, QuNet.purify, user_pairs, 1)
            # Pathset is an array of vectors containing edges describing paths between end-user pairs
            # Objective: find the largest timedepth used in the pathsets

            max_depth_mem = QuNet.max_timedepth(pathset_mem, T_mem)
            max_depth = QuNet.max_timedepth(pathset, T_copy)

            # Get the bandwidth of this quantity
            push!(raw_data, max_depth / max_depth_mem)
        end

        # Average the raw data, add it to plot data:
        push!(plot_data, mean(raw_data))
        # Get standard error
        push!(error_data, std(raw_data)/sqrt(num_trials - 1))
    end

    open("data/bandwidth_with_memory_rate.txt", "w") do io
        writedlm(io, ["grid_size = $grid_size",
        "num_trials = $num_trials",
        "time_depth = $time_depth",
        "num_pairs = $num_pairs",
        "perc_range = $perc_range",
        "asynchronus_weight = $asynchronus_weight",
        "bandwidth_ratios = $plot_data",
        "bandwidth_error = $error_data"])
    end

    # Plot
    x = collect(perc_range[1]:perc_range[2]:perc_range[3])
    plot(x, plot_data, yerr = error_data, legend = false)
    xaxis!(L"$\textrm{Proportion of Nodes with Quantum Memory}$")
    yaxis!(L"$\textrm{R}$")

    savefig("plots/bandwidth_with_memory_rate.png")
    savefig("plots/bandwidth_with_memory_rate.pdf")
end

# MAIN
"""
Uncomment functions to reproduce plots from the paper / create your own
"""
# Usage : (max_pairs::Int64, num_trials::Int64)
# plot_with_userpairs(50, 5000)

# Usage : (perc_range::Tuple{Float64, Float64, Float64}, num_trials::Int64)
# plot_with_percolations((0.0, 0.01, 0.7), 5000)

# Usage : (num_trials::Int64, max_depth::Int64)
# plot_with_timedepth(1000, 15)

# Usage : (num_trials::Int64, num_pairs::Int64, min_size::Int64, max_size::Int64)
# plot_with_gridsize(100, 40, 10, 150)

# Usage : (num_trials::Int64, num_pairs::Int64, min_size::Int64, max_size::Int64)
# plot_maxpaths_with_gridsize(5000, 10, 50)

# Usage : (num_trials::Int64, max_pairs::Int64)
# plot_bandwidth_ratio_with_userpairs(1000, 50)

# Usage : num_trials::Int64, perc_range::Tuple{Float64, Float64, Float64}
# plot_bandwidth_ratio_with_memory_rate(5000, (0.0, 0.05, 1.0))

# Usage : None
# draw_network_routing()

# Usage: (num_trials::Int64, max_depth::Int64)
# plot_nomultipath_with_timedepth(10, 10)

ext: jl

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
using Pkg; Pkg.activate("c:/git/SortingLab")

using Weave

weave("c:/git/SortingLab/README.jmd", out_path="c:/git/SortingLab/", doctype="github")


text:
__precompile__(true)

module InvariantSets

using Requires

using LazySets
import Distributions
import LazySets.Interval

# ==== SET TYPES ====
# Exported set types from LazySets.jl
export LazySet
# Balls
export Ball1, Ball2, BallInf, Ballp
# Polytopic Sets
export HPolygon, HPolytope, HPolyhedron, VPolygon, VPolytope
# others
export  Singleton,
        Interval,
        LineSegment,
        SingleEntryVector,
        Hyperplane,
        HalfSpace,
        Hyperrectangle,
        Ellipsoid,
        Zonotope,
        Universe

# ==== SET OPERATIONS ====
#  Exported lazy set operators from LazySets.jl
export +, ⊕, *, ∩
# Exported functionalities from LazySets.jl
export  dim,
        isbounded,
        tosimplehrep,
        constraint_list,
        ρ, support_function,
        σ, support_vector,
        is_interior_point,
        ∈, ≈, ⊆, ==,
        isequivalent,
        rand, an_element,
        chebyshev_center, intersection,
        sample # needs Distributions to be loaded
        # TODO: distance_point_to_set
# Convenience operator for set equality check (⫓), because == checks for type and field equality
⫓(X::LazySet, Y::LazySet) = X ⊆ Y && Y ⊆ X
# TODO: Check if ⫓(X::LazySet, Y::LazySet) = isequivalent(X,Y) works for all set types
export ⫓

# TODO: add docstring to +ᶜ and -ᶜ
export +ᶜ, ⊕ᶜ, -ᶜ, ⊖ᶜ, *ᶜ, ∩ᶜ, reflect, translate
include(raw"concrete_operator.jl")

export  ρ_exact, σ_exact, ρ_matrix,
        canonical_length, marginal_enlargment
include(raw"helper\lazy_sets.jl")
include(raw"helper\control_systems.jl")

abstract type AbstractAlgorithm end
export preset,
       state_constraints,
       maximum_invariant_set,
       maximum_control_invariant_set,
       maximal_RPI_set
# TODO: add maximal_RPI_set
include(raw"maximal_invariance/main.jl")
export minimal_RPI_set,
       implicit_rakovic
# TODO: add implicit_rakovic once paper is published
# TODO: add tests for F_lazy (with the straightforward definition)
# TODO: Document  Schulze Algorithm
# TODO: add minimal_RPI_set
include(raw"minimal_invariance/main.jl")

# Derived sets used in MPC and control systems
export terminal_set,
       feasible_set,
       tightened_state_constraint,
       tightened_input_constraint,
       tightened_constraints
# TODO: documentation
# TODO: add `tightened_constraints`once paper is published
include(raw"applied_sets/control_systems.jl")
export refine
include(raw"approximation.jl")

# ===================================================
# Load external packages on-demand (using 'Requires')
# ===================================================

#TODO: fix add_constraint! for x and X (array and vecotr) add write docu new!
include(raw"init.jl")

end # module

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
module GridapRunTests

using Test

@time @testset "Helpers" begin include("HelpersTests/runtests.jl") end

@time @testset "Io" begin include("IoTests/runtests.jl") end

@time @testset "Algebra" begin include("AlgebraTests/runtests.jl") end

@time @testset "Arrays" begin include("ArraysTests/runtests.jl") end

@time @testset "TensorValues" begin include("TensorValuesTests/runtests.jl") end

@time @testset "Fields" begin include("FieldsTests/runtests.jl") end

@time @testset "Polynomials" begin include("PolynomialsTests/runtests.jl") end

@time @testset "ReferenceFEs" begin include("ReferenceFEsTests/runtests.jl") end

@time @testset "Geometry" begin include("GeometryTests/runtests.jl") end

@time @testset "CellData" begin include("CellDataTests/runtests.jl") end

@time @testset "Visualization" begin include("VisualizationTests/runtests.jl") end

@time @testset "FESpaces (1/2)" begin include("FESpacesTests/runtests_1.jl") end

@time @testset "FESpaces (2/2)" begin include("FESpacesTests/runtests_2.jl") end

@time @testset "MultiField" begin include("MultiFieldTests/runtests.jl") end

end # module


text:
module ErrorTypes

struct Unsafe end
const unsafe = Unsafe()

"""
    Ok{T}

The success state of a `Result{T, E}`, carrying an object of type `T`.
For convenience, `Ok(x)` creates a dummy value that can be converted to the
appropriate `Result type`.
"""
struct Ok{T}
    x::T
    Ok{T}(::Unsafe, x::T) where T = new{T}(x)
end

"""
    Err

The error state of a `Result{O, E}`, carrying an object of type `E`.
For convenience, `Err(x)` creates a dummy value that can be converted to the
appropriate `Result type`.
"""
struct Err{E}
    x::E
    Err{E}(::Unsafe, x::E) where E = new{E}(x)
end

struct ResultConstructor{T, K <: Union{Ok, Err}}
    x::T
end

Ok{T}(x::T2) where {T, T2 <: T} = ResultConstructor{T, Ok}(convert(T, x))
Err{T}(x::T2) where {T, T2 <: T} = ResultConstructor{T, Err}(convert(T, x))
Ok(x::T) where T = Ok{T}(x)
Err(x::T) where T = Err{T}(x)

"""
    Result{O, E}

A sum type of either `Ok{O}` or `Err{E}`. Used as return value of functions that
can error with an informative error object of type `E`.
"""
struct Result{O, E}
    x::Union{Ok{O}, Err{E}}
end

function Result{T1, E1}(x::Result{T2, E2}) where {T1, T2 <: T1, E1, E2 <: E1}
    inner = x.x
    return Result{T1, E1}(inner isa Err ? Err(inner.x) : Ok(inner.x))
end

function Base.convert(::Type{Result{T1, E1}}, x::Result{T2, E2}) where {T1, T2 <: T1, E1, E2 <: E1}
    Result{T1, E1}(x)
end

# Bizarrely, we need this method to avoid recursion when creating a Result
Base.convert(::Type{T}, x::T) where {T <: Result} = x

# Strict supertype accepted
Result{O, E}(x::ResultConstructor{O2, Ok}) where {O, E, O2 <: O} = Result{O, E}(Ok{O}(unsafe, x.x))

# We have <:Err here to also cover ResultConstructor{E2, Union{}} propagated with @?
Result{O, E}(x::ResultConstructor{E2, <:Err}) where {O, E, E2 <: E} = Result{O, E}(Err{E}(unsafe, x.x))

# ResultConstructor propagated with @? can convert to any Option
Result{O, Nothing}(x::ResultConstructor{E, Union{}}) where {O, E} = Result{O, Nothing}(Err{Nothing}(unsafe, nothing))
Base.convert(::Type{T}, x::ResultConstructor) where {T <: Result} = T(x)

function Base.show(io::IO, x::Result)
    print(io, typeof(x), '(', Base.typename(typeof(x.x)).name, '(', repr(x.x.x), "))")
end

"""
    Option{T}

Alias for `Result{T, Nothing}`
"""
const Option{T} = Result{T, Nothing}

"""
	is_error(x::Result)

Check if `x` contains an error value.

# Example

```jldoctest
julia> is_error(none(Int)), is_error(some(5))
(true, false)
```
"""
is_error(x::Result) = x.x isa Err

Option(x::Result{T, E}) where {T, E} = is_error(x) ? none(T) : Option{T}(x.x)
function Result(x::Option{T}, v::E) where {T, E}
    data = x.x
    return Result{T, E}(data isa Err ? Err(v) : x.x)
end

some(x::T) where T = Option{T}(Ok{T}(unsafe, x))

const none = ResultConstructor{Nothing, Err}(nothing)
none(x::Type{T}) where T = Option{T}(Err{Nothing}(unsafe, nothing))

function Base.show(io::IO, x::Option{T}) where T
    inner = x.x
    if inner isa Err
        print(io, "none(", T, ')')
    else
        # If the inner value is exactly T, elide printing it for brevity.
        value = inner.x
        if typeof(value) === T
            print(io, "some(", repr(x.x.x), ')')
        else
            print(io, "Option{", T, "}(some(", repr(x.x.x), "))")
        end
    end
end

"""
    @?(expr)

Propagate a `Result` with `Err` value to the outer function.
Evaluate `expr`, which should return a `Result`. If it contains an `Ok` value `x`,
evaluate to the unwrapped value `x`. Else, evaluates to `return Err(x)`.

# Example

```jldoctest
julia> (f(x::Option{T})::Option{T}) where T = Ok(@?(x) + one(T));

julia> f(some(1.0)), f(none(Int))
(some(2.0), none(Int64))
```
"""
macro var"?"(expr)
    quote
        local res = $(esc(expr))
        isa(res, Result) || throw(TypeError(Symbol("@?"), Result, res))
        local data = res.x
        data isa Ok ? data.x : return ResultConstructor{typeof(data.x), Union{}}(data.x)
    end
end

macro unwrap_t_or(expr, T, exec)
    quote
        local res = $(esc(expr))
        isa(res, Result) || throw(TypeError(Symbol("@unwrap_or"), Result, res))
        local data = res.x
        data isa $(esc(T)) ? data.x : $(esc(exec))
    end
end

# Note: jldoctests crashes on this ATM, even if it works correctly.
"""
    @unwrap_or(expr, exec)

Evaluate `expr` to a `Result`. If `expr` is a error value, evaluate
`exec` and return that. Else, return the wrapped value in `expr`.

# Examples

```
julia> safe_inv(x)::Option{Float64} = iszero(x) ? none : Ok(1/x);

julia> function skip_inv_sum(it)
    sum = 0.0
    for i in it
        sum += @unwrap_or safe_inv(i) continue
    end
    sum
end;

julia> skip_inv_sum([2,1,0,1,2])
3.0
```
"""
macro unwrap_or(expr, exec)
    return :(@unwrap_t_or $(esc(expr)) Ok $(esc(exec)))
end

"""
    @unwrap_error_or(expr, exec)

Same as `@unwrap_or`, but unwraps errors.
"""
macro unwrap_error_or(expr, exec)
    return :(@unwrap_t_or $(esc(expr)) Err $(esc(exec)))
end

@noinline throw_unwrap_err() = error("unwrap on unexpected type")

"""
    unwrap(x::Result)

If `x` is of the associated error type, throw an error. Else, return the contained
result type.
"""
unwrap(x::Result) = @unwrap_or x throw_unwrap_err()

"""
    unwrap_error(x::Result)

If `x` contains an `Err`, return the content of the `Err`. Else, throw an error.
"""
unwrap_error(x::Result) = @unwrap_error_or x throw_unwrap_err()

"""
    expect(x::Result, s::AbstractString)

If `x` is of the associated error type, error with message `s`. Else, return
the contained result type.
"""
expect(x::Result, s::AbstractString) = @unwrap_or x error(s)

"""
    expect_error(x::Result, s::AbstractString)

If `x` contains an `Err`, return the content of the `Err`. Else, throw an error
with message `s`.
"""
expect_error(x::Result, s::AbstractString) = @unwrap_error_or x error(s)

"""
    and_then(f, ::Type{T}, x::Result{O, E})

If `is` a result value, apply `f` to `unwrap(x)`, else return the error value. Always returns a `Result{T, E}`.

__WARNING__
If `f(unwrap(x))` is not a `T`, this functions throws an error.
"""
function and_then(f, ::Type{T}, x::Result{O, E})::Result{T, E} where {T, O, E}
    data = x.x
    Result{T, E}(data isa Ok ? Ok(f(data.x)) : Err(data.x))
end

"""
    map_or(f, x::Result, v)

If `x` is a result value, return `f(unwrap(x))`. Else, return `v`.
"""
map_or(f, x::Result{O, E}, v) where {O, E} = f(@unwrap_or x return v)

"""
    unwrap_or(x::Result, v)

If `x` is an error value, return `v`. Else, unwrap `x` and return its content.
"""
unwrap_or(x::Result, v) = @unwrap_or x v

"""
    unwrap_error_or(x::Result, v)

Like `unwrap_or`, but unwraps an error.
"""
unwrap_error_or(x::Result, v) = @unwrap_error_or x v

"""
    flatten(x::Option{Option{T}})

Convert an `Option{Option{T}}` to an `Option{T}`.

# Examples

```jldoctest
julia> flatten(some(some("x")))
some("x")

julia> flatten(some(none(Int)))
none(Int)
```
"""
flatten(x::Option{Option{T}}) where T = unwrap_or(x, none(T))

"""
    base(x::Option{T})

Convert an `Option{T}` to a `Union{Some{T}, Nothing}`.
"""
base(x::Option{T}) where T = Some(@unwrap_or x (return nothing))

export Ok, Err, Result, Option,
    is_error, some, none,
    unwrap, unwrap_error, expect, expect_error,
    and_then, unwrap_or, map_or, unwrap_error_or, flatten, base,
    @?, @unwrap_or, @unwrap_error_or

end


text:
"""
Vector of a judiVector and a judiWeight
"""
mutable struct judiVStack{vDT<:Number}
    m::Integer
    n::Integer
    components::Array{Any, 1}
end

function vcat(x::judiVector, y::judiWeights)
    components = Array{Any}(undef, 2)
    components[1] = x
    components[2] = y
    m = length(x)+length(y)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiWeights, y::judiVector)
    components = Array{Any}(undef, 2)
    components[1] = x
    components[2] = y
    m = length(x)+length(y)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiVStack, y::Union{judiWeights, judiVector})
    components = Array{Any}(undef, length(x.components) + 1)
    for i=1:length(x.components)
        components[i] = x.components[i]
    end
    components[end] = y
    m = x.m+length(y)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::Union{judiWeights, judiVector}, y::judiVStack)
    components = Array{Any}(undef, length(y.components) + 1)
    components[1] = x
    for i=2:length(y.components)+1
        components[i] = y.components[i-1]
    end
    m = y.m + length(x)
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function vcat(x::judiVStack, y::judiVStack)
    nx = length(x.components)
    ny = length(y.components)
    components = Array{Any}(undef, nx+ny)
    for i=1:nx
        components[i] = x.components[i]
    end
    for i=nx+1:nx+ny
        components[i] = y.components[i-nx]
    end
    m = x.m + y.m
    n = 1
    return judiVStack{Float32}(m, n, components)
end

function *(F::joAbstractLinearOperator, v::judiVStack)
    return sum(F.fop[i]*v[i] for i=1:length(v.components))
end

############################################################
## overloaded Base functions

# conj(jo)
conj(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.m,a.n,a.components)

# transpose(jo)
transpose(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.n,a.m,a.components)

# adjoint(jo)
adjoint(a::judiVStack{vDT}) where vDT =
    judiVStack{vDT}(a.n,a.m,a.components)

##########################################################
# Utilities

size(x::judiVStack) = (x.m, x.n)
size(x::judiVStack, ind::Integer) = (x.m, x.n)[ind]

length(x::judiVStack) = x.m

eltype(v::judiVStack{vDT}) where {vDT} = vDT

similar(x::judiVStack) = judiVStack{Float32}(x.m, x.n, 0f0 .* x.components)

similar(x::judiVStack,  element_type::DataType, dims::Union{AbstractUnitRange, Integer}...) = similar(x)

getindex(x::judiVStack, a) = x.components[a]

firstindex(x::judiVStack) = 1

lastindex(x::judiVStack) = length(x.components)

dot(x::judiVStack, y::judiVStack) = sum(dot(x[i],y[i]) for i=1:length(x.components))

function norm(x::judiVStack, order::Real=2)
    if order == Inf
        return max([norm(x[i], Inf) for i=1:length(x.components)]...)
    end
    sum(norm(x[i], order)^order for i=1:length(x.components))^(1/order)
end

iterate(S::judiVStack, state::Integer=1) = state > length(S.components) ? nothing : (S.components[state], state+1)

isfinite(S::judiVStack) = all(isfinite(c) for c in S)

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

# minus
-(a::judiVStack) = -1*a
+(a::judiVStack, b::judiVStack) = judiVStack{Float32}(a.m, a.n, a.components + b.components)
-(a::judiVStack, b::judiVStack) = judiVStack{Float32}(a.m, a.n, a.components - b.components)
+(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n,a.components .+ b)
-(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n,a.components .- b)
-(a::Number, b::judiVStack) = judiVStack{Float32}(b.m, b.n, a .- b.components)

*(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n, b .* a.components)
*(a::Number, b::judiVStack) = b * a
+(a::Number, b::judiVStack) = b + a

/(a::judiVStack, b::Number) = judiVStack{Float32}(a.m, a.n, a.components ./ b)

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

BroadcastStyle(::Type{judiVStack}) = Base.Broadcast.DefaultArrayStyle{1}()

broadcasted(::typeof(+), x::judiVStack, y::judiVStack) = x + y
broadcasted(::typeof(-), x::judiVStack, y::judiVStack) = x - y
broadcasted(::typeof(+), x::judiVStack, y::Number) = x + y
broadcasted(::typeof(-), x::judiVStack, y::Number) = x - y
broadcasted(::typeof(+), y::Number, x::judiVStack) = x + y
broadcasted(::typeof(-), y::Number, x::judiVStack) = x - y

function broadcasted(::typeof(*), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] .* y.components[j]
    end
    return z
end

function broadcasted!(::typeof(*), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] .* y.components[j]
    end
    return z
end

function broadcasted(::typeof(/), x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] = x.components[j] ./ y.components[j]
    end
    return z
end

function broadcasted(::typeof(*), x::judiVStack, y::Number)
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] .*= y
    end
    return z
end

broadcasted(::typeof(*), y::Number, x::judiVStack) = x .* y

function broadcasted(::typeof(/), x::judiVStack, y::Number)
    z = deepcopy(x)
    for j=1:length(x.components)
        z.components[j] ./= y
    end
    return z
end

function materialize!(x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    for j=1:length(x.components)
        try
            x.components[j].data .= y.components[j].data
        catch e
            x.components[j].weights .= y.components[j].weights
        end
    end
    return x
end

function broadcast!(identity, x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    copy!(x,y)
end

function broadcasted(identity, x::judiVStack)
    return x
end

function copy!(x::judiVStack, y::judiVStack)
    size(x) == size(y) || throw(judiWeightsException("dimension mismatch"))
    for j=1:length(x.components)
        try
            x.components[j].data .= y.components[j].data
        catch e
            x.components[j].weights .= y.components[j].weights
        end
    end
end

function isapprox(x::judiVStack, y::judiVStack; rtol::Real=sqrt(eps()), atol::Real=0)
    x.m == y.m || throw("Shape error")
    all(isapprox(xx, yy; rtol=rtol, atol=atol) for (xx, yy)=zip(x.components, y.components))
end

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

function A_mul_B!(x::judiWeights, F::joCoreBlock, y::judiVStack)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    x.weights = z.weights
end

function A_mul_B!(x::judiVStack, F::joCoreBlock, y::judiWeights)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.components)
        try
            x.components[j].data .= z.components[j].data
        catch e
            x.components[j].weights .= z.components[j].weights
        end
    end
end

mul!(x::judiWeights, J::joCoreBlock, y::judiVStack) = A_mul_B!(x, J, y)
mul!(x::judiVStack, J::joCoreBlock, y::judiWeights) = A_mul_B!(x, J, y)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
#=

This file contains RT elemental-related functions

=#

"Elemental single-scattering layer"
function elemental!(pol_type, SFI::Bool,
                            τ_sum::AbstractArray,#{FT2,1}, #Suniti
                            dτ_λ::AbstractArray{FT,1},  # dτ_λ: total optical depth of elemental layer (per λ)
                            dτ::FT,                     # dτ:   scattering optical depth of elemental layer (scalar)
                            ϖ_λ::AbstractArray{FT,1},   # ϖ: single scattering albedo of elemental layer (per λ, absorptions by gases included)
                            ϖ::FT,                      # ϖ: single scattering albedo of elemental layer (no trace gas included)
                            Z⁺⁺::AbstractArray{FT,2},   # Z matrix
                            Z⁻⁺::AbstractArray{FT,2},
                            m::Int,                     # m: fourier moment
                            ndoubl::Int,                # ndoubl: number of doubling computations needed
                            scatter::Bool,              # scatter: flag indicating scattering
                            quad_points::QuadPoints{FT2}, # struct with quadrature points, weights,
                            added_layer::AddedLayer{FT},
                            I_static,
                            architecture) where {FT<:Union{AbstractFloat, ForwardDiff.Dual},FT2}

    @unpack r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, J₀⁺, J₀⁻ = added_layer
    @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, iμ₀ = quad_points
    arr_type = array_type(architecture)
    # Need to check with paper nomenclature. This is basically eqs. 19-20 in vSmartMOM

    # Later on, we can have Zs also vary with index, pretty easy here:
    # Z⁺⁺_ = repeat(Z⁺⁺, 1, 1, 1)
    Z⁺⁺_ = reshape(Z⁺⁺, (size(Z⁺⁺,1), size(Z⁺⁺,2),1))
    # Z⁻⁺_ = repeat(Z⁻⁺, 1, 1, 1)
    Z⁻⁺_ = reshape(Z⁻⁺, (size(Z⁺⁺,1), size(Z⁺⁺,2),1))

    D = Diagonal(arr_type(repeat(pol_type.D, size(qp_μ,1))))
    I₀_NquadN = arr_type(zeros(FT,size(qp_μN,1))); #incident irradiation
    i_start   = pol_type.n*(iμ₀-1) + 1
    i_end     = pol_type.n*iμ₀
    I₀_NquadN[iμ₀Nstart:i_end] = pol_type.I₀

    device = devi(architecture)

    # If in scattering mode:
    if scatter

        NquadN = length(qp_μN)

        # Needs explanation still, different weigths:
        wct0  = m == 0 ? FT(0.50) * ϖ * dτ     : FT(0.25) * ϖ * dτ
        wct02 = m == 0 ? FT(0.50)              : FT(0.25)
        wct   = m == 0 ? FT(0.50) * ϖ * wt_μN  : FT(0.25) * ϖ * wt_μN
        wct2  = m == 0 ? wt_μN/2               : wt_μN/4

        # Get the diagonal matrices first
        d_qp  = Diagonal(1 ./ qp_μN)
        d_wct = Diagonal(wct)

        # Calculate r⁻⁺ and t⁺⁺

        # Version 1: no absorption in batch mode (initiation of a single scattering layer with no or low absorption)
        if false #maximum(dτ_λ) < 0.0001
            # R⁻⁺₀₁(λ) = M⁻¹(0.5ϖₑ(λ)Z⁻⁺C)δ (See Eqs.7 in Raman paper draft)
            r⁻⁺[:,:,:] .= d_qp * Z⁻⁺ * (d_wct * dτ)
            # T⁺⁺₀₁(λ) = {I-M⁻¹[I - 0.5*ϖₑ(λ)Z⁺⁺C]}δ (See Eqs.7 in Raman paper draft)
            t⁺⁺[:,:,:] .= I_static - (d_qp * ((I_static - Z⁺⁺ * d_wct) * dτ))
            if SFI
                # Reminder: Add equation here what it does
                expk = exp.(-τ_sum/qp_μ[iμ₀]) #exp(-τ(z)/μ₀)
                # J₀⁺ = 0.5[1+δ(m,0)]M⁻¹ϖₑ(λ)Z⁺⁺τI₀exp(-τ(z)/μ₀)
                J₀⁺[:,1,:] .= ((d_qp * Z⁺⁺ * I₀_NquadN) * wct0) .* expk'
                # J₀⁻ = 0.5[1+δ(m,0)]M⁻¹ϖₑ(λ)Z⁻⁺τI₀exp(-τ(z)/μ₀)
                J₀⁻[:,1,:] .= ((d_qp * Z⁻⁺ * I₀_NquadN) * wct0) .* expk'

            end
        else
            #Version 2: More computationally intensive definition of a single scattering layer with variable (0-∞) absorption
            # Version 2: with absorption in batch mode, low tau_scatt but higher tau_total, needs different equations
            kernel! = get_elem_rt!(device)
            event = kernel!(r⁻⁺, t⁺⁺, ϖ_λ, dτ_λ, Z⁻⁺, Z⁺⁺, qp_μN, wct, ndrange=size(r⁻⁺));
            wait(device, event)
            synchronize_if_gpu()

            if SFI
                kernel! = get_elem_rt_SFI!(device)
                event = kernel!(J₀⁺, J₀⁻, ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, qp_μN, ndoubl, wct02, pol_type.n, arr_type(pol_type.I₀), iμ₀, D, ndrange=size(J₀⁺))
                wait(device, event)
            end
            #ii = pol_type.n*(iμ0-1)+1
            #@show 'B',iμ0,  r⁻⁺[1,ii,1]/(J₀⁻[1,1,1]*wt_μ[iμ0]), r⁻⁺[1,ii,1], J₀⁻[1,1,1]*wt_μ[iμ0], J₀⁺[1,1,1]*wt_μ[iμ0]
            synchronize_if_gpu()
        end
        # Apply D Matrix
        apply_D_matrix_elemental!(ndoubl, pol_type.n, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻)

        if SFI
            apply_D_matrix_elemental_SFI!(ndoubl, pol_type.n, J₀⁻)
        end
    else
        # Note: τ is not defined here
        t⁺⁺[:] = Diagonal{exp(-τ ./ qp_μN)}
        t⁻⁻[:] = Diagonal{exp(-τ ./ qp_μN)}
    end
    #@pack! added_layer = r⁺⁻, r⁻⁺, t⁻⁻, t⁺⁺, J₀⁺, J₀⁻
end

@kernel function get_elem_rt!(r⁻⁺, t⁺⁺, ϖ_λ, dτ_λ, Z⁻⁺, Z⁺⁺, qp_μN, wct2)
    i, j, n = @index(Global, NTuple)

    if (wct2[j]>1.e-8)
        # 𝐑⁻⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁻⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ+μⱼ)) ̇(1 - exp{-τ ̇(1/μᵢ + 1/μⱼ)}) ̇𝑤ⱼ
        r⁻⁺[i,j,n] = ϖ_λ[n] * Z⁻⁺[i,j] * (qp_μN[j] / (qp_μN[i] + qp_μN[j])) * (1 - exp(-dτ_λ[n] * ((1 / qp_μN[i]) + (1 / qp_μN[j])))) * (wct2[j])

        if (qp_μN[i] == qp_μN[j])
            # 𝐓⁺⁺(μᵢ, μᵢ) = (exp{-τ/μᵢ} + ϖ ̇𝐙⁺⁺(μᵢ, μᵢ) ̇(τ/μᵢ) ̇exp{-τ/μᵢ}) ̇𝑤ᵢ
            if i == j
                t⁺⁺[i,j,n] = exp(-dτ_λ[n] / qp_μN[i])*(1 + ϖ_λ[n] * Z⁺⁺[i,i] * (dτ_λ[n] / qp_μN[i]) * wct2[i])
            else
                t⁺⁺[i,j,n] = 0.0
            end
        else

            # 𝐓⁺⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁺⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ-μⱼ)) ̇(exp{-τ/μᵢ} - exp{-τ/μⱼ}) ̇𝑤ⱼ
            # (𝑖 ≠ 𝑗)
            t⁺⁺[i,j,n] = ϖ_λ[n] * Z⁺⁺[i,j] * (qp_μN[j] / (qp_μN[i] - qp_μN[j])) * (exp(-dτ_λ[n] / qp_μN[i]) - exp(-dτ_λ[n] / qp_μN[j])) * wct2[j]
        end
    else
        r⁻⁺[i,j,n] = 0.0
        if i==j
            t⁺⁺[i,j,n] = exp(-dτ_λ[n] / qp_μN[i]) #Suniti
        else
            t⁺⁺[i,j,n] = 0.0
        end
    end

end

@kernel function get_elem_rt_SFI!(J₀⁺, J₀⁻, ϖ_λ, dτ_λ, τ_sum, Z⁻⁺, Z⁺⁺, qp_μN, ndoubl, wct02, nStokes ,I₀, iμ0, D)
    i_start  = nStokes*(iμ0-1) + 1
    i_end    = nStokes*iμ0

    i, _, n = @index(Global, NTuple) ##Suniti: What are Global and Ntuple?
    FT = eltype(I₀)
    J₀⁺[i, 1, n]=0
    J₀⁻[i, 1, n]=0

    Z⁺⁺_I₀ = FT(0.0);
    Z⁻⁺_I₀ = FT(0.0);
    for ii = i_start:i_end
        Z⁺⁺_I₀ += Z⁺⁺[i,ii] * I₀[ii-i_start+1]
        Z⁻⁺_I₀ += Z⁻⁺[i,ii] * I₀[ii-i_start+1]
    end

    if (i>=i_start) && (i<=i_end)
        ctr = i-i_start+1
        # J₀⁺ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁺⁺ * I₀ * (dτ(λ)/μ₀) * exp(-dτ(λ)/μ₀)
        J₀⁺[i, 1, n] = wct02 * ϖ_λ[n] * Z⁺⁺_I₀ * (dτ_λ[n] / qp_μN[i]) * exp(-dτ_λ[n] / qp_μN[i])
    else
        # J₀⁺ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁺⁺ * I₀ * [μ₀ / (μᵢ - μ₀)] * [exp(-dτ(λ)/μᵢ) - exp(-dτ(λ)/μ₀)]
        J₀⁺[i, 1, n] = wct02 * ϖ_λ[n] * Z⁺⁺_I₀ * (qp_μN[i_start] / (qp_μN[i] - qp_μN[i_start])) * (exp(-dτ_λ[n] / qp_μN[i]) - exp(-dτ_λ[n] / qp_μN[i_start]))
    end
    #J₀⁻ = 0.25*(1+δ(m,0)) * ϖ(λ) * Z⁻⁺ * I₀ * [μ₀ / (μᵢ + μ₀)] * [1 - exp{-dτ(λ)(1/μᵢ + 1/μ₀)}]
    J₀⁻[i, 1, n] = wct02 * ϖ_λ[n] * Z⁻⁺_I₀ * (qp_μN[i_start] / (qp_μN[i] + qp_μN[i_start])) * (1 - exp(-dτ_λ[n] * ((1 / qp_μN[i]) + (1 / qp_μN[i_start]))))

    J₀⁺[i, 1, n] *= exp(-τ_sum[n]/qp_μN[i_start])
    J₀⁻[i, 1, n] *= exp(-τ_sum[n]/qp_μN[i_start])

    if ndoubl >= 1
        J₀⁻[i, 1, n] = D[i,i]*J₀⁻[i, 1, n] #D = Diagonal{1,1,-1,-1,...Nquad times}
    end
end

@kernel function apply_D_elemental!(ndoubl, pol_n, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻)
    i, j, n = @index(Global, NTuple)

    if ndoubl < 1
        ii = mod(i, pol_n)
        jj = mod(j, pol_n)
        if ((ii <= 2) & (jj <= 2)) | ((ii > 2) & (jj > 2))
            r⁺⁻[i,j,n] = r⁻⁺[i,j,n]
            t⁻⁻[i,j,n] = t⁺⁺[i,j,n]
        else
            r⁺⁻[i,j,n] = -r⁻⁺[i,j,n]
            t⁻⁻[i,j,n] = -t⁺⁺[i,j,n]
        end
    else
        if mod(i, pol_n) > 2
            r⁻⁺[i,j,n] = - r⁻⁺[i,j,n]
        end
    end
end

@kernel function apply_D_elemental_SFI!(ndoubl, pol_n, J₀⁻)
    i, _, n = @index(Global, NTuple)

    if ndoubl>1
        if mod(i, pol_n) > 2
            J₀⁻[i, 1, n] = - J₀⁻[i, 1, n]
        end
    end
end

function apply_D_matrix_elemental!(ndoubl::Int, n_stokes::Int, r⁻⁺::CuArray{FT,3}, t⁺⁺::CuArray{FT,3}, r⁺⁻::CuArray{FT,3}, t⁻⁻::CuArray{FT,3}) where {FT}
    device = devi(Architectures.GPU())
    applyD_kernel! = apply_D_elemental!(device)
    event = applyD_kernel!(ndoubl,n_stokes, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, ndrange=size(r⁻⁺));
    wait(device, event);
    synchronize_if_gpu();
    return nothing
end

function apply_D_matrix_elemental!(ndoubl::Int, n_stokes::Int, r⁻⁺::Array{FT,3}, t⁺⁺::Array{FT,3}, r⁺⁻::Array{FT,3}, t⁻⁻::Array{FT,3}) where {FT}
    device = devi(Architectures.CPU())
    applyD_kernel! = apply_D_elemental!(device)
    event = applyD_kernel!(ndoubl,n_stokes, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻, ndrange=size(r⁻⁺));
    wait(device, event);
    return nothing
end

function apply_D_matrix_elemental_SFI!(ndoubl::Int, n_stokes::Int, J₀⁻::CuArray{FT,3}) where {FT}
    if ndoubl > 1
        return nothing
    else
        device = devi(Architectures.GPU())
        applyD_kernel! = apply_D_elemental_SFI!(device)
        event = applyD_kernel!(ndoubl,n_stokes, J₀⁻, ndrange=size(J₀⁻));
        wait(device, event);
        synchronize();
        return nothing
    end
end

function apply_D_matrix_elemental_SFI!(ndoubl::Int, n_stokes::Int, J₀⁻::Array{FT,3}) where {FT}
    if ndoubl > 1
        return nothing
    else
        device = devi(Architectures.CPU())
        applyD_kernel! = apply_D_elemental_SFI!(device)
        event = applyD_kernel!(ndoubl,n_stokes, J₀⁻, ndrange=size(J₀⁻));
        wait(device, event);
        return nothing
    end
end


text:
using DiscretePredictors
using Test

macro test_nothrow(ex)
    quote
        try
            $(esc(ex))
            true
        catch e
            print( "ERROR: " )
            if isa( e, KeyError )
                println( "KeyError: key ", e.key, " not found" )
            else
                println( e )
            end
            false
        end
    end
end

tests = [
            "AdaptiveMPP",
            "ALZ",
            "DG",
            "dHedgePPM",
            "ewPPM",
            "KOM",
            "LeZiUpdate",
            "LZ78"
        ]

if length(ARGS) > 0
    tests = ARGS
end

for t in tests
    fp = joinpath( dirname(@__FILE__), "$t.jl" )
    println( "Testing : ", t )
    include( fp )
end


text:
@inline get_data(A::AbstractArray) = A
@inline get_data(A::SizedArray) = A.data

"""
	AbstractConstraintValues{C<:AbstractConstraint}

An abstract type for working with and storing constraint values, such as current
constraint values, Jacobians, dual variables, penalty parameters, etc.
The information that is actually store, and the way it is stored, is up to the
child type. However, at a minimum, it should store the following fields:
* `con::AbstractConstraint`: the actual constraint
* `vals::AbstractVector{<:AbstractVector}`: stores the constraint value for all time indices.
* `jac::AbstractMatrix{<:AbstractMatrix}`: stores the constraint Jacobian for all time indices.
* `inds::AbstractVector{Int}`: stores the time step indices.
The first dimension of all of these data fields should be the same
(the number time indices).

With these fields, the following methods are implemented:
* `evaluate!(::AbstractConstraintValues, ::AbstractTrajectory)`
* `jacobian!(::AbstractConstraintValues, ::AbstractTrajectory)`
* `max_violation(::AbstractConstraintValues)`
"""
abstract type AbstractConstraintValues{C<:AbstractConstraint} end

"""
	ConVal{C,V,M,W}

Holds information about a constraint of type `C`. Allows for any type of
vector (`V`) or matrix (`M`) storage for constraint values and Jacobians
(allowing StaticArrays or views into a large, sparse matrix).
"""
struct ConVal{C,V,M,W} <: AbstractConstraintValues{C}
    con::C
    inds::Vector{Int}
    vals::Vector{V}
	vals2::Vector{V}
    jac::Matrix{M}
    ∇x::Matrix{W}
    ∇u::Matrix{W}
    c_max::Vector{Float64}
	is_const::BitArray{2}
	const_hess::BitVector
	iserr::Bool  # are the Jacobians on the error state
	function ConVal(n::Int, m::Int, con::AbstractConstraint, inds::Union{Vector{Int},UnitRange},
			jac, vals, iserr::Bool=false)
		if !iserr && size(gen_jacobian(con)) != size(jac[1])
			throw(DimensionMismatch("size of jac[i] $(size(jac[1])) does not match the expected size of $(size(gen_jacobian(con)))"))
		end
		vals2 = deepcopy(vals)
        p = length(con)
        P = length(vals)
        ix = 1:n
        iu = n .+ (1:m)
		views = [gen_views(∇c, con, n, m) for ∇c in jac]
		∇x = [v[1] for v in views]
		∇u = [v[2] for v in views]
        c_max = zeros(P)
		is_const = BitArray(undef, size(jac))
		const_hess = BitVector(undef, P)
        new{typeof(con), eltype(vals), eltype(jac), eltype(∇x)}(con,
			collect(inds), vals, vals2, jac, ∇x, ∇u, c_max, is_const, const_hess, iserr)
    end
end

function ConVal(n::Int, m::Int, cval::ConVal)
	# create a ConVal for the "raw" Jacobians, if needed
	# 	otherwise return the same ConVal
	if cval.iserr
		p = length(cval.con)
		ws = widths(cval.con, n, m)
		jac = [SizedMatrix{p,w}(zeros(p,w)) for k in cval.inds, w in ws]
		ConVal(n, m, cval.con, cval.inds, jac, cval.vals, false)
	else
		return cval
	end
end

function ConVal(n::Int, m::Int, con::AbstractConstraint, inds::UnitRange{Int}, iserr::Bool=false)
	C,c = gen_convals(n,m,con,inds)
	ConVal(n, m, con, inds, C, c)
end

function _index(cval::AbstractConstraintValues, k::Int)
	if k ∈ cval.inds
		return k - cval.inds[1] + 1
	else
		return 0
	end
end

Base.length(cval::AbstractConstraintValues) = length(cval.con)

function evaluate!(cval::AbstractConstraintValues, Z::AbstractTrajectory)
	evaluate!(cval.vals, cval.con, Z, cval.inds)
end

function jacobian!(cval::AbstractConstraintValues, Z::AbstractTrajectory, init::Bool=false)
	if cval.iserr
		throw(ErrorException("Can't evaluate Jacobians directly on the error state Jacobians"))
	else
		if init
			cval.is_const .= true
		elseif all(cval.is_const)
			return
		end
		jacobian!(cval.jac, cval.con, Z, cval.inds, cval.is_const)
	end
end

function ∇jacobian!(G, cval::AbstractConstraintValues, Z::AbstractTrajectory, λ, init::Bool=false)
	∇jacobian!(G, cval.con, Z, λ, cval.inds, cval.const_hess, init)
end

violation(::Equality, x) = x
violation(cone::Conic, x) = projection(cone, x) - x
∇violation!(::Equality, ∇v, ∇c, x, tmp=zeros(length(x), length(x))) = ∇v .= ∇c
function ∇violation!(cone::Conic, ∇v, ∇c, x, tmp=zeros(length(x), length(x)))
	∇projection!(cone, tmp, x)
	for i = 1:length(x)
		tmp[i,i] -= 1
	end
	mul!(∇v, tmp, ∇c)
end

@inline max_violation(::Equality, v) = norm(v,Inf)
@inline max_violation(::Inequality, v) = max(0,maximum(v))
@inline max_violation(::Equality, v::Real) = abs(v)
@inline max_violation(::Inequality, v::Real) = v > 0 ? v : 0.0

function max_violation(cone::SecondOrderCone, x)
	proj = projection(cone, x)
	return norm(x - proj,Inf)
end

function max_violation(cval::AbstractConstraintValues)
	max_violation!(cval)
    return maximum(cval.c_max)
end

function max_violation!(cval::AbstractConstraintValues)
	s = sense(cval.con)
	P = length(cval)
    for i in eachindex(cval.inds)
        cval.c_max[i] = max_violation(s, SVector{P}(cval.vals[i]))
    end
end

@inline norm_violation(::Equality, v, p=2) = norm(v,p)

@inline function norm_violation(::Inequality, v, p=2)
	# TODO: try this with LazyArrays?
	if p == 1
		a = zero(eltype(v))
		for x in v
			a += max(x,0)
		end
		return a
	elseif p == 2
		a = zero(eltype(v))
		for x in v
			a += max(x, 0)^2
		end
		return sqrt(a)
	elseif p == Inf
		return maximum(v)
	else
		throw(ArgumentError("$p is not a valid norm value. Must be 1,2 or Inf"))
	end
end

function norm_violation(cval::AbstractConstraintValues, p=2)
	norm_violation!(cval, p)
	return norm(cval.c_max, p)
end

function norm_violation!(cval::AbstractConstraintValues, p=2)
	s = sense(cval.con)
	for i in eachindex(cval.inds)
		cval.c_max[i] = norm_violation(s, cval.vals[i], p)
	end
end

function norm_dgrad!(cval::AbstractConstraintValues, Z::AbstractTrajectory, p=1)
	for (i,k) in enumerate(cval.inds)
		zs = RobotDynamics.get_z(cval.con, Z, k)
		mul!(cval.vals2[i], cval.jac[i,1], zs[1])
		if length(zs) > 1
			mul!(cval.vals2[i], cval.jac[i,2], zs[2], 1.0, 1.0)
		end
		cval.c_max[i] = norm_dgrad(cval.vals[i], cval.vals2[i], p)
	end
	return nothing
end
"""
	dgrad(x, dx, p=1)
Directional derivative of `norm(x, p)` in the direction `dx`
"""
function norm_dgrad(x, dx, p=1)
	g = zero(eltype(x))
	if p == 1
		@assert length(x) == length(dx)
		g = zero(eltype(x))
		for i in eachindex(x)
			if x[i] < 0
				g += -dx[i]
			elseif x[i] > 0
				g += dx[i]
			else
				g += abs(dx[i])
			end
		end
	else
		throw("Directional derivative of $p-norm isn't implemented yet")
	end
	return g
end

function norm_residual!(res, cval::AbstractConstraintValues, λ::Vector{<:AbstractVector}, p=2)
	for (i,k) in enumerate(cval.inds)
		mul!(res[i], cval.jac[i,1], λ[i])
		if size(cval.jac,2) > 1
			mul!(res[i], cval.jac[i,2], λ[i], 1.0, 1.0)
		end
		cval.c_max[i] = norm(res[i], p)
	end
	return nothing
end

function error_expansion!(errval::AbstractConstraintValues{C}, conval::AbstractConstraintValues,
		model::AbstractModel, G) where C
	if errval.jac !== conval.jac
		n,m = size(model)
		ix = 1:n
		iu = 1:m
		for (i,k) in enumerate(conval.inds)
			if C <: StateConstraint
				mul!(errval.jac[i], conval.jac[i], get_data(G[k]))
			elseif C <: ControlConstraint
				errval.jac[i] .= conval.jac[i]
			else
				∇x  = view(errval.jac[i], :, ix)
				∇u  = view(errval.jac[i], :, iu)
				∇x0 = view(conval.jac[i], :, ix)
				∇u0 = view(conval.jac[i], :, iu)
				mul!(∇x, ∇x0, get_data(G[k]))
				∇u .= ∇u0
			end
		end
	end
end

function error_expansion!(con::AbstractConstraint, err, jac, G)
	mul!(err, jac, G)
end

function gen_convals(n̄::Int, m::Int, con::AbstractConstraint, inds)
    # n is the state diff size
    p = length(con)
	ws = widths(con, n̄,m)
    C = [SizedMatrix{p,w}(zeros(p,w)) for k in inds, w in ws]
	c = [@MVector zeros(p) for k in inds]
    return C, c
end

function gen_convals(D::AbstractMatrix, d::AbstractVector, cinds, zinds, con::AbstractConstraint, inds)
    P = length(inds)
    p = length(con)
	n,m = get_dims(con, length(zinds[1]))
    ws = widths(con, n, m)

    C = [begin
		view(D, cinds[i], zinds[k+(j-1)][1:ws[j]])
	end for (i,k) in enumerate(inds), j = 1:length(ws)]
    c = [view(d, cinds[i]) for i in 1:P]
    return C,c
end

function gen_convals(blocks::Vector, cinds, con::AbstractConstraint, inds)
	# assumes all cinds are contiguous indices (i.e. can be represented as a UnitRange)
    C1 = map(enumerate(inds)) do (i,k)
        nm = size(blocks[k].Y,2)
		if con isa StateConstraint
			iz = 1:width(con)
		elseif con isa ControlConstraint
			m = control_dim(con)
			n = nm - m
			iz = n .+ (1:m)
		else
			iz = 1:nm
		end
		ic = cinds[i][1]:cinds[i][end]
		n1 = size(blocks[k].D2, 1)
        view(blocks[k].Y, n1 .+ (ic), iz)
    end
	C2 = map(enumerate(inds)) do (i,k)
		if con isa StageConstraint
			w = size(blocks[k].Y,2)
			view(blocks[k].Y,1:0,1:w)
		else
			w = size(blocks[k+1].Y,2)
			n = state_dim(con)
			view(blocks[k+1].Y,1:n,1:w)
		end
	end
	C = [C1 C2]
    c = map(enumerate(inds)) do (i,k)
		ic = cinds[i][1]:cinds[i][end]
        view(blocks[k].y, ic)
    end
    return C,c
end


text:
# conclusion the array_board moves are very slow compared to bitboard

using Game2048
using Game2048: Bitboard, initboard
using BenchmarkTools

@benchmark move(board, LEFT) setup=(board=Bitboard(rand(UInt64)))
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     2.000 ns (0.00% GC)
#   median time:      2.100 ns (0.00% GC)
#   mean time:        2.514 ns (0.00% GC)
#   maximum time:     69.200 ns (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1000

@benchmark move(board, left) setup=(board=Bitboard(rand(UInt64)))
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     3.600 ns (0.00% GC)
#   median time:      3.700 ns (0.00% GC)
#   mean time:        4.308 ns (0.00% GC)
#   maximum time:     30.900 ns (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1000

@benchmark move!(board, left) setup=(board=initboard())
# BenchmarkTools.Trial:
#   memory estimate:  48 bytes
#   allocs estimate:  2
#   --------------
#   minimum time:     245.036 ns (0.00% GC)
#   median time:      279.661 ns (0.00% GC)
#   mean time:        292.479 ns (0.71% GC)
#   maximum time:     5.629 μs (93.97% GC)
#   --------------
#   samples:          10000
#   evals/sample:     413

@benchmark move(board, right) setup=(board=Bitboard(rand(UInt64)))
# BenchmarkTools.Trial:
#   memory estimate:  0 bytes
#   allocs estimate:  0
#   --------------
#   minimum time:     3.500 ns (0.00% GC)
#   median time:      3.500 ns (0.00% GC)
#   mean time:        4.061 ns (0.00% GC)
#   maximum time:     35.900 ns (0.00% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1000

@benchmark move(board, up) setup=(board=Bitboard(rand(UInt64)))
# BenchmarkTools.Trial:
#   memory estimate:  192 bytes
#   allocs estimate:  2
#   --------------
#   minimum time:     75.154 ns (0.00% GC)
#   median time:      85.421 ns (0.00% GC)
#   mean time:        98.508 ns (8.79% GC)
#   maximum time:     2.728 μs (95.36% GC)
#   --------------
#   samples:          10000
#   evals/sample:     974

@benchmark move(board, down) setup=(board=Bitboard(rand(UInt64)))
# BenchmarkTools.Trial:
#   memory estimate:  192 bytes
#   allocs estimate:  2
#   --------------
#   minimum time:     75.334 ns (0.00% GC)
#   median time:      88.181 ns (0.00% GC)
#   mean time:        101.071 ns (8.87% GC)
#   maximum time:     2.872 μs (96.13% GC)
#   --------------
#   samples:          10000
#   evals/sample:     973


text:
using Multilane
using StatsBase
using MCTS
using JLD
using POMCP

using Gallium

cors = linspace(0,1,5)

behaviors = Dict{String,Any}()
for cor in cors
    behaviors[@sprintf("cor_%03d", 100*cor)] = standard_uniform(1.0, correlation=cor)
end

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 estimate_value=RolloutEstimator(SimpleSolver()))

pomcps = POMCPDPWSolver(
    eps=0.001,
    max_depth=20,
    c=50.0,
    tree_queries=2500,
    k_observation=4.0,
    alpha_observation=1/8,
    rollout_solver=SimpleSolver()
)

agg_up = AggressivenessUpdater(nothing, 500, 0.1, 0.1, WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5), MersenneTwister(123))

solvers = Dict{String, Any}(
    "dpw"=>dpws,
    "assume_normal"=>SingleBehaviorSolver(dpws, Multilane.NORMAL),
    "pomcp"=>MLPOMDPSolver(pomcps, agg_up),
    "mlmpc"=>MLMPCSolver(dpws, agg_up)
)

point = TestSet(rmodel_type="TargetLaneReward", vel_sigma=0.25, p_appear=1.0, brake_threshold=4.0, brake_terminate_thresh=4.0, N=500)

tests = []
for cor in cors
    bg = @sprintf("cor_%03d", 100*cor)
    push!(tests, TestSet(point, solver_key="dpw", key=@sprintf("upper_bound_%03d", 100*cor), behaviors=bg))
    push!(tests, TestSet(point, solver_key="assume_normal", key=@sprintf("assume_normal_%03d", 100*cor), behaviors=bg))
    push!(tests, TestSet(point, solver_key="pomcp", key=@sprintf("pomcp_%03d", 100*cor), behaviors=bg))
    push!(tests, TestSet(point, solver_key="mlmpc", key=@sprintf("mlmpc_%03d", 100*cor), behaviors=bg))
end

objects = gen_initials(tests, behaviors=behaviors, generate_physical=true)

println(objects["param_table"])
objects["solvers"] = solvers
objects["note"] = "pomcp vs mlmpc with various levels of correlation"

files = sbatch_spawn(tests, objects,
                     batch_size=50,
                     time_per_batch="1:00:00",
                     submit_command="sbatch",
                     template_name="sherlock.sh")

# files = sbatch_spawn(tests, objects,
#                      batch_size=4,
#                      time_per_batch="1:00:00",
#                      submit_command="bash",
#                      template_name="theresa.sh")

# results = evaluate(tests, objects, parallel=true)
#
# filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
# results["histories"] = nothing
# save(filename, results)
# println("results saved to $filename")


text:
abstract type ActivationFunction end

"""
    ReLU <: ActivationFunction

    (ReLU())(x) -> max.(x, 0)
"""
struct ReLU <: ActivationFunction end

"""
    Max <: ActivationFunction

    (Max())(x) -> max(maximum(x), 0)
"""
struct Max <: ActivationFunction end

"""
    Id <: ActivationFunction
Identity operator

    (Id())(x) -> x
"""
struct Id <: ActivationFunction end

"""
    Sigmoid <: ActivationFunction

    (Sigmoid())(x) -> 1 ./ (1 .+ exp.(-x))
"""
struct Sigmoid <: ActivationFunction end

"""
    Tanh <: ActivationFunction

    (Tanh())(x) -> tanh.(x)
"""
struct Tanh <: ActivationFunction end

"""
    GeneralAct <: ActivationFunction
Wrapper type for a general activation function.

### Usage
```julia
act = GeneralAct(tanh)

act(0) == tanh(0)           # true
act(10.0) == tanh(10.0)     # true
```
```julia
act = GeneralAct(x->tanh.(x))

julia> act(-2:2)
5-element Array{Float64,1}:
 -0.9640275800758169
 -0.7615941559557649
  0.0
  0.7615941559557649
  0.9640275800758169
```
"""
struct GeneralAct <: ActivationFunction
    f::Function
end

#=
TODO: consider writing our own interpolation scheme to avoid a dependency for this one thing.
Should only require a handful of functions.
Also NOTE: inherently not type stable unless parameterized.
=#
"""
    PiecewiseLinear <: ActivationFunction
Activation function that uses linear interpolation between supplied `knots`.
An extrapolation condition can be set for values outside the set of knots. Default is `Linear`.

    PiecewiseLinear(knots_x, knots_y, [extrapolation = Line()])

### Usage
```julia
kx = [0.0, 1.2, 1.7, 3.1]
ky = [0.0, 0.5, 1.0, 1.5]
act = PiecewiseLinear(kx, ky)

act(first(kx)) == first(ky) == 0.0
act(last(kx))  == last(ky)  == 1.5

act(1.0)    # 0.4166666666666667
act(-102)   # -42.5
```
```julia
act = PiecewiseLinear(kx, ky, Flat())

act(-102)   # 0.0
act(Inf)    # 1.5
```

### Extrapolations
- Flat()
- Line()
- constant (supply a number as the argument)
- Throw() (throws bounds error)

`PiecewiseLinear` uses [Interpolations.jl](http://juliamath.github.io/Interpolations.jl/latest/).
"""
struct PiecewiseLinear <: ActivationFunction
    f::Interpolations.Extrapolation
end

# default extrapolation is Line(). Can also do Flat() or supply a constant,
function PiecewiseLinear(knots_x::AbstractVector,
                         knots_y::AbstractVector,
                         extrapolation = Interpolations.Line())
    PiecewiseLinear(LinearInterpolation(knots_x, knots_y, extrapolation_bc = extrapolation))
end

# the type stable definitions probably don't need to go in the paper as-is
(f::ReLU)(x) = max.(x, zero(eltype(x)))
(f::Max)(x) = max(maximum(x), zero(eltype(x)))
(f::Id)(x) = x
(f::Sigmoid)(x) = @. 1 / (1 + exp(-x))
(f::Tanh)(x) = tanh.(x)
(G::GeneralAct)(x) = G.f(x)
(PL::PiecewiseLinear)(x) = PL.f(x)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletrectangle

"""
    ''_leafletrectangle(;kwargs...)

A LeafletRectangle component.

Keyword arguments:
- `id` (String; optional)
- `bounds` (required): . bounds has the following type: lists containing elements 'minx', 'miny', 'maxx', 'maxy'.
Those elements have the following types:
  - `minx` (Real; required)
  - `miny` (Real; required)
  - `maxx` (Real; required)
  - `maxy` (Real; required)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `pathOptions` (optional)
"""
function ''_leafletrectangle(; kwargs...)
        available_props = Symbol[:id, :bounds, :loading_state, :pathOptions]
        wild_props = Symbol[]
        return Component("''_leafletrectangle", "LeafletRectangle", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
function generateTimeToAttrition( attrition::Attrition, nVals::Int )

    urand = rand( attrition.timeRNG, nVals )
    nSections = map( u -> findlast( u .<= attrition.gammas ), urand )
    result = map( 1:nVals ) do ii
        nSection = nSections[ii]

        if attrition.lambdas[nSection] == 0
            return +Inf
        end  # if attrition.lambdas[nSection] == 0

        attrTime = attrition.curvePoints[nSection]
        attrTime -= log( urand[ii] / attrition.gammas[nSection] ) /
            attrition.lambdas[nSection]
        return attrTime #* attrition.period
    end  # map( 1:nVals ) do ii

    return result #./ attrition.period

end  # generateTimeToAttrition( attrition, nVals )

function generateTimeToAttrition( attrition::Attrition,
    excessTimes::Vector{T} ) where T <: Real

    if isempty( excessTimes )
        return Vector{Float64}()
    end  # if isempty( excessTimes )

    eSections = map( et -> findlast( attrition.curvePoints .<= et ),
        excessTimes )
    urand = rand( attrition.condtimeRNG, length( excessTimes ) )
    urand .*= attrition.gammas[eSections]

    result = map( eachindex( excessTimes ) ) do ii
        eSection = eSections[ii]
        condprob = exp( -attrition.lambdas[eSection] *
            ( excessTimes[ii] - attrition.curvePoints[eSection] ) )
        nSection = findlast( urand[ii] * condprob .<= attrition.gammas )

        if attrition.lambdas[nSection] == 0
            return +Inf
        end  # if attrition.lambdas[nSection] == 0

        attrTime = attrition.curvePoints[nSection]
        attrTime += ( excessTimes[ii] - attrition.curvePoints[eSection] ) *
            attrition.lambdas[eSection] / attrition.lambdas[nSection]
        attrTime -= log( urand[ii] / attrition.gammas[nSection] ) /
            attrition.lambdas[nSection]
        return attrTime
    end  # map( eachindex( excessTimes ) ) do ii

    return result

end  # generateTimeToAttrition( attrition, excessTimes )

@resumable function checkAttritionProcess( sim::Simulation, mpSim::MPsim )

    # Timing of the process.
    processTime = Dates.Millisecond( 0 )
    tStart = now()
    attName = string( "Attrition check process " )

    # Preparatory steps.
    timeOfNextCheck = 0.0
    checkPriority = typemin(Int)
    execPriority = typemax(Int)
    # queryCmd = string( "SELECT `", mpSim.idKey, "`, expectedAttritionTime ",
    queryCmd = string( "SELECT `", mpSim.idKey, "`, expectedAttritionTime, currentNode ",
        "FROM `", mpSim.persDBname, "` WHERE",
        "\n    status IS 'active' AND expectedAttritionTime <= " )
    idSymb = Symbol( mpSim.idKey )

    while now(sim) < mpSim.simLength
        processTime += now() - tStart
        @yield( timeout( sim, timeOfNextCheck - now( sim ),
            priority=checkPriority ) )
        tStart = now()

        # Find who'll undergo attrition in the next cycle.
        # timeOfNextCheck += mpSim.attritionTimeSkip
        timeOfNextCheck = getNextEventTime(sim)
        timeOfNextCheck = timeOfNextCheck > mpSim.simLength ? mpSim.simLength :
            timeOfNextCheck
        # now(sim) == timeOfNextCheck && continue
        result = DataFrame( DBInterface.execute( mpSim.simDB,
            string( queryCmd, timeOfNextCheck ) ) )

        # Execute attrition process for these people.
        processTime += now() - tStart

        if !isempty(result)
            tStart = now()
            removePersons( result[:, idSymb], result[:, :currentNode],
                result[:, :expectedAttritionTime], "attrition", mpSim )
            mpSim.attritionExecTime += now() - tStart
        end  # if !isempty(result)
#=
        for ii in 1:size( result, 1 )
            @process executeAttritionProcess( sim, result[ii, idSymb],
                result[ii, :expectedAttritionTime], mpSim, execPriority )
        end  # for ii in 1:size( result, 1 )
=#
        tStart = now()
    end  # while now( sim ) < mpSim.simLength

    # Timing of the process.
    processTime += now() - tStart

    if mpSim.showInfo
        println( attName, "took ", processTime.value / 1000,
            " seconds." )
    end  # if mpSim.showInfo

end  # checkAttritionProcess( sim, mpSim )

@resumable function executeAttritionProcess( sim::Simulation, id::String,
    expAttrTime::Float64, mpSim::MPsim, priority::Int )

    # Wait until the time attrition should take place.
    @yield( timeout( sim, expAttrTime - now( sim ), priority=priority ) )
    tStart = now()

    # Get the personnel record.
    queryCmd = string( "SELECT expectedAttritionTime, status, currentNode ",
        "FROM `", mpSim.persDBname, "`  WHERE",
        "\n    `", mpSim.idKey, "` IS '", id, "'" )
    persRecord = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

    # Only perform the attrition if the person is still active and the time of
    #   attrition hasn't changed.
    if ( persRecord[1, :status] == "active" ) &&
        ( persRecord[1, :expectedAttritionTime] == expAttrTime )
        removePerson( id, persRecord[1, :currentNode], "attrition", mpSim )
    end  # if ( persRecord[1, :status] == "active" ) && ...

    mpSim.attritionExecTime += now() - tStart

end  # executeAttritionProcess( sim, id, expAttrTime, mpSim )

getNextEventTime( sim::Simulation ) = first(sim.heap).second.time


text:
"""
    minimalsetup1d(N)

Creates a minimal setup of all important variables.

Such it generates the moments `height`, `velocity` as well as a dummy distribution with its speeds `f0`, `f1` and `f2`.
Further a set of forcing is supplied, they can be used with `foces.slip`, `.h∇p`, `.bathymetry` and `.thermal`.

# Example
```jldoctest
julia> using JuSwalbe

julia> constants, mom, f, distribution = minimalsetup1d(10)
(JuSwalbe.Inputconstants
  lx: Int64 10
  ly: Int64 512
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0
, JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  velocity: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  pressure: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
, JuSwalbe.Forces{Array{Float64,1}}
  slip: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  h∇p: Array{Float64}((10,)) [-0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1, -0.1]
  bathymetry: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  thermal: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
, JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  f1: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  f2: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)

julia> mom.height
10-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
```
See also: [`minimalsetup2d`](@ref)
"""
function minimalsetup1d(N::Int)
    input = Inputconstants(lx=N)
    mom = Macroquant(height = ones(N), velocity=fill(0.1, N), pressure=zeros(N), energy=zeros(N))
    forces = Forces(slip=fill(0.1, N), thermal=zeros(N), h∇p=fill(-0.1, N), bathymetry=zeros(N))
    dist = DistributionD1Q3(f0 = ones(N), f1 = zeros(N), f2 = zeros(N))
    return input, mom, forces, dist
end

"""
    minimalsetup2d(N,M)

Creates a minimal working example set of macroscopic moments, forces and a distribution.

# Example
```jldoctest
julia> using JuSwalbe

julia> input, mom, f, dist = minimalsetup2d(5,5)
(JuSwalbe.Inputconstants
  lx: Int64 5
  ly: Int64 5
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0
, JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
, JuSwalbe.Forces{JuSwalbe.Twovector{Array{Float64,2}}}
  slip: JuSwalbe.Twovector{Array{Float64,2}}
  h∇p: JuSwalbe.Twovector{Array{Float64,2}}
  bathymetry: JuSwalbe.Twovector{Array{Float64,2}}
  thermal: JuSwalbe.Twovector{Array{Float64,2}}
, JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f2: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f3: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f4: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f5: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f6: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f7: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  f8: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
)

```

See also: [`minimalsetup1d`](@ref)
"""
function minimalsetup2d(N::Int, M::Int)
    input = Inputconstants(lx=N, ly=M)
    vel = Twovector(x=fill(0.02, (N,M)), y=fill(-0.02, (N,M)))
    zerovec = Twovector(x=zeros(N,M), y=zeros(N,M))
    mom = Macroquant(height = ones(N,M), velocity=vel, pressure=zeros(N,M), energy=zeros(N,M))
    forces = Forces(slip=vel, thermal=zerovec, h∇p=vel, bathymetry=zerovec)
    dist = DistributionD2Q9(f0 = ones(N,M), f1 = zeros(N,M), f2 = zeros(N,M),
                            f3 = zeros(N,M), f4 = zeros(N,M), f5 = zeros(N,M),
                            f6 = zeros(N,M), f7 = zeros(N,M), f8 = zeros(N,M))
    return input, mom, forces, dist
end

"""
  simplemoment2d(n,m,type)

Creates moments of type JuSwalbe.Macroquant with dimensions (n,m).

# Example
```jldoctest
julia> using JuSwalbe

julia> mom = simplemoment2d(5,5)
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> mom.height
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

```
"""
function simplemoment2d(n::Int, m::Int; T=Float64)
  mom = JuSwalbe.Macroquant(height=ones(T, (n,m)),
                            velocity=JuSwalbe.Twovector(x=fill(T(0.1),(n,m)), y=fill(T(-0.1),(n,m))),
                            pressure=zeros(T, (n,m)),
                            energy=zeros(T,(n,m)))
  return mom
end

"""
  simplemoment1d(n,type)

Creates moments of type JuSwalbe.Macroquant with length n.

# Example
```jldoctest
julia> using JuSwalbe

julia> mom = simplemoment1d(5)
JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((5,)) [1.0, 1.0, 1.0, 1.0, 1.0]
  velocity: Array{Float64}((5,)) [0.1, 0.1, 0.1, 0.1, 0.1]
  pressure: Array{Float64}((5,)) [0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((5,)) [0.0, 0.0, 0.0, 0.0, 0.0]

julia> mom.height
5-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0

```
"""
function simplemoment1d(n::Int; T=Float64)
  mom = JuSwalbe.Macroquant(height=ones(T, n), velocity=fill(T(0.1),n), pressure=zeros(T, n), energy=zeros(T, n))
  return mom
end

"""
  simpleTwovector(n,m,T)

Generates a JuSwalbe.Twovector of type T and dimension (n,m).

# Example
```jldoctest
julia> using JuSwalbe

julia> twovec = simpleTwovector(5,5)
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  y: Array{Float64}((5, 5)) [-0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1; … ; -0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1]

julia> twovec.x
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
```
"""
function simpleTwovector(n::Int,m::Int; T=Float64)
  xy = JuSwalbe.Twovector(x=fill(T(0.1),(n,m)), y=fill(T(-0.1),(n,m)))
  return xy
end

function zeroTwovector(n::Int,m::Int; T=Float64)
  xy = JuSwalbe.Twovector(x=zeros(T ,(n,m)), y=zeros(T, (n,m)))
  return xy
end

function zeroTwovectorCU(n::Int,m::Int; T=Float32)
  jaja = zeros(n,m)
  xy = JuSwalbe.Twovector(x=CuArray(jaja), y=CuArray(jaja))
  return xy
end

"""
  simpledistD2Q9(n, m, T)

Generates a `D2Q9` distribution function of dimension `n`,`m` with defined fill statements.

# Example
```jldoctest
julia> using JuSwalbe

julia> dist = simpledistD2Q9(5,5)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  f1: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f2: Array{Float64}((5, 5)) [-0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1; … ; -0.1 -0.1 … -0.1 -0.1; -0.1 -0.1 … -0.1 -0.1]
  f3: Array{Float64}((5, 5)) [0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01; … ; 0.01 0.01 … 0.01 0.01; 0.01 0.01 … 0.01 0.01]
  f4: Array{Float64}((5, 5)) [-0.01 -0.01 … -0.01 -0.01; -0.01 -0.01 … -0.01 -0.01; … ; -0.01 -0.01 … -0.01 -0.01; -0.01 -0.01 … -0.01 -0.01]
  f5: Array{Float64}((5, 5)) [0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2; … ; 0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2]
  f6: Array{Float64}((5, 5)) [0.02 0.02 … 0.02 0.02; 0.02 0.02 … 0.02 0.02; … ; 0.02 0.02 … 0.02 0.02; 0.02 0.02 … 0.02 0.02]
  f7: Array{Float64}((5, 5)) [0.02 0.02 … 0.02 0.02; 0.02 0.02 … 0.02 0.02; … ; 0.02 0.02 … 0.02 0.02; 0.02 0.02 … 0.02 0.02]
  f8: Array{Float64}((5, 5)) [0.5 0.5 … 0.5 0.5; 0.5 0.5 … 0.5 0.5; … ; 0.5 0.5 … 0.5 0.5; 0.5 0.5 … 0.5 0.5]

julia> dist.f0
5×5 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> typeof(dist)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
```

See also: [`simpleTwovector`](@ref), [`simplemoment2d`](@ref)

"""
function simpledistD2Q9(n::Int, m::Int; T=Float64)
  dist = JuSwalbe.DistributionD2Q9(f0 = ones(T, (n,m)), f1 = fill(T(0.1), (n,m)), f2 = fill(T(-0.1), (n,m)),
                                   f3 = fill(T(0.01), (n,m)), f4 = fill(T(-0.01), (n,m)), f5 = fill(T(0.2), (n,m)),
                                   f6 = fill(T(0.02), (n,m)), f7 = fill(T(0.02), (n,m)), f8 = fill(T(0.5), (n,m)))

  return dist
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    ydn2md(year, day) -> date

### Purpose ###

Convert from year and day number of year to a date.

### Explanation ###

Returns the date corresponding to the `day` of `year`.

### Arguments ###

* `year`: the year, as a scalar integer.
* `day`: the day of `year`, as an integer.  It is can be either a scalar or
  array of integers.

### Output ###

The date, of `Date` type, of \$\\text{day} - 1\$ days after January 1st of
`year`.

### Example ###

Find the date of the 60th and 234th days of the year 2016.

``` julia
ydn2md(2016, [60, 234])
# => 2-element Array{Date,1}:
#     2016-02-29
#     2016-08-21
```

### Note ###

`ymd2dn` converts from a date to day of the year.
"""
function ydn2md(year::Integer, day::Integer)
    return Dates.firstdayofyear(Date(year)) + Dates.Day(day - 1)
end

function ydn2md{D<:Integer}(year::Integer, days::AbstractArray{D})
    dates = similar(days, Date)
    for d in eachindex(days)
        dates[d] = ydn2md(year, days[d])
    end
    return dates
end


text:
################################################################################
# convenient function for CategoricalArrays
################################################################################
import SortingLab:sorttwo!
using SortingLab
import StatsBase: rle
using CategoricalArrays

SortingLab.sorttwo!(x::CategoricalVector, y) = begin
    SortingLab.sorttwo!(x.refs, y)
    x, y
end

pooltype(::CategoricalPool{T,S}) where {T, S} = T,S

rle(x::CategoricalVector) = begin
   	refrle = rle(x.refs)
   	T,S = pooltype(x.pool)
   	(CategoricalArray{T, 1}(S.(refrle[1]), x.pool), refrle[2])
end

"""
    StringVector(v::CategoricalVector{String})

Convert `v::CategoricalVector` efficiently to WeakRefStrings.StringVector

## Example
```julia
using DataFrames
a  = categorical(["a","c", "a"])
a.refs
a.pool.index

# efficiently convert
sa = StringVector(a)

sa.buffer
sa.lengths
sa.offsets
```
"""
StringVector(v::CategoricalVector{S}) where S<:AbstractString = begin
    sa = StringVector(v.pool.index)
    StringVector{S}(sa.buffer, sa.offsets[v.refs], sa.lengths[v.refs])
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
include("builder_defs.jl")
include("builder_tools.jl")

# Ignore deploy flags
filter!(arg -> !startswith(arg, "--deploy"), ARGS)

# Allow to build specific versions
requested_version = last(BB.extract_flag!(ARGS, "--version"))
wants_version = vn -> (
    requested_version === nothing || VersionNumber(requested_version) == vn
)

# Allow to specify subsets of platforms (useful for Travis CI)
requested_targets = extract_targets!(ARGS)
wants_target = t -> (isempty(requested_targets) || t in requested_targets)

#==============================#
#  Collect selected platforms  #
#==============================#
platforms = Platform[]
if wants_target("x86_64-linux-gnu")
    push!(platforms, Linux(:x86_64, libc=:glibc))
end
if wants_target("aarch64-linux-gnu")
    push!(platforms, Linux(:aarch64, libc=:glibc))
end
if wants_target("powerpc64le-linux-gnu")
    push!(platforms, Linux(:powerpc64le, libc=:glibc))
end
platforms = expand_cxxstring_abis(platforms)

#======================#
#  Build the tarballs  #
#======================#
output = Dict()

for version in keys(versions_info)
    if wants_version(version)
        tag, hash = versions_info[version]
        sources = [ArchiveSource(source_url(tag), hash)]
        output[version] = build_tarballs(
            ARGS, name, version, sources, script, platforms, products, dependencies
        )
    end
end

save(output)

max_stars_repo_name: pabloferz/LAMMPSBuilder

max_issues_repo_name: pabloferz/LAMMPSBuilder

max_forks_repo_name: pabloferz/LAMMPSBuilder