text:
#///////////////////////////////////////
#// File Name: utils.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/02/25
#// Description: Utility functions for Risk Sensitive SAC
#///////////////////////////////////////

import StatsBase: fit, Histogram
import PlotUtils: plot_color

using Distributions
using FileIO
using JLD2
using Plots
using Printf
using ProgressMeter
using Random
using RobotOS

# convert u_nominal to u_schedule (ordereddict)
function convert_to_schedule(t_init::Time, u_nominal::Vector{Float64},
                             sim_param::SimulationParameter)
    u_schedule = OrderedDict{Time, Vector{Float64}}();
    plan_horizon = sim_param.prediction_steps*sim_param.dto;
    t = t_init;
    while t < Time(to_sec(t_init) + plan_horizon)
        u_schedule[t] = u_nominal;
        t += Duration(sim_param.dtc);
    end
    return u_schedule
end

# compute nominal trajectory
@inline function get_nominal_trajectory(e_init::RobotState,
                                e_pos_goal_vec::Vector{Float64},
                                target_speed::Float64,
                                sim_horizon::Float64,
                                prediction_horizon::Float64)
    target_horizon = norm(e_pos_goal_vec - get_position(e_init))/target_speed;
    points = [WayPoint2D(get_position(e_init), e_init.t),
              WayPoint2D(e_pos_goal_vec, e_init.t + Duration(target_horizon))];
    if target_horizon <= sim_horizon + prediction_horizon
        push!(points, WayPoint2D(e_pos_goal_vec,
                                 e_init.t + Duration(sim_horizon) +
                                            Duration(prediction_horizon)));
    end
    return Trajectory2D(points)
end

function init_condition_setup(;# Ego Initial Conditions
                              ego_pos_init_vec::Vector{Float64},
                              ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                              ego_pos_goal_vec::Vector{Float64},
                              t_init::Time=Time(0.0),
                              # Ado Initial Positions
                              ado_positions::Dict{String, Vector{Float64}},
                              # Other Parameters
                              target_speed::Union{Nothing, Float64}=nothing,
                              sim_horizon::Float64,
                              sim_param::SimulationParameter,
                              target_trajectory_required::Bool=true)
    # Initial condition setup
    if isnothing(ego_vel_init_vec)
        @assert !isnothing(target_speed)
        ego_vel_init_vec = (ego_pos_goal_vec - ego_pos_init_vec)./
                           norm(ego_pos_goal_vec - ego_pos_init_vec).*
                           target_speed;
    end
    ego_state_init_vec = vcat(ego_pos_init_vec, ego_vel_init_vec);
    e_init = RobotState(ego_state_init_vec, t_init);
    w_init = WorldState(e_init, ado_positions, t_init);

    # Measurement schedule setup
    measurement_schedule =
        get_measurement_schedule(w_init, sim_horizon, sim_param);

    if target_trajectory_required
        @assert !isnothing(target_speed)
        # Target Trajectory setup
        prediction_horizon = sim_param.dto*sim_param.prediction_steps;
        target_trajectory = get_nominal_trajectory(e_init,
                                                   ego_pos_goal_vec,
                                                   target_speed,
                                                   sim_horizon,
                                                   prediction_horizon);
        return w_init, measurement_schedule, target_trajectory
    else
        return w_init, measurement_schedule
    end
end

function controller_setup(# Scene Loader parameters
                          scene_param::TrajectronSceneParameter,
                          # Predictor Parameters
                          predictor_param::TrajectronPredictorParameter;
                          prediction_device::String,
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::ControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          # Ego Initial Conditions
                          ego_pos_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Union{Nothing, Vector{Float64}}=nothing,
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Union{Nothing, Float64}=nothing,
                          sim_horizon::Float64,
                          ado_id_to_replace::Union{Nothing, String}=nothing,
                          verbose=true)
    if verbose
        println("Scene Mode: data")
        println("Prediction Mode: trajectron")
        println("Deterministic Prediction: $(predictor_param.deterministic)")
    end
    if prediction_device == "cpu"
        device = py"torch".device("cpu");
    elseif prediction_device == "cuda"
        device = py"torch".device("cuda");
    end
    scene_loader = TrajectronSceneLoader(scene_param, verbose=verbose,
                                         ado_id_removed=ado_id_to_replace);
    if !isnothing(ado_id_to_replace)
        @assert isnothing(ego_pos_init_vec) && isnothing(ego_vel_init_vec) &&
                isnothing(ego_pos_goal_vec) && isnothing(target_speed)
        target_trajectory = get_trajectory_for_ado(scene_loader,
                                                   t_init,
                                                   ado_id_to_replace,
                                                   sim_horizon);
        t_init_new = minimum(target_trajectory)[1];
        t_end = maximum(target_trajectory)[1];
        timesteps_forward = Int64(round(to_sec(t_init_new - t_init)/scene_loader.dto, digits=5));
        println("Found $(ado_id_to_replace). start_time_idx updated to: $(scene_param.start_time_idx + timesteps_forward) Re-loading Scene...")
        sim_horizon -= to_sec(t_init_new - t_init);
        ego_pos_init_vec = minimum(target_trajectory)[2];
        ego_vel_init_vec = (collect(values(target_trajectory))[2] - ego_pos_init_vec)./scene_loader.dto;
        ego_pos_goal_vec = maximum(target_trajectory)[2];
        scene_param_new = TrajectronSceneParameter(scene_param.conf_file_name,
                                               scene_param.test_data_name,
                                               scene_param.test_scene_id,
                                               scene_param.start_time_idx + timesteps_forward,
                                               scene_param.incl_robot_node);
        scene_loader = TrajectronSceneLoader(scene_param_new, verbose=verbose,
                                             ado_id_removed=ado_id_to_replace);
        target_trajectory = get_trajectory_for_ado(scene_loader,
                                                   t_init_new,
                                                   ado_id_to_replace,
                                                   sim_horizon);
        prediction_horizon = scene_loader.dto*predictor_param.prediction_steps;
        if to_sec(t_end - t_init_new) < sim_horizon + prediction_horizon
            target_trajectory[t_init_new + Duration(sim_horizon + prediction_horizon)] = ego_pos_goal_vec;
        end
    end
    predictor = TrajectronPredictor(predictor_param,
                                    scene_loader.model_dir,
                                    scene_loader.param.conf_file_name,
                                    device, verbose=verbose);
    initialize_scene_graph!(predictor, scene_loader);
    sim_param = SimulationParameter(scene_loader, predictor, dtc, cost_param);

    ado_inputs = fetch_ado_positions!(scene_loader, return_full_state=true);
    ado_positions = reduce_to_positions(ado_inputs);

    if !isnothing(ado_id_to_replace)
        key_to_remove = nothing
        for key in keys(ado_positions)
            if pybuiltin("str")(key) == ado_id_to_replace
                key_to_remove = key
            end
        end
        delete!(ado_positions, key_to_remove)
        delete!(ado_inputs, key_to_remove)

        println("Note initial time is set to $(to_sec(t_init_new)) [s].")
        w_init, measurement_schedule =
            init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                                 ego_vel_init_vec=ego_vel_init_vec,
                                 ego_pos_goal_vec=ego_pos_goal_vec,
                                 t_init=t_init_new,
                                 ado_positions=convert_nodes_to_str(ado_positions),
                                 target_speed=target_speed,
                                 sim_horizon=sim_horizon,
                                 sim_param=sim_param,
                                 target_trajectory_required=false);
        target_speed = norm((ego_pos_goal_vec - ego_pos_init_vec)./(to_sec(t_end - t_init_new)));
    else
        @assert !isnothing(ego_pos_init_vec) && !isnothing(ego_pos_goal_vec) &&
                !isnothing(target_speed)
        w_init, measurement_schedule, target_trajectory =
            init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                                 ego_vel_init_vec=ego_vel_init_vec,
                                 ego_pos_goal_vec=ego_pos_goal_vec,
                                 t_init=t_init,
                                 ado_positions=convert_nodes_to_str(ado_positions),
                                 target_speed=target_speed,
                                 sim_horizon=sim_horizon,
                                 sim_param=sim_param,
                                 target_trajectory_required=true);
    end
    # Controller setup
    u_schedule = convert_to_schedule(w_init.t, cnt_param.u_nominal_base, sim_param);
    controller = RSSACController(predictor, u_schedule, sim_param, cnt_param);
    if predictor.param.use_robot_future
        schedule_prediction!(controller, ado_inputs, nothing, w_init.e_state);
    else
        schedule_prediction!(controller, ado_inputs);
    end
    wait(controller.prediction_task);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed
end

function controller_setup(# Scene Loader parameters
                          scene_param::TrajectronSceneParameter,
                          # Predictor Parameters
                          predictor_param::OraclePredictorParameter;
                          prediction_device::String,
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::ControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: data")
        println("Prediction Mode: oracle")
        println("Deterministic Prediction: true")
    end
    if prediction_device == "cuda"
        @warn "Prediction on CUDA devices is not currently supported. Using CPU."
    end
    device = py"torch".device("cpu");
    scene_loader = TrajectronSceneLoader(scene_param, verbose=verbose);
    predictor = OraclePredictor(predictor_param, scene_loader);
    sim_param = SimulationParameter(scene_loader, predictor, dtc, cost_param);

    ado_positions = fetch_ado_positions!(scene_loader);

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=convert_nodes_to_str(ado_positions),
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);

    # Controller setup
    u_schedule = convert_to_schedule(w_init.t, cnt_param.u_nominal_base, sim_param);
    controller = RSSACController(predictor, u_schedule, sim_param, cnt_param);
    schedule_prediction!(controller, ado_positions);
    wait(controller.prediction_task);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed
end

function controller_setup(# Scene Loader parameters
                          scene_param::TrajectronSceneParameter,
                          # Predictor Parameters
                          predictor_param::GaussianPredictorParameter;
                          prediction_device::String,
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::ControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          ado_vel_dict::Dict{String, <:MvNormal},
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: data")
        println("Prediction Mode: gaussian")
        println("Deterministic Prediction: $(predictor_param.deterministic)")
    end
    if prediction_device == "cuda"
        @warn "Prediction on CUDA devices is not currently supported. Using CPU."
    end
    device = py"torch".device("cpu");
    scene_loader = TrajectronSceneLoader(scene_param, verbose=verbose);
    @assert collect(keys(ado_vel_dict)) == ["Any"] "Set ado_vel_dict key to `Any` in this mode."
    predictor = GaussianPredictor(predictor_param, scene_loader.dto, ado_vel_dict);
    sim_param = SimulationParameter(predictor, dtc, cost_param);

    ado_positions = fetch_ado_positions!(scene_loader);

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=convert_nodes_to_str(ado_positions),
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);

    # Controller setup
    u_schedule = convert_to_schedule(w_init.t, cnt_param.u_nominal_base, sim_param);
    controller = RSSACController(predictor, u_schedule, sim_param, cnt_param);
    schedule_prediction!(controller, convert_nodes_to_str(ado_positions));
    wait(controller.prediction_task);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed
end

function controller_setup(# Scene Loader parameters
                          scene_param::SyntheticSceneParameter,
                          # Predictor Parameters
                          predictor_param::GaussianPredictorParameter;
                          prediction_device::String,
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::ControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          dto::Float64,
                          ado_pos_init_dict::Dict{String, Vector{Float64}},
                          ado_vel_dict::Dict{String, <:MvNormal},
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: synthetic")
        println("Prediction Mode: gaussian")
        println("Deterministic Prediction: $(predictor_param.deterministic)")
    end
    if prediction_device == "cuda"
        @warn "Prediction on CUDA devices is not currently supported. Using CPU."
    end
    device = py"torch".device("cpu");

    scene_loader = SyntheticSceneLoader(scene_param,
                                        deepcopy(ado_pos_init_dict));
    predictor = GaussianPredictor(predictor_param, dto, ado_vel_dict);
    sim_param = SimulationParameter(predictor, dtc, cost_param);

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=ado_pos_init_dict,
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);
    # Controller setup
    u_schedule = convert_to_schedule(w_init.t, cnt_param.u_nominal_base, sim_param);
    controller = RSSACController(predictor, u_schedule, sim_param, cnt_param);
    schedule_prediction!(controller, ado_pos_init_dict);
    wait(controller.prediction_task);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed
end

function controller_setup(# Scene Loader parameters
                          scene_param::TrajectronSceneParameter;
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::BICControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          prediction_steps::Int64,
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: data")
        println("BIC Controller");
    end
    scene_loader = TrajectronSceneLoader(scene_param, verbose=true);
    num_samples = 0;
    sim_param = SimulationParameter(dtc, scene_loader.dto, prediction_steps,
                                    num_samples, cost_param)
    ado_positions = fetch_ado_positions!(scene_loader);

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=convert_nodes_to_str(ado_positions),
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);

    # Controller setup
    controller = BICController(sim_param, cnt_param);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed
end

function controller_setup(# Scene Loader parameters
                          scene_param::SyntheticSceneParameter;
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Control Parameters
                          cnt_param::BICControlParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          dto::Float64,
                          prediction_steps::Int64, # for LQ tracker
                          num_samples::Int64, # for stochastic simulation, not for prediction
                          sim_rng::MersenneTwister, # for stochastic simulation, not for prediction
                          ado_pos_init_dict::Dict{String, Vector{Float64}},
                          ado_vel_dict::Dict{String, <:MvNormal},
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: synthetic")
        println("BIC Controller");
    end
    scene_loader = SyntheticSceneLoader(scene_param,
                                        deepcopy(ado_pos_init_dict));
    deterministic = false;
    predictor_param = GaussianPredictorParameter(prediction_steps,
                                                 num_samples, deterministic, sim_rng);
    predictor = GaussianPredictor(predictor_param, dto, ado_vel_dict);
    sim_param = SimulationParameter(predictor, dtc, cost_param);

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=ado_pos_init_dict,
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);

    # Controller setup
    controller = BICController(sim_param, cnt_param);

    return scene_loader, controller, w_init, measurement_schedule,
           target_trajectory, target_speed, predictor
end

function controller_setup(# Scene Loader parameters
                          scene_param::TrajectronSceneParameter,
                          # Control Parameters
                          cnt_param::CrowdNavControlParameter;
                          # Cost Parameters
                          cost_param::CostParameter,
                          # Simulation Parameters
                          dtc::Float64,
                          prediction_steps::Int64,
                          # Ego Initial Conditions
                          ego_pos_init_vec::Vector{Float64},
                          ego_vel_init_vec::Union{Nothing, Vector{Float64}}=nothing,
                          ego_pos_goal_vec::Vector{Float64},
                          t_init::Time=Time(0.0),
                          # Other Parameters
                          target_speed::Float64,
                          sim_horizon::Float64,
                          verbose=true)
    if verbose
        println("Scene Mode: data")
        println("CrowdNav Controller");
    end
    scene_loader = TrajectronSceneLoader(scene_param, verbose=true);
    num_samples = 0;
    sim_param = SimulationParameter(dtc, scene_loader.dto, prediction_steps,
                                    num_samples, cost_param)
    ado_inputs = fetch_ado_positions!(scene_loader, return_full_state=true);
    ado_positions = reduce_to_positions(ado_inputs)

    w_init, measurement_schedule, target_trajectory =
        init_condition_setup(ego_pos_init_vec=ego_pos_init_vec,
                             ego_vel_init_vec=ego_vel_init_vec,
                             ego_pos_goal_vec=ego_pos_goal_vec,
                             t_init=t_init,
                             ado_positions=convert_nodes_to_str(ado_positions),
                             target_speed=target_speed,
                             sim_horizon=sim_horizon,
                             sim_param=sim_param);

    # Controller setup
    controller = CrowdNavController(sim_param, cnt_param);
    @assert controller.rl_robot.v_pref == target_speed "target_speed needs to match v_pref of the CrowdNav Robot"

    return scene_loader, controller, w_init, ado_inputs, measurement_schedule,
           target_trajectory, target_speed
end

# display log
function display_log(log::Vector{Tuple{Time, String}})
    for l in log
        time = @sprintf "Time %.2f" round(to_sec(l[1]), digits=5)
        println("$(time) [sec]: $(l[2])");
    end
end

# plot helper functions
#=
function visualize!(color_dict::Dict, ap_dict::Dict{String, Vector{Float64}},
                    outputs_dict::Dict{String, Array{Float64, 3}};
                    xlim=(0., 20.), ylim=(0., 20.), markersize=8.0)
    plt = plot(xlim=xlim, ylim=ylim, aspect_ratio=1.0, xlabel="x[m]",
               ylabel="y[m]");
    for key in keys(ap_dict)
        if !haskey(color_dict, key)
            color_dict[key] = palette(:default)[length(color_dict) + 1];
        end
        scatter!((ap_dict[key][1], ap_dict[key][2]),
                 color=color_dict[key], label="", #label=key,
                 markersize=markersize/1.5, markershape=:rect);
        if haskey(outputs_dict, key)
            for jj = 1:size(outputs_dict[key], 1)
                plot!(outputs_dict[key][jj, :, 1], outputs_dict[key][jj, :, 2],
                      color=color_dict[key], label="", alpha=0.3)
            end
        end
        #jj = 5
        #plot!(outputs_array[ii, jj, :, 1], outputs_array[ii, jj, :, 2],
        #      color=color_dict[node_names_and_pos[ii][1]], label="", alpha=0.3)
    end
    return plt
end
=#

function visualize!(color_dict::Dict, w::WorldState,
                    target_trajectory::Trajectory2D,
                    outputs_dict::Dict{String, Array{Float64, 3}},
                    num_samples::Int64,
                    #nominal_control_idx::Int64,
                    nominal_trajectory::Union{Nothing, Vector{Vector{Float64}}}=nothing;
                    figsize=(600, 400), legend=:topright,
                    legendfontsize=8.,
                    xlim=(0., 20.), ylim=(0., 20.), markersize=8.0,
                    show_velocity=true,
                    show_prediction=true,
                    dummy_pos=nothing)
    plt = plot(legendfontsize=legendfontsize,
               legend=legend, size=figsize);
    plot!(xlim=xlim, ylim=ylim, aspect_ratio=1.0, xlabel="x[m]", ylabel="y[m]");
    robot_key = "Ego Robot"
    # palette = get_color_palette(:auto, plot_color(:white), 30)
    palette = get_color_palette(:auto, plot_color(:white), 150)
    if !haskey(color_dict, robot_key)
        color_dict[robot_key] = palette[1];
    end
    t = @sprintf "Time: %.2f [s]" round(to_sec(w.t), digits=2)
    scatter!((w.e_state.x[1], w.e_state.x[2]), color=color_dict[robot_key],
             label=robot_key,
             markersize=markersize,
             title=t)
    if show_velocity
        pos = get_position(w.e_state);
        vel = get_velocity(w.e_state);
        plot!([pos[1], pos[1] + vel[1]],
              [pos[2], pos[2] + vel[2]],
              color=color_dict[robot_key], label="Velocity", arrow=:arrow,
              alpha=0.5)
    end
    plot!([p[1] for p in values(target_trajectory)],
          [p[2] for p in values(target_trajectory)],
          color=color_dict[robot_key], label="Target Trajectory", linestyle=:dash, alpha=0.5)
    goal_pos = collect(values(target_trajectory))[end];
    scatter!((goal_pos[1], goal_pos[2]), color=color_dict[robot_key],
             label="Ego Goal",
             markershape=:star5,
             markersize=markersize)
    if !isnothing(nominal_trajectory)
        plot!([p[1] for p in nominal_trajectory],
              [p[2] for p in nominal_trajectory],
              color=color_dict[robot_key], label="Ego Nominal Trajectory",
              alpha=1.0, linestyle=:dashdot, linewidth=1.5)
    end
    for key in keys(w.ap_dict)
        if !haskey(color_dict, key)
            color_dict[key] = palette[length(color_dict) + 1];
        end
        scatter!((w.ap_dict[key][1], w.ap_dict[key][2]),
                 color=color_dict[key], label="", #label=key,
                 markersize=markersize/1.5, markershape=:rect);
        if show_prediction && haskey(outputs_dict, key)
            # note that prediction is dependent on the nominal control choice.
            prediction_array = outputs_dict[key]
            if isnothing(dummy_pos)
                for jj = 1:size(prediction_array, 1)
                    plot!(prediction_array[jj, :, 1], prediction_array[jj, :, 2],
                          color=color_dict[key], label="", alpha=0.3)
                end
            else
                for jj = 1:size(prediction_array, 1)
                    prediction_array_filtered = [hcat(filter(x -> x != dummy_pos, prediction_array[jj, :, 1])...);
                                                 hcat(filter(y -> y != dummy_pos, prediction_array[jj, :, 2])...)];
                    if size(prediction_array_filtered, 1) == 2 &&
                       length(prediction_array_filtered[1, :]) == length(prediction_array_filtered[2, :])
                       plot!(prediction_array_filtered[1, :], prediction_array_filtered[2, :],
                             color=color_dict[key], label="", alpha=0.3)
                    end
                end
            end
        end
        #jj = 5
        #plot!(outputs_array[ii, jj, :, 1], outputs_array[ii, jj, :, 2],
        #      color=color_dict[node_names_and_pos[ii][1]], label="", alpha=0.3)
    end
    return plt
end

function make_gif(result::Union{EvaluationResult, BICEvaluationResult, CrowdNavEvaluationResult};
                  dtplot::Float64,
                  fps::Int64,
                  figsize::Tuple{Int64, Int64},
                  legendfontsize::Int64,
                  legend::Symbol,
                  xlim::Tuple{Float64,Float64},
                  ylim::Tuple{Float64,Float64},
                  markersize::Float64,
                  filename::String,
                  show_prediction=true,
                  show_nominal_trajectory=false,
                  show_past_ego_trajectory=false,
                  dummy_pos=nothing)
    if (typeof(result) == BICEvaluationResult || typeof(result) == CrowdNavEvaluationResult) && show_prediction
        @warn "No prediction is available with BIC or CrowdNav Controller"
        show_prediction=false;
    end
    anim = Animation();
    color_dict = Dict();
    prediction_steps = result.sim_param.prediction_steps;
    dtc = result.sim_param.dtc;
    dto = result.sim_param.dto;
    plan_horizon = dto*prediction_steps;
    @assert dtplot >= dtc;
    dtplot = dtplot;
    if typeof(result) == EvaluationResult
        prediction_dict_last = nothing;
    end
    @showprogress for ii = 1:length(result.w_history)
        w = result.w_history[ii]
        if round(to_sec(w.t)/dtplot) ≈ to_sec(w.t)/dtplot
            if typeof(result) == EvaluationResult
                if isnothing(result.prediction_dict_history[ii])
                    prediction_dict = prediction_dict_last;
                else
                    prediction_dict = result.prediction_dict_history[ii];
                    prediction_dict_last = prediction_dict;
                end
                #if ii == length(result.w_history)
                #    u_nominal_idx = result.u_nominal_idx_history[ii - 1]
                #else
                #    u_nominal_idx = result.u_nominal_idx_history[ii]
                #end
            else
                prediction_dict = Dict{String, Array{Float64, 3}}();
                #u_nominal_idx = 0;
            end
            if show_nominal_trajectory
                nominal_trajectory = result.nominal_trajectory_history[ii]
            else
                nominal_trajectory = nothing
            end
            plt = visualize!(color_dict, w,
                             result.target_trajectory_history[ii],
                             prediction_dict,
                             result.sim_param.num_samples,
                             #u_nominal_idx,
                             nominal_trajectory,
                             figsize=figsize, legend=legend,
                             legendfontsize=legendfontsize,
                             xlim=xlim, ylim=ylim, markersize=markersize,
                             show_velocity=true,
                             show_prediction=show_prediction,
                             dummy_pos=dummy_pos)
            if show_past_ego_trajectory
                ego_traj_x = [get_position(w.e_state)[1] for w in result.w_history[1:ii]]
                ego_traj_y = [get_position(w.e_state)[2] for w in result.w_history[1:ii]]
                ego_color = color_dict["Ego Robot"];
                plot!(ego_traj_x, ego_traj_y, color=ego_color, label="Ego Trajectory", linewidth=1.5)
            end
            frame(anim, plt)
        end
    end
    gif(anim, filename, fps=fps);
end

function fetch_stats_filtered(func::Function...;
                              file_dir::String,
                              filtering_str::String)
    files = readdir(file_dir)
    valid_files = filter(file -> (occursin(".jld2", file) && occursin(filtering_str, file)), files)
    println("Loading $(length(valid_files)) jld2 files with '$(filtering_str)'.")

    return_array = [Real[] for ii = 1:length(func)]

    @showprogress for file in valid_files
        file_path = joinpath(file_dir, file);
        file_content = load(file_path);
        result = file_content["result"]

        for ii = 1:length(func)
            push!(return_array[ii], (func[ii])(result))
        end

        close(file_path)
        file_content = nothing;
    end
    return return_array
end

function plot_histogram(data_array...; min_val::Real, max_val::Real, num_bins::Int64,
                        color::Vector{Symbol}, alpha::Real,
                        xlabel="", ylabel="", label = ["" for ii = 1:length(data_array)],
                        title="")
    plt = plot(xlabel=xlabel, ylabel=ylabel, title=title);

    edges = range(min_val, stop=max_val, length=num_bins+1);
    offset = 0.0
    for ii = 1:length(data_array)
        hist = fit(Histogram, data_array[ii], edges .+ offset)
        delta = (edges[2] - edges[1])/length(data_array);
        bar!(hist.edges[1][1:end-1], hist.weights, bar_width=delta, alpha=alpha,
             color=color[ii], label=label[ii])
        offset += delta
    end
    #=
    for ii = 1:length(data_array)
        histogram!(data_array[ii], xlim=(min_val, max_val), bins=min_val:(max_val-min_val)/num_bins:max_val,
                    title=title, xlabel=xlabel, label=label[ii], alpha=alpha);
    end
    =#
    return plt
end


text:
using BinDeps: LibraryDependency

import BinDeps: PackageManager
import BinDeps: pkg_name
import BinDeps: provides

type Ports <: PackageManager
    package::AbstractString
end

pkg_name(p::Ports) = p.package

doc"""
    provides{S<:AbstractString}(p::Type{Ports},
                                pkgs::Vector{S},
                                deps::Vector{LibraryDependency})
"""
function provides{S<:AbstractString}(p::Type{Ports},
                                     pkgs::Vector{S},
                                     deps::Vector{LibraryDependency};
                                     opts...)
    if length(pkgs) ≠ length(deps)
        error("length of pkgs and deps mismatch")
    end

    map(zip(pkgs, deps)) do pair
        provides(p, pair...; opts...)
    end
end

ext: jl


text:
function average_linear_density_x(file, species)
    n = file["number_density/"*species]

    ω = 2π*c_0/λ
    unit_length = c_0 / ω

    x = uconvert.(NoUnits, n.grid[1]./unit_length)
    y = uconvert.(NoUnits, n.grid[2] ./ unit_length)
    z = uconvert.(NoUnits, n.grid[3] ./ unit_length)
    n_normed = uconvert.(NoUnits, n.data .* unit_length^3)
    nlx = ThreadsX.map(axes(n.data, 1)) do i
        integrate((y, z), view(n_normed, i,:,:))
    end

    x, nlx
end

function average_radial_density(file, λ)
    blocks = file_summary(file)

    n = open(file) do f
        SDFReader.read_scalar_field(f, blocks, "number_density/electron")
    end

    ω = 2π*c_0/λ
    unit_length = c_0 / ω

    x = uconvert.(NoUnits, n.grid[1]./unit_length)
    y = uconvert.(NoUnits, n.grid[2] ./ unit_length)
    z = uconvert.(NoUnits, n.grid[3] ./ unit_length)
    r = SVector.(Iterators.product(x,y,z))
    n_normed = uconvert.(NoUnits, n.data .* unit_length^3)

    coords = Cylindrical.(SVector{3}.(y,z,x))
    # n_cyl =

    nlx = ThreadsX.map(axes(n.data, 1)) do i
        integrate((y, z), view(n_normed, i,:,:))
    end

    x, nlx
end


text:
export Ellipsoid

struct Ellipsoid{K,K²} <: AbstractSimpleShape{K}
    c::SFloat{K}  # center of ellipsoid
    r::SFloat{K}  # semiaxes ("radii") in axis directions
    p::S²Float{K,K²}  # projection matrix to Ellipsoid coordinates; must be orthonormal (see surfpt_nearby)
    Ellipsoid{K,K²}(c,r,p) where {K,K²} = new(c,r,p)  # suppress default outer constructor
end

function Ellipsoid(c::AbsVecReal,  # center of ellipsoid
                   r::AbsVecReal,  # semiaxes ("radii")
                   axes::AbsMatReal=S²Float{length(c)}(I))  # columns are axes vector; assumed orthogonal
    K = length(c)
    l_ax = SVec(ntuple(k->norm(@view(axes[:,k])), Val(K)))
    p = axes' ./ l_ax

    return Ellipsoid{K,K*K}(c,r,p)
end

# Ellipsoid(s::Cuboid{K,K²}) where {K,K²} = Ellipsoid{K,K²}(s.c, (s.r).^-2, s.p)

Base.:(==)(s1::Ellipsoid, s2::Ellipsoid) = s1.c==s2.c && s1.r==s2.r && s1.p==s2.p
Base.isapprox(s1::Ellipsoid, s2::Ellipsoid) = s1.c≈s2.c && s1.r≈s2.r && s1.p≈s2.p
Base.hash(s::Ellipsoid, h::UInt) = hash(s.c, hash(s.r, hash(s.p, hash(:Ellipsoid, h))))

level(x::SReal{K}, s::Ellipsoid{K}) where {K} = norm((s.p * (x-s.c)) ./ s.r) - 1.0
# level(x::SReal{K}, s::Ellipsoid{K}) where {K} = √sum(abs2,((s.p * (x-s.c)) ./ s.r)) - 1.0
center(s::Ellipsoid) = s.c

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
# ----------------------------------------------------------------
# The required packages
using Printf
using Plots
using ProgressBars
#----------------------------------------------------------------

# ------------------------------------------------------------------------------
# To download packages from this scripts directly uncomment the following lines:
using Pkg
Pkg.add("Plots")
Pkg.add("ProgressBars")
# ------------------------------------------------------------------------------

#------------------------------------------------------------------------------
# lets define the constants
const L       = 50             # Define size of lattice
const n_sweep = 20             # Define number of sweeps between sampling
const n_therm = 100       # Define number of sweeps
const n_data  = 10000            # Define number of data samples per temperature
const temps   = 4.0:-0.3:0.1   # Define temperatures range
# ------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------------
# This function measures the i'th sample of energy and magnetization
function measure(i, energy, magnetization, s)
    en = 0
    m = 0
    for x = 1:L
        for y = 1:L
            u = 1+mod(y,L) # up
            r = 1+mod(x,L) # right
            en -= s[x,y]*(s[x,u]+s[r,y]) # energy
            m  += s[x,y]                 # magnetization
        end
    end
    energy[i] = en
    magnetization[i] = abs(m)
end
# ---------------------------------------------------------------------------------------------------------------------

#
# ----------------------------------------------------------------------------------------
 # implementing flip function to apply metropolis spin flip algorithm to site (x,y) / temp T
function flip(x, y, T, s)
    u = 1+mod(y,L)   # up
    d = 1+mod(y-2,L) # down
    r = 1+mod(x,L)   # right
    l = 1+mod(x-2,L) # left
    de = 2*s[x,y]*(s[x,u]+s[x,d]+s[l,y]+s[r,y])
    if (de < 0)
        s[x,y] = -s[x,y]
    else
        p = rand()
        if (p < exp(-de/T))
            s[x,y] = -s[x,y]
        end
    end
end
# ----------------------------------------------------------------------------------------

# --------------------------------------------------------------------
# Implementing sweeb function to apply  flip function to every site on the lattice
function sweep(n, T, s)
    for i = 1:n
        for x = 1:L
            for y = 1:L
                flip(x,y,T, s)
            end
        end
    end
end
# ----------------------------------------------------------------------------------------

# -----------------------------------------------------------------------------------------------
# Implement main julia function for algorithm opreating and calculation of physical quantities
function main()
    e1 = Array(1:n_data)     # This array purpose is to hold energy measurements (fixed T)
    m1 = Array(1:n_data)     # This array purpose is to hold magnetization measurements (fixed T)
    en = []                  # This array purpose is to append average energy at each T
    mz = []                  # This is magnetization array
    s  = ones(Int32,L,L)     # lattice of Ising spins (+/-1)
    for T in tqdm(temps)              # Lets loop over temperatures range
        sweep(n_therm, T, s)    # Sweebs over the lattice to assign tempreatures
        energy        = e1      # Reset energy measurement array
        magnetization = m1      # Reset magnetization measurement array
        for i = 1:n_data        # Take n_data measurements w/ n_sweep
            sweep(n_sweep, T, s)
            measure(i, energy, magnetization, s)
        end
        avg_en = sum(energy)/n_data           # calculate average energy
        avg_ma = sum(magnetization)/n_data    # calculate average magnetization
        push!(en,avg_en/(L*L))                # Add those to the list
        push!(mz,avg_ma/(L*L))
        # @printf("%8.3f  %8.3f \n", avg_en/(L*L), avg_ma/(L*L)) # using this for debugging code purposes to see values.
    end
# ----------------------------------------------------------------------------------------------------------------

# -------------------------------------------------------------------------------------------------------------------------------------------------
# Plotting section
    gr() # Set the backend to GR
    plot(temps,mz,title = "Magnetization vs Temperature",xlabel = "Temperature",ylabel = "Magnetization",label = "Magnetization") # plot Magnetization vs Temperature
    savefig("plots/Magnetization.pdf")
    plot(temps,en,title = "Energy vs Temperature",xlabel = "Temperature",ylabel = "Energy",label="Energy") # plot Energy vs. Temperature
    savefig("plots/Energy.pdf")
end
# ------------------------------------------------------------------------------------------------------------------------------------------------

main()

max_stars_repo_name: MohamedElashri/IsingModel

max_issues_repo_name: MohamedElashri/IsingModel

max_forks_repo_name: MohamedElashri/IsingModel


text:
using Merlin

include("bioes.jl")
include("ner.jl")
include("model.jl")

const wordembeds_file = ".data/glove.6B.100d.h5"

function main()
    argdict = Dict()
    for arg in ARGS
        kv = Vector{String}(split(arg,"="))
        length(kv) == 2 || throw("Invalid arguments: $arg")
        argdict[kv[1]] = kv[2]
    end

    if parse(Bool, argdict["train"])
        nepochs = parse(Int, argdict["nepochs"])
        #savefile = argdict["savefile"]
        ner = NER()
        traindata = readdata!(ner, ".data/eng.train.BIOES")
        testdata = readdata!(ner, ".data/eng.testb.BIOES")
        train(ner, traindata, testdata)
    else
        #seg = load("NER.jld2")
        #println(seg["a"].char2id)
        #seg = Merlin.load("NER.merlin")
    end
end
main()

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

using QuadGK #for numerical integration

"""
    C = AMFMtriplet(a, ω, φ)
    C = AMFMtriplet(a, ω)

Create a *canonical triplet* 'AMFMtriplet' consisting of an *instantenouse amplitude function* `a`, an *instantaneous frequency function* `ω`, and a *phase reference* `φ`.

Called with two inputs `a, ω`, this is equivalent to `AMFMtriplet(a, ω, 0.0)`.

# Examples
```@example
using ISA
a₀(t) = exp(-t^2)
ω₀(t) = 2.0
φ₀ = 0.0
𝐶₀ = AMFMtriplet(a₀,ω₀,φ₀)
```
"""
struct AMFMtriplet
  a::Function #real-valued function of time
  ω::Function #real-valued function of time
  φ::Real     #real-value
  function AMFMtriplet(a, ω, φ)
    if !(all(Base.return_types(a, (Int64,)).<:Real)) ; error("IA must be a real-valued function of one real variable"); end
    if !(all(Base.return_types(ω, (Int64,)).<:Real)) ; error("IF must be a real-valued function of one real variable"); end
    return new(a, ω, φ)
  end
end

# CONSTRUCTIONS
AMFMtriplet(a, ω) = AMFMtriplet(a, ω, zero(Float64))

# DISPLAY
Base.show(io::IO, x::AMFMtriplet) = print(io, "cannonical triplet")

"""
    ψ = AMFMcomp(C)
    ψ = AMFMcomp(a, ω, φ)
    ψ = AMFMcomp(a, ω)

Create a 'AMFMcomp' parameterized by a single 'AMFMtriplet'.

# Examples
```@example
 using ISA
 𝐶₀ = AMFMtriplet(t->exp(-t^2),t->2.0,0.0)
 ψ₀ = AMFMcomp(𝐶₀)
```
Another convenient way to create a 'AMFMcomp' is by providing
an *instantenouse amplitude function* `a`,
an *instantaneous frequency function* `ω`, and a *phase reference* `φ`.

Called with two inputs `a, ω`, this is equivalent to `AMFMcomp(a, ω, 0.0)`.

```@example
using ISA
a₀(t) = exp(-t^2)
ω₀(t) = 2.0
φ₀ = 0.0
ψ₀ = AMFMcomp(a₀,ω₀,φ₀)
```
"""
struct AMFMcomp
  C::AMFMtriplet
end

# CONSTRUCTIONS
AMFMcomp(a, ω, φ) = AMFMcomp(AMFMtriplet(a, ω, φ))
AMFMcomp(a, ω) = AMFMcomp(AMFMtriplet(a, ω, 0.0))

# METHODS
function (ψ::AMFMcomp)(t::Vector{<:Real})
  a, ω, φ = ψ.C.a, ψ.C.ω, ψ.C.φ
  temp = quadgk.(ω, 0.0, t)
  θ = [ temp[i][1] for i in 1:length(t) ]
  a.(t) .* exp.( 1im .* ( θ .+ φ ) )
end
function (ψ::AMFMcomp)(a::Function,ω::Function,φ::Real)
   return AMFMcomp(AMFMtriplet(a, ω, φ))
end
function (ψ::AMFMcomp)(t::Real)
  return (ψ::AMFMcomp)([t::Real])[1]
end
function (ψ::AMFMcomp)(t::StepRangeLen)
  return (ψ::AMFMcomp)(collect(t))
end
function (ψ::AMFMcomp)(t::UnitRange)
  return (ψ::AMFMcomp)(collect(t))
end

# DISPLAY
Base.show(io::IO, x::AMFMcomp) = print(io, "AM--FM component")


text:
using FreudenthalTriangulation
using Test

@testset "FreudenthalTriangulation.jl" begin
    # test freudenthal_vertices
    n, m = rand(1:10, 2)
    V = freudenthal_vertices(n, m)
    @test size(V, 1) == binomial(n + m - 1, m)
    v = rand(V)
    @test size(v, 1) == n && v[1] == m
    @test issorted(v, rev=true)

    # test freudenthal_simplex
    n = rand(5:20)
    x = rand(Float64, n)
    V = freudenthal_simplex(x)
    @test size(V, 1) == n + 1
    x = rand(Int64, n)
    x = convert(Vector{Float64}, x)
    V = freudenthal_simplex(x)
    distances = [abs(sum(x - V[i])) for i in 1:n]
    @test issorted(distances)

    # test barycentric_coordinates
    n = rand(5:20)
    x = rand(Float64, n)
    V = freudenthal_simplex(x)
    λ = barycentric_coordinates(x, V)
    @test size(λ, 1) == n + 1
    @test abs(sum(λ) - 1.0) < 1e-8

    # test freudenthal_simplex_and_coords!
    n = rand(5:20)
    x = rand(Float64, n)
    V = freudenthal_simplex(x)
    λ = barycentric_coordinates(x, V)
    V2, λ2 = freudenthal_simplex_and_coords(x)
    i = rand(1:n)
    v = V[i]
    v2 = V2[i]
    @test abs(sum(v - v2)) < 1e-8
    @test abs(sum(λ - λ2)) < 1e-8

    # test freudenthal_simplex_and_coords!
    n = rand(5:20)
    x = rand(Float64, n)
    V = freudenthal_simplex(x)
    λ = barycentric_coordinates(x, V)
    V2 = [rand(Int64, n) for i=1:n+1]
    λ2 = rand(Float64, n+1)
    V2, λ2 = freudenthal_simplex_and_coords!(x, V2, λ2)
    i = rand(1:n)
    v = V[i]
    v2 = V2[i]
    @test abs(sum(v - v2)) < 1e-8
    @test abs(sum(λ - λ2)) < 1e-8
end


text:
# ========== Transforms On A Group Of Signals ==========
# ----- SDWT on a set of signals -----
"""
    sdwtall(x, wt[, L])

Computes the stationary discrete wavelet transform (SDWT) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SDWT on all signals in x
xw = sdwtall(x, wt)
```

**See also:** [`sdwt`](@ref)
"""
function sdwtall(x::AbstractArray{T},
                 wt::OrthoFilter,
                 L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                 T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? L+1 : 3*L+1
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            sdwt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    isdwtall(xw[, wt])
    isdwtall(xw, wt, sm)

Computes the inverse stationary discrete wavelet transform (ISDWT) on each slice of signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: SDWT-transformed signal.
- `wt::OrthoFilter`: (Default: `nothing`) Orthogonal wavelet filter.
- `sm::Integer`: If `sm` is included as an argument, the `sm`-shifted inverse transform will
  be computed. This results in significantly faster computation, but fails to fully utilize
  the strength of redundant wavelet transforms.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SDWT on all signals in x
xw = sdwtall(x, wt)

# ISDWT on all signals in xw
x̂ = isdwtall(xw)
```

**See also:** [`isdwt`](@ref)
"""
function isdwtall(xw::AbstractArray{T}, wt::OrthoFilter) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            isdwt!(xᵢ, xwᵢ, wt)
        end
    end
    return x
end

function isdwtall(xw::AbstractArray{T}, wt::OrthoFilter, sm::Integer) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            isdwt!(xᵢ, xwᵢ, wt, sm)
        end
    end
    return x
end

# ----- SWPT on a set of signals -----
"""
    swptall(x, wt[, L])

Computes the stationary wavelet packet transform (SWPT) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SWPT on all signals in x
xw = swptall(x, wt)
```

**See also:** [`swpt`](@ref)
"""
function swptall(x::AbstractArray{T},
                 wt::OrthoFilter,
                 L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                 T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? 1<<L : Int(4^L)
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            swpt!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iswptall(xw[, wt])
    iswptall(xw, wt, sm)

Computes the inverse stationary wavelet packet transform (ISWPT) on each slice of signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: SWPT-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.
- `sm::Integer`: If `sm` is included as an argument, the `sm`-shifted inverse transform will
  be computed. This results in significantly faster computation, but fails to fully utilize
  the strength of redundant wavelet transforms.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SWPT on all signals in x
xw = swptall(x, wt)

# ISWPT on all signals in xw
x̂ = iswptall(xw)
```

**See also:** [`iswpt`](@ref)
"""
function iswptall(xw::AbstractArray{T}, wt::OrthoFilter) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iswpt!(xᵢ, xwᵢ, wt)
        end
    end
    return x
end

function iswptall(xw::AbstractArray{T}, wt::OrthoFilter, sm::Integer) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iswpt!(xᵢ, xwᵢ, wt, sm)
        end
    end
    return x
end

# ----- SWPD on a set of signals -----
"""
    swpdall(x, wt[, L])

Computes the stationary wavelet packet decomposition (SWPD) on each slice of signal.

# Arguments
- `x::AbstractArray{T} where T<:Number`: Input `N-1`-D signals, where each signal is sliced
  at dimension `N`.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `Wavelets.maxtransformlevels(xᵢ)`) Number of levels of wavelet
  transforms.

# Returns
- `::Array{T}`: Slices of transformed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SWPD on all signals in x
xw = swpdall(x, wt)
```

**See also:** [`swpd`](@ref)
"""
function swpdall(x::AbstractArray{T},
                  wt::OrthoFilter,
                  L::Integer = minimum(size(x)[1:end-1]) |> maxtransformlevels) where
                  T<:Number
    @assert 2 ≤ ndims(x) ≤ 3
    # Allocate space for transforms
    sz = size(x)[1:(end-1)]
    k = ndims(x)==2 ? 1<<(L+1)-1 : sum(4 .^(0:L))
    N = size(x)[end]
    xw = Array{T}(undef, (sz...,k,N))
    # Dimension to slice
    dim_xw = ndims(xw)
    dim_x = ndims(x)
    # Compute transforms
    @inbounds begin
        @views for (xwᵢ, xᵢ) in zip(eachslice(xw, dims=dim_xw), eachslice(x, dims=dim_x))
            swpd!(xwᵢ, xᵢ, wt, L)
        end
    end
    return xw
end

"""
    iswpdall(xw, wt, L, sm)
    iswpdall(xw, wt[, L])
    iswpdall(xw, wt, tree, sm)
    iswpdall(xw, wt, tree)

Computes the inverse autocorrelation wavelet packet decomposition (ISWPD) on each slice of
signal.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: SWPD-transformed signal.
- `wt::Union{OrthoFilter, Nothing}`: (Default: `nothing`) Orthogonal wavelet filter.
- `L::Integer`: (Default: `minimum(size(xw)[1:end-2]) |> maxtransformlevels`) Number of
  levels of wavelet transforms.
- `tree::BitVector`: Binary tree for inverse transform to be computed accordingly.
- `sm::Integer`: If `sm` is included as an argument, the `sm`-shifted inverse transform will
  be computed. This results in significantly faster computation, but fails to fully utilize
  the strength of redundant wavelet transforms.

# Returns
- `::Array{T}`: Slices of reconstructed signals.

# Examples
```julia
using Wavelets, WaveletsExt

# Generate random signals
x = randn(32, 5)
# Create wavelet
wt = wavelet(WT.db4)

# SWPD on all signals in x
xw = swpdall(x, wt)

# ISWPD on all signals in xw
x̂ = iswpdall(xw)
x̂ = iswpdalll(xw, maketree(x))
x̂ = iswpdall(xw, 5)
```

**See also:** [`iswpd`](@ref)
"""
function iswpdall(xw::AbstractArray{T},
                   wt::OrthoFilter,
                   L::Integer = minimum(size(xw)[1:end-2]) |> maxtransformlevels) where
                   T<:Number
    return iswpdall(xw, wt, maketree(size(xw)[1:(end-2)]..., L))
end

function iswpdall(xw::AbstractArray{T}, wt::OrthoFilter, L::Integer, sm::Integer) where
                  T<:Number
    return iswpdall(xw, wt, maketree(size(xw)[1:(end-2)]..., L), sm)
end

function iswpdall(xw::AbstractArray{T}, wt::OrthoFilter, tree::BitVector) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iswpd!(xᵢ, xwᵢ, wt, tree)
        end
    end
    return x
end

function iswpdall(xw::AbstractArray{T},
                  wt::OrthoFilter,
                  tree::BitVector,
                  sm::Integer) where T<:Number
    @assert 3 ≤ ndims(xw) ≤ 4
    # Allocate space for transforms
    sz = size(xw)[1:(end-2)]
    N = size(xw)[end]
    x = Array{T}(undef, (sz...,N))
    # Dimension to slice
    dim_x = ndims(x)
    dim_xw = ndims(xw)
    # Compute transforms
    @inbounds begin
        @views for (xᵢ, xwᵢ) in zip(eachslice(x, dims=dim_x), eachslice(xw, dims=dim_xw))
            iswpd!(xᵢ, xwᵢ, wt, tree, sm)
        end
    end
    return x
end


text:
module DataStructures

    using Base: HasEltype, HasLength, IteratorEltype, IteratorSize, SizeUnknown,
                lt, Ordering, ForwardOrdering, Forward, ReverseOrdering, Reverse, Lt,
                isbitsunion, isiterable, dict_with_eltype, KeySet, Callable, _tablesz,
                findnextnot, unsafe_getindex, unsafe_setindex!, peek

    using Compat # Provides Base.Order.ReverseOrdering(). May remove this line with julia 1.4
    using OrderedCollections
    using OrderedCollections: isordered
    export OrderedDict, OrderedSet, LittleDict
    export DefaultDict, DefaultOrderedDict

    export complement, complement!

    export Deque, Stack, Queue, CircularDeque
    export enqueue!, dequeue!, dequeue_pair!, update!, reverse_iter
    export capacity, num_blocks, top_with_handle, sizehint!

    export Accumulator, counter, reset!, inc!, dec!
    export IntDisjointSets, DisjointSets, num_groups, find_root!, in_same_set, root_union!
    export FenwickTree, length, inc!, dec!, incdec!, prefixsum

    export AbstractHeap, compare, extract_all!, extract_all_rev!
    export BinaryHeap, BinaryMinHeap, BinaryMaxHeap, nlargest, nsmallest
    export MutableBinaryHeap, MutableBinaryMinHeap, MutableBinaryMaxHeap
    export heapify!, heapify, heappop!, heappush!, isheap
    export BinaryMinMaxHeap, popmin!, popmax!, popall!

    export Trie, subtrie, keys_with_prefix, partial_path

    export LinkedList, Nil, Cons, nil, cons, head, tail, list, filter, cat,
           reverse
    export MutableLinkedList
    export SortedDict, SortedMultiDict, SortedSet
    export SDToken, SDSemiToken, SMDToken, SMDSemiToken
    export SetToken, SetSemiToken
    export startof
    export pastendsemitoken, beforestartsemitoken
    export searchsortedafter, searchequalrange
    export packcopy, packdeepcopy
    export exclusive, inclusive, semitokens
    export orderobject, ordtype, Lt, compare, onlysemitokens

    export MultiDict, enumerateall
    export RobinDict
    export OrderedRobinDict, isordered
    export SwissDict

    export DiBitVector

    export RBTree, search_node, minimum_node
    export AVLTree, sorted_rank

    export findkey

    include("delegate.jl")

    include("deque.jl")
    include("circ_deque.jl")
    include("stack.jl")
    include("queue.jl")
    include("accumulator.jl")
    include("disjoint_set.jl")
    include("heaps.jl")

    include("default_dict.jl")
    include("dict_support.jl")
    include("trie.jl")

    include("int_set.jl")

    include("fenwick.jl")

    include("list.jl")
    include("mutable_list.jl")
    include("balanced_tree.jl")
    include("tokens.jl")

    import .Tokens: IntSemiToken

    include("multi_dict.jl")
    include("sorted_dict.jl")
    include("sorted_multi_dict.jl")
    include("sorted_set.jl")
    include("tokens2.jl")
    include("container_loops.jl")
    include("robin_dict.jl")
    include("ordered_robin_dict.jl")
    include("swiss_dict.jl")
    export
        CircularBuffer,
        capacity,
        isfull
    include("circular_buffer.jl")

    export status
    export deref_key, deref_value, deref, advance, regress

    export PriorityQueue, peek

    include("priorityqueue.jl")
    include("sparse_int_set.jl")
    export SparseIntSet

    include("dibit_vector.jl")
    include("avl_tree.jl")
    include("red_black_tree.jl")

    include("deprecations.jl")
end


text:
using FluxReconstruction, OrdinaryDiffEq, Plots, LinearAlgebra, KitBase
using KitBase.ProgressMeter: @showprogress

ps0 = KitBase.PSpace2D(0.0, 1.0, 20, 0.0, 1.0, 20)
deg = 2
ps = FRPSpace2D(ps0, deg)
a = [0.1, 0.1]

u0 = zeros(ps.nx, ps.ny, deg + 1, deg + 1)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    u0[i, j, k, l] = max(
        exp(-300 * ((ps.xpg[i, j, k, l, 1] - 0.5)^2 + (ps.xpg[i, j, k, l, 2] - 0.5)^2)),
        1e-2,
    )
end

function dudt!(du, u, p, t)
    J, ll, lr, dhl, dhr, lpdm, ax, ay = p

    nx = size(u, 1)
    ny = size(u, 2)
    nsp = size(u, 3)

    f = zeros(nx, ny, nsp, nsp, 2)
    for i in axes(f, 1), j in axes(f, 2), k = 1:nsp, l = 1:nsp
        fg, gg = ax * u[i, j, k, l], ay * u[i, j, k, l]
        f[i, j, k, l, :] .= inv(J[i, j][k, l]) * [fg, gg]
    end

    u_face = zeros(nx, ny, 4, nsp)
    f_face = zeros(nx, ny, 4, nsp, 2)
    for i in axes(u_face, 1), j in axes(u_face, 2), l = 1:nsp
        u_face[i, j, 1, l] = dot(u[i, j, l, :], ll)
        u_face[i, j, 2, l] = dot(u[i, j, :, l], lr)
        u_face[i, j, 3, l] = dot(u[i, j, l, :], lr)
        u_face[i, j, 4, l] = dot(u[i, j, :, l], ll)

        for m = 1:2
            f_face[i, j, 1, l, m] = dot(f[i, j, l, :, m], ll)
            f_face[i, j, 2, l, m] = dot(f[i, j, :, l, m], lr)
            f_face[i, j, 3, l, m] = dot(f[i, j, l, :, m], lr)
            f_face[i, j, 4, l, m] = dot(f[i, j, :, l, m], ll)
        end
    end

    fx_interaction = zeros(nx + 1, ny, nsp)
    for i = 2:nx, j = 1:ny, k = 1:nsp
        au =
            (f_face[i, j, 4, k, 1] - f_face[i-1, j, 2, k, 1]) /
            (u_face[i, j, 4, k] - u_face[i-1, j, 2, k] + 1e-6)
        fx_interaction[i, j, k] = (
            0.5 * (f_face[i, j, 4, k, 1] + f_face[i-1, j, 2, k, 1]) -
            0.5 * abs(au) * (u_face[i, j, 4, k] - u_face[i-1, j, 2, k])
        )
    end
    fy_interaction = zeros(nx, ny + 1, nsp)
    for i = 1:nx, j = 2:ny, k = 1:nsp
        au =
            (f_face[i, j, 1, k, 2] - f_face[i, j-1, 3, k, 2]) /
            (u_face[i, j, 1, k] - u_face[i, j-1, 3, k] + 1e-6)
        fy_interaction[i, j, k] = (
            0.5 * (f_face[i, j, 1, k, 2] + f_face[i, j-1, 3, k, 2]) -
            0.5 * abs(au) * (u_face[i, j, 1, k] - u_face[i, j-1, 3, k])
        )
    end

    rhs1 = zeros(nx, ny, nsp, nsp)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp
        rhs1[i, j, k, l] = dot(f[i, j, :, l, 1], lpdm[k, :])
    end
    rhs2 = zeros(nx, ny, nsp, nsp)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp
        rhs2[i, j, k, l] = dot(f[i, j, k, :, 2], lpdm[l, :])
    end

    for i = 2:nx-1, j = 2:ny-1, k = 1:nsp, l = 1:nsp
        du[i, j, k, l] = -(
            rhs1[i, j, k, l] +
            rhs2[i, j, k, l] +
            (fx_interaction[i, j, l] - f_face[i, j, 4, l, 1]) * dhl[k] +
            (fx_interaction[i+1, j, l] - f_face[i, j, 2, l, 1]) * dhr[k] +
            (fy_interaction[i, j, k] - f_face[i, j, 1, k, 2]) * dhl[l] +
            (fy_interaction[i, j+1, k] - f_face[i, j, 3, k, 2]) * dhr[l]
        )
    end

    return nothing
end

tspan = (0.0, 1.0)
p = (ps.J, ps.ll, ps.lr, ps.dhl, ps.dhr, ps.dl, a[1], a[2])
prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.01
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Euler(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    step!(itg)
end

contourf(ps.x[:, 1], ps.y[1, :], itg.u[:, :, 2, 2])


text:
function getTableauLieA(T=Float64)
    a = Array{T}([ 1 ])
    b = Array{T}([ 0 ])
    TableauSplittingNS(:LieTrotterSplittingA, 1, a, b)
end

function getTableauLieB(T=Float64)
    a = Array{T}([ 0 ])
    b = Array{T}([ 1 ])
    TableauSplittingNS(:LieTrotterSplittingB, 1, a, b)
end

function getTableauStrang(T=Float64)
    a = Array{T}([ 1//2 ])
    b = Array{T}([ 1//2 ])
    TableauSplittingNS(:StrangSplitting, 2, a, b)
end

function getTableauMcLachlan2(T=Float64; α=0.1932)
    a = Array{T}([ α, 0.5 - α ])
    b = Array{T}([ 0.5 - α, α ])
    TableauSplittingNS(:McLachlanSplitting, 2, a, b)
end

function getTableauMcLachlan4(T=Float64)
    a = Array{T}(@dec128 [ (146 +  5*√19) / 540,
                           ( -2 + 10*√19) / 135,
                           1/5,
                           (-23 - 20*√19) / 270,
                           ( 14 -    √19) / 108])
    TableauSplittingNS(:McLachlanSplitting, 4, a, a[end:-1:1])
end

function getTableauTripleJump(T=Float64)
    fac = @dec128 2^(1/3)
    den = @dec128 1/(2-fac)
    a = Array{T}([ den, -fac*den ])
    TableauSplittingSS(:TripleJumpSplitting, 4, a)
end

function getTableauSuzukiFractal(T=Float64)
    fac = @dec128 4^(1/3)
    den = @dec128 1/(4-fac)
    a = Array{T}([ den, den, -fac*den ])
    TableauSplittingSS(:SuzukiFractalSplitting, 4, a)
end


text:
"""
  - `run(sim::Simulation, until::Event)`
  - `run(sim::Simulation, until::TimeType)`
  - `run(sim::Simulation, period::Period)`
  - `run(sim::Simulation, period::Number)`
  - `run(sim::Simulation)`

Executes [`step`](@ref) until the given criterion `until` is met:

- if it is not specified, the method will return when there are no further events to be processed
- if it is an `Event`, the method will continue stepping until this event has been triggered and will return its value
- if it is a `TimeType`, the method will continue stepping until the simulation’s time reaches until
- if it is a `Period`, the method will continue stepping until the simulation’s time has passed until periods
- if it is a `Number`, the method will continue stepping until the simulation’s time has passed until elementary periods

In the last two cases, the simulation can prematurely stop when there are no further events to be processed.
"""
function run(sim::Simulation, until::Event) :: Any
  append_callback(until, stop_simulation)
  try
    while true
      step(sim)
    end
  catch exc
    if isa(exc, StopSimulation)
      return exc.value
    else
      rethrow(exc)
    end
  end
end

function run(sim::Simulation, period::Period) :: Any
  run(sim, timeout(sim, period))
end

function run(sim::Simulation, period::Number) :: Any
  run(sim, timeout(sim, sim.granularity(period)))
end

function run{T}(sim::Simulation{T}, until::T) :: Any
  run(sim, timeout(sim, now(sim)-until))
end

function run(sim::Simulation) :: Any
  run(sim, timeout(sim, typemax(sim.granularity)))
end

function stop_simulation(sim::Simulation, ev::Event)
  throw(StopSimulation(ev.value))
end

"""
  `now(sim::Simulation) :: TimeType`

Returns the current simulation time.
"""
function now(sim::Simulation) :: TimeType
  return sim.time
end

"""
  - `schedule!(sim::Simulation, ev::Event, delay::Period; priority::Bool=false, value::Any=nothing) :: Event`
  - `schedule!(sim::Simulation, ev::Event, delay::Number=0; priority::Bool=false, value::Any=nothing) :: Event`

Schedules an event at time `sim.time + delay` with a `priority` and a `value`.

If the event is already scheduled, the key is updated with the new `delay` and `priority`. The new `value` is also set.

If the event is being processed, an [`EventProcessing`](@ref) exception is thrown.
"""
function schedule!(sim::Simulation, ev::Event, delay::Period; priority::Bool=false, value::Any=nothing) :: Event
  if ev.state == processed
    throw(EventProcessed())
  end
  ev.value = value
  if ev.state == triggered
    id = sim.heap[ev].id
  else
    id = sim.sid+=one(UInt)
  end
  sim.heap[ev] = EventKey(sim.time + delay, priority, id)
  ev.state = triggered
  return ev
end

function schedule!(sim::Simulation, ev::Event, delay::Number=0; priority::Bool=false, value::Any=nothing) :: Event
  schedule!(sim, ev, sim.granularity(delay), priority=priority, value=value)
end

"""
  - `schedule(sim::Simulation, ev::Event, delay::Period; priority::Bool=false, value::Any=nothing) :: Event`
  - `schedule(sim::Simulation, ev::Event, delay::Number=0; priority::Bool=false, value::Any=nothing) :: Event`

Schedules an event at time `sim.time + delay` with a `priority` and a `value`.

If the event is already scheduled or is being processed, an [`EventNotIdle`](@ref) exception is thrown.
"""
function schedule(sim::Simulation, ev::Event, delay::Period; priority::Bool=false, value::Any=nothing) :: Event
  if !(ev.state == idle)
    throw(EventNotIdle())
  end
  ev.value = value
  sim.heap[ev] = EventKey(sim.time + delay, priority, sim.sid+=one(UInt))
  ev.state = triggered
  return ev
end

function schedule(sim::Simulation, ev::Event, delay::Number=0; priority::Bool=false, value::Any=nothing) :: Event
  schedule(sim, ev, sim.granularity(delay), priority=priority, value=value)
end

"""
  `step(sim::Simulation) :: Bool`

Does a simulation step and processes the next event.

Only used internally.
"""
function step(sim::Simulation)
  if isempty(sim.heap)
    throw(EmptySchedule())
  end
  (ev, key) = peek(sim.heap)
  dequeue!(sim.heap)
  sim.time = key.time
  ev.state = processed
  while !isempty(ev.callbacks)
    cb = dequeue!(ev.callbacks)
    cb(sim, ev)
  end
end


text:
abstract type AbstractCorrelator{T, TA} end
"""
$(SIGNATURES)

EarlyPromptLateCorrelator holding a user defined number of correlation values.
"""
struct EarlyPromptLateCorrelator{T, TA <: AbstractVector{T}} <: AbstractCorrelator{T, TA}
    accumulators::TA
    EarlyPromptLateCorrelator{T, TA}(accumulators) where {T, TA <: AbstractVector{T}} = length(accumulators) < 3 ?
        error("Early-Prompt-Late-Correlator needs at least 3 accumulators") : new{T, TA}(accumulators)
end

function EarlyPromptLateCorrelator(accumulators::AbstractVector)
    EarlyPromptLateCorrelator{typeof(first(accumulators)), typeof(accumulators)}(accumulators)
end

type_for_num_ants(num_ants::NumAnts{1}) = ComplexF64
type_for_num_ants(num_ants::NumAnts{N}) where N = SVector{N, ComplexF64}

"""
$(SIGNATURES)

EarlyPromptLateCorrelator constructor without parameters and some default parameters.
"""
function EarlyPromptLateCorrelator(
    num_ants::NumAnts = NumAnts(1),
    num_accumulators::NumAccumulators{M} = NumAccumulators(3),
) where M
    EarlyPromptLateCorrelator(
        zero(SVector{M, type_for_num_ants(num_ants)})
    )
end

"""
$(SIGNATURES)

EarlyPromptLateCorrelator constructor for large number of accumulators.
"""
function EarlyPromptLateCorrelator(
    num_ants::NumAnts,
    num_accumulators::Integer
)
    EarlyPromptLateCorrelator(
        [zero(type_for_num_ants(num_ants)) for i = 1:num_accumulators]
    )
end

"""
$(SIGNATURES)

Get number of antennas from correlator
"""
get_num_ants(correlator::AbstractCorrelator{Complex{T}}) where {T} = 1
get_num_ants(correlator::AbstractCorrelator{SVector{N,T}}) where {N,T} = N

"""
$(SIGNATURES)

Get number of accumulators
"""
get_num_accumulators(correlator::AbstractCorrelator) = size(correlator.accumulators, 1)

"""
$(SIGNATURES)

Get early correlator index
"""
get_early_index(correlator_sample_shifts, early_late_index_shift) =
    get_prompt_index(correlator_sample_shifts) + early_late_index_shift

"""
$(SIGNATURES)

Get prompt correlator index
"""
get_prompt_index(correlator_sample_shifts) = findfirst(iszero, correlator_sample_shifts)

"""
$(SIGNATURES)

Get late correlator index
"""
get_late_index(correlator_sample_shifts, early_late_index_shift) =
    get_prompt_index(correlator_sample_shifts) - early_late_index_shift

"""
$(SIGNATURES)

Get all correlator accumulators
"""
get_accumulators(correlator::AbstractCorrelator) = correlator.accumulators

"""
$(SIGNATURES)

Get a specific accumulator with `index` counted negative for late and
positive for early accumulators.
"""
function get_accumulator(correlator::AbstractCorrelator, correlator_sample_shifts, index::Integer)
    correlator.accumulators[index + get_prompt_index(correlator_sample_shifts)]
end

"""
$(SIGNATURES)

Get early correlator
"""
function get_early(
    correlator::AbstractCorrelator,
    correlator_sample_shifts,
    early_late_index_shift
)
    correlator.accumulators[
        get_early_index(correlator_sample_shifts, early_late_index_shift)
    ]
end

"""
$(SIGNATURES)

Get prompt correlator
"""
function get_prompt(correlator::AbstractCorrelator, correlator_sample_shifts)
    correlator.accumulators[get_prompt_index(correlator_sample_shifts)]
end

"""
$(SIGNATURES)

Get late correlator
"""
function get_late(
    correlator::AbstractCorrelator,
    correlator_sample_shifts,
    early_late_index_shift
)
    correlator.accumulators[
        get_late_index(correlator_sample_shifts, early_late_index_shift)
    ]
end

"""
$(SIGNATURES)

Zero the Correlator
"""
function zero(correlator::T) where T <: AbstractCorrelator
    T(zero(correlator.accumulators))
end

"""
$(SIGNATURES)

Filter the Correlator by the function `post_corr_filter`
"""
function filter(post_corr_filter, correlator::T) where T <: AbstractCorrelator
    (T.name.wrapper)(map(x -> post_corr_filter(x), get_accumulators(correlator)))
end

"""
$(SIGNATURES)

Calculate the replica phase offset required for the correlator with
respect to the prompt correlator, expressed in samples. The shifts are
ordered from latest to earliest replica.
"""
function get_correlator_sample_shifts(
    system::AbstractGNSS,
    correlator::AbstractCorrelator{T, <:SVector{M}},
    sampling_frequency,
    preferred_code_shift
) where {T,M}
    num_corrs = floor(Int, M / 2)
    sample_shift = preferred_code_shift_to_sample_shift(
        preferred_code_shift,
        sampling_frequency,
        system
    )
    SVector{M}(-num_corrs:num_corrs) .* sample_shift
end
function get_correlator_sample_shifts(
    system::AbstractGNSS,
    correlator::AbstractCorrelator,
    sampling_frequency,
    preferred_code_shift
)
    num_corrs = floor(Int, get_num_accumulators(correlator) / 2)
    sample_shift = preferred_code_shift_to_sample_shift(
        preferred_code_shift,
        sampling_frequency,
        system
    )
    (-num_corrs:num_corrs) .* sample_shift
end
function preferred_code_shift_to_sample_shift(
    preferred_code_shift,
    sampling_frequency,
    system
)
    sample_shift = round(Int, preferred_code_shift * sampling_frequency / get_code_frequency(system))
    max(1, sample_shift)
end

"""
$(SIGNATURES)

Calculates the index for the early and late samples
"""
function get_early_late_index_shift(
    system,
    correlator_sample_shifts,
    correlator,
    sampling_frequency,
    preferred_code_shift
)
    preferred_code_phase_distance = Inf
    early_late_index_shift = 1
    for i = get_prompt_index(correlator_sample_shifts) + 1:length(correlator_sample_shifts)
        code_phase_shift = correlator_sample_shifts[i] * get_code_frequency(system) / sampling_frequency
        preferred_code_phase_distance_temp = abs(code_phase_shift - preferred_code_shift)
        if preferred_code_phase_distance_temp < preferred_code_phase_distance
            early_late_index_shift = i - get_prompt_index(correlator_sample_shifts)
            preferred_code_phase_distance = preferred_code_phase_distance_temp
        end
    end
    early_late_index_shift
end

"""
$(SIGNATURES)

Calculate the total spacing between early and late correlator in samples.
"""
function get_early_late_sample_spacing(
    correlator_sample_shifts,
    early_late_index_shift
)
    correlator_sample_shifts[get_early_index(correlator_sample_shifts, early_late_index_shift)] -
    correlator_sample_shifts[get_late_index(correlator_sample_shifts, early_late_index_shift)]
end

"""
$(SIGNATURES)

Normalize the correlator
"""
function normalize(correlator::AbstractCorrelator, integrated_samples)
    filter(x -> x / integrated_samples, correlator)
end
"""
$(SIGNATURES)

Perform a correlation for multi antenna systems
"""
function correlate(
    correlator::T,
    downconverted_signal::AbstractVector,
    code,
    correlator_sample_shifts,
    start_sample,
    num_samples
) where {T <: AbstractCorrelator}
    accumulators = zero_accumulators(get_accumulators(correlator), downconverted_signal)
    @inbounds @fastmath for i = start_sample:num_samples + start_sample - 1
        for j = 1:length(accumulators)
            sample_shift = correlator_sample_shifts[j] - correlator_sample_shifts[1]
            accumulators[j] += downconverted_signal[i] * code[i + sample_shift]
        end
    end

    T(map(+, get_accumulators(correlator), accumulators))
end

function zero_accumulators(accumulators::SVector, signal)
    zeros(MVector{length(accumulators), eltype(signal)})
end
function zero_accumulators(accumulators::Vector, signal)
    zeros(eltype(signal), length(accumulators))
end

"""
$(SIGNATURES)

Perform a correlation for multi antenna systems
"""
function correlate(
    correlator::AbstractCorrelator{<: SVector{N}},
    downconverted_signal::AbstractMatrix,
    code,
    correlator_sample_shifts,
    start_sample,
    num_samples,
) where {N}

    accumulators = map(correlator_sample_shifts) do correlator_sample_shift
        correlate_single_tap(
            NumAnts(N),
            correlator_sample_shift - correlator_sample_shifts[1],
            start_sample,
            num_samples,
            downconverted_signal,
            code
        )
    end

    typeof(correlator)(map(+, get_accumulators(correlator), accumulators))
end

function correlate_single_tap(
    ::NumAnts{N},
    offset,
    start_sample,
    num_samples,
    downconverted_signal,
    code
) where N
    accumulator = zero(MVector{N, eltype(downconverted_signal)})
    @inbounds @fastmath for i = start_sample:num_samples + start_sample - 1, j = 1:length(accumulator)
        accumulator[j] += downconverted_signal[i,j] * code[i + offset]
    end
    SVector(accumulator)
end


text:
srand(1)

# --- generate a random (n,ne,np) network.
# --- d = matrix diagonal
# --- ne = number of edges
# --- np = number of positive edges
function randomnkp (d,ne,np)
    n = length(d)                       # int matrix diagonal - never change
    nn = ne-np                          # negative edges
    nz = n*(n-1)-ne                     # number of zero-edges
    intlist = [ones(Int64,np);-ones(Int64,nn);zeros(Int64,nz)]
    rect = reshape(shuffle(intlist),n,n-1)
    rrect = [rect zeros(Int64,n)]       # now insert back the diagonal
    lrect = [zeros(Int64,n) rect]
    interrand = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
    return interrand
end

# --- return a (n,ne,np) preserved shuffled version of a network
function shufflenet (intr,slf)
    n = length(slf)
    nw = similar(intr)                  # this will be the randomized version
    shf = shuffle(slf)                  # keep number of self degredating genes
    tl = tril(intr)                     # remove diag and join upper/lower triang
    tu = triu(intr)
    rect = tu[:,2:end]+tl[:,1:end-1]    # shuffle this
    rect = reshape(shuffle(reshape(rect,length(rect),1)[:,1]),n,n-1)
    rrect = [rect zeros(Int64,n)]       # now insert back the shuffled diagonal
    lrect = [zeros(Int64,n) rect]
    nw = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
    return (nw,shf)
end

# --- generate the symmetric adjacency matrix with elements 0,1 ---
function symadj(m)
    n = length(m[:,1])
    ad = zeros(Int64,n,n)
    for (i=1:n) for (j=i+1:n)       # adjacency matrix - without diagonals
        ad[i,j] = m[i,j]==0 && m[j,i]==0 ? 0 : 1
        ad[j,i] = ad[i,j]
    end end
    return ad
end

## --- check connectivity -----------
function isconnected (adj)
    n = length(adj[:,1])                # number of nodes
    nsq = n*n
    dum = adj
    (sumprev,sumdum) = (0,0)
    for (i=1:n)
        dum = sign(dum + dum*adj)
        sumdum = sum(dum)
        if (sumdum == sumprev)
            break
        end
        sumprev = sumdum
    end
    return (sumdum==nsq)
end

## --- node update function (integrator): yeast's majority rule -----------
function upd(s,intr,slf,n)
    sss = similar(s)
    inp = intr*s
    for (i=1:n)
        sss[i] = (inp[i]>0)?1:(inp[i]<0)?0:(slf[i]==1)?0:s[i]
    end
    return sss
end

## ---- output is the attractor list + basin sizes ---
## ---- calculated from a subset of initial conditions ---
function getcyclessample(intr,slf,ninit)

    n = length(slf)
    attrlist = zeros(Int64,n)           # s=0 is always an attractor
    keyzeros = join(zeros(Int64,n))
    attdict = [keyzeros=>0]    # values are the basin sizes
    for (i=1:ninit)       # calculate basins using ninit initial states
        s = rand(0:1,n)
        skey = join(s)
        pathdict = Dict()
        while (!haskey(attdict,skey) && !haskey(pathdict,skey))
            pathdict[skey] = "-"
            s = upd(s,intr,slf,n)
            skey=join(s)
        end
        if (haskey(attdict,skey)) attdict[skey] += 1
        else                        # this is a new cycle or fp
            attdict[skey] = 1
            attrlist = [attrlist s]
        end
    end
    basinlist = [attdict[join(attrlist[:,i])]/ninit for i in 1:length(attrlist[1,:])]
    return attrlist,basinlist
end

## ---- find number of coop nodes in a state ---
function ncoop (s,intr)
    n = length(s)
    act = intr*diagm(s) # keep only the columns for active nodes
    nc=0
    for (i=1:n)
#        if (abs(sum(act[i,:])) == sum(abs(act[i,:])))
        if (!in(1,act[i,:]) || !in(-1,act[i,:]))
            nc += 1
        end
    end
    return nc
end

## ---- find coop coef for sampled cycles of a network ---
function coopattr (intr,slf,ninit)
    n = length(slf)
    (listatt,listb) = getcyclessample(intr,slf,ninit)
    listperiod=Int64[];
    listcop=Float64[];
    listact=Float64[];
    npointattr = 0;
    for (i=1:length(listatt[1,:]))        # go over each attractor
        s = listatt[:,i]
#        if (sum(s) != 0)
        ncp = ncoop(s,intr)
#            nac = sum(abs(intr)*s)
        nac = sum(s)
        period = 1
        s=upd(s,intr,slf,n)
        while (s != listatt[:,i])
            ncp += ncoop(s,intr)
#            nac += sum(abs(intr)*s)
            nac += sum(s)
            period += 1
            s=upd(s,intr,slf,n)
        end
        push!(listperiod,period)
        push!(listcop,ncp/period)
        push!(listact,nac/period)
        if (period==1) npointattr += 1 end
    #        end
    end
#    if (length(listperiod)==1)  return -1,-1,-1,-1 # no fixed point other than (0,..0)
#    else                    # also returns the fracion of active nodes
#        return (1.0*mean(listcop)/n,mean(listact)/n,length(listc[1,:]),npointattr)
    return listatt,listb,listperiod,listcop,listact
#    end
end

# returns the node-id permuted version of the square matrix a
function motifpermute (a,perm)
    b = similar(a)
    n = length(a[:,1])
    for i=1:n
        for j=1:n
            b[i,j]=a[perm[i],perm[j]]
        end
    end
    return b
end

# returns the dumb id of a motif matrix (3x3)
function dumbid (a)
    id=0;
    for i=1:3
        for j=1:3
            digit = 3*(i-1)+(j-1) # (1,1)=0, (1,2)=1,..,(2,1)=3..
            id += (a[i,j]+1)*(3^digit)
        end
    end
    return id
end

function motifmap ()
    a=zeros(Int64,3,3)
    motifid = Dict()                    # values are the true ids
    motifidrev = Dict()                 # values are the dumb ids
    idx = 1 # this will be used as an array index, so I don't start from 0.
    for a[1,1]=-1:1
        for a[1,2]=-1:1
            for a[1,3]=-1:1
                for a[2,1]=-1:1
                    for a[2,2]=-1:1
                        for a[2,3]=-1:1
                            for a[3,1]=-1:1
                                for a[3,2]=-1:1
                                    for a[3,3]=-1:1
                                        perm = [1,2,3]
                                        id = dumbid(motifpermute(a,perm))
                                        if (!haskey(motifid,id))
                                            motifidrev[idx]=id
                                            for perm in collect(permutations([1:3]))
                                                motifid[dumbid(motifpermute(a,perm))]=idx
                                            end
                                            idx += 1
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return motifid,motifidrev
end

function motifmapnoloop ()
    a=zeros(Int64,3,3)
    motifid = Dict()                    # values are the true ids
    motifidrev = Dict()                 # values are the dumb ids
    idx=1 # this will be used as an array index, so I don't start from 0.
    for a[1,2]=-1:1
        for a[1,3]=-1:1
            for a[2,1]=-1:1
                for a[2,3]=-1:1
                    for a[3,1]=-1:1
                        for a[3,2]=-1:1
                            perm = [1,2,3]
                            id = dumbid(motifpermute(a,perm))
                            if (!haskey(motifid,id))
                                motifidrev[idx]=id
                                for perm in collect(permutations([1:3]))
                                    motifid[dumbid(motifpermute(a,perm))]=idx
                                end
                                idx += 1
                            end
                        end
                    end
                end
            end
        end
    end
    return motifid,motifidrev
end

function getmotifstats(nw,map,mc)

    n=length(nw[:,1])
    a=zeros(Int64,3,3)
    for i=1:n
        for j=i+1:n
            for k=j+1:n
                a = nw[[i,j,k],[i,j,k]]
                mc[map[dumbid(a)]] += 1
            end
        end
    end
end

# --- calculate average cooperativity over random networks
n = 11
ntotal = 1000000;           # ensemble size
numinit = 2^11              # num of initial cond for fp determination

# generate motif id maps
(motifid,motifidrev) = motifmap();

a = zeros(Int64,0)
for x in values(motifid) push!(a,x) end
## motifcount=zeros(Int64,sort(a)[length(a)])
## motifcountperf=zeros(Int64,sort(a)[length(a)])
## motifcountbio=zeros(Int64,sort(a)[length(a)])

(motifidnoloop,motifidnolooprev) = motifmapnoloop();

a = zeros(Int64,0)
for x in values(motifidnoloop) push!(a,x) end
## motifcountnoloop=zeros(Int64,sort(a)[length(a)])
## motifcountnoloopperf=zeros(Int64,sort(a)[length(a)])
## motifcountnoloopbio=zeros(Int64,sort(a)[length(a)])

## ## for f=["a.thaliana","d.melanogaster","lymphocyte","mammalian","myeloid","s.cerevisiae","s.pombe"]
## for f=["a.thaliana","lymphocyte","mammalian","myeloid","s.cerevisiae","s.pombe"]
##     path="7org_nw"
##     nw = readdlm("$path/$f");
##     getmotifstats(nw,motifid,motifcountbio)
##     intersh = nw-diagm(diag(nw));
##     getmotifstats(intersh,motifidnoloop,motifcountnoloopbio)
## end
## for ne=29:29
##     for np=15:15
##         ## dirname = "$n.$ne.$np"
##         ## if (!isdir(dirname))
##         ##     run(`mkdir $dirname`)
##         ## end
##         nself = 5
##         selfsh = [ones(Int64,nself),zeros(Int64,n-nself)]
##         listc=zeros(Float64,ntotal)    # list of coherence values
##         lista=zeros(Float64,ntotal)    # list of active node fractions
##         listl=zeros(Float64,ntotal) # list of fraction of point attractors
##         intersh = randomnkp(selfsh,ne,np) # generate one sample with given (nkp)
##         i=0
##         while (i<ntotal)
##             (intersh,selfsh) = shufflenet(intersh,selfsh) # shuffle
##             if (isconnected(symadj(intersh))) # check if connected
##                 (attlist,basinlist,listperiod,listcop,listact) = coopattr(intersh,selfsh,numinit)
##                 if (length(attlist) > 1)
##                     i += 1              # increase network count
##                     nw = intersh - diagm(selfsh)
##                     ## edges = zeros(Int64,ne+sum(selfsh),3)
##                     ## idx=1
##                     ## for j=1:n
##                     ##     for k=1:n
##                     ##         if nw[j,k] != 0
##                     ##             edges[idx,:] = [j k nw[j,k]]
##                     ##             idx += 1
##                     ##         end end end
##                     ## writedlm("$dirname/$i.dat",edges,' ')
##                     ## writedlm("$dirname/$i.att",[attlist' basinlist listperiod listcop/n listact],' ')
##                     getmotifstats(nw,motifid,motifcount)
##                     getmotifstats(intersh,motifidnoloop,motifcountnoloop)
##                     if (mean(listcop/n)==1) # perfect coherence
##                         print("$i\n")
##                         getmotifstats(nw,motifid,motifcountperf)
##                         getmotifstats(intersh,motifidnoloop,motifcountnoloopperf)
##                     end
##                 end
##             end
##         end
##     end
## end

max_stars_repo_name: akabakcioglu/network-motifs

max_issues_repo_name: akabakcioglu/network-motifs

max_forks_repo_name: akabakcioglu/network-motifs


text:
@doc raw"""
    usinc(θ::Real)

Unnormalized version of `sinc` function, i.e. $\operatorname{usinc}(θ) = \frac{\sin(θ)}{θ}$.
This is equivalent to `sinc(θ/π)`.
"""
@inline usinc(θ::Real) = θ == 0 ? one(θ) : isinf(θ) ? zero(θ) : sin(θ) / θ

@doc raw"""
    usinc_from_cos(x::Real)

Unnormalized version of `sinc` function, i.e. $\operatorname{usinc}(θ) = \frac{\sin(θ)}{θ}$,
computed from $x = cos(θ)$.
"""
@inline function usinc_from_cos(x::Real)
    return if x >= 1
        one(x)
    elseif x <= -1
        zero(x)
    else
        sqrt(1 - x^2) / acos(x)
    end
end

@doc raw"""
    nzsign(z[, absz])

Compute a modified `sign(z)` that is always nonzero, i.e. where
````math
\operatorname(nzsign)(z) = \begin{cases}
    1 & \text{if } z = 0\\
    \frac{z}{|z|} & \text{otherwise}
\end{cases}
````
"""
@inline function nzsign(z, absz=abs(z))
    psignz = z / absz
    return ifelse(iszero(absz), one(psignz), psignz)
end

allocate(p, s::Size{S}) where {S} = similar(p, S...)
allocate(p::StaticArray, s::Size{S}) where {S} = similar(p, maybesize(s))
allocate(p, ::Type{T}, s::Size{S}) where {S,T} = similar(p, T, S...)
allocate(p::StaticArray, ::Type{T}, s::Size{S}) where {S,T} = similar(p, T, maybesize(s))

"""
    eigen_safe(x)

Compute the eigendecomposition of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before the decomposition.
"""
@inline eigen_safe(x; kwargs...) = eigen(x; kwargs...)
@inline function eigen_safe(x::StaticMatrix; kwargs...)
    s = size(x)
    E = eigen!(Matrix(parent(x)); kwargs...)
    return Eigen(SizedVector{s[1]}(E.values), SizedMatrix{s...}(E.vectors))
end

"""
    log_safe(x)

Compute the matrix logarithm of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before computing the log.
"""
@inline log_safe(x) = log(x)
@inline function log_safe(x::StaticMatrix)
    s = Size(x)
    return SizedMatrix{s[1],s[2]}(log(Matrix(parent(x))))
end

@generated maybesize(s::Size{S}) where {S} = prod(S) > 100 ? S : :(s)

"""
    select_from_tuple(t::NTuple{N, Any}, positions::Val{P})

Selects elements of tuple `t` at positions specified by the second argument.
For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns
`("c", "a", "a")`.
"""
@generated function select_from_tuple(t::NTuple{N,Any}, positions::Val{P}) where {N,P}
    for k in P
        (k < 0 || k > N) && error("positions must be between 1 and $N")
    end
    return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...)
end

@doc raw"""
    vec2skew!(X, v, k)

create a skew symmetric matrix inplace in `X` of size $k\times k$ from a vector `v`,
for example for `v=[1,2,3]` and `k=3` this
yields
````julia
[  0  1  2;
  -1  0  3;
  -2 -3  0
]
````
"""
function vec2skew!(X, v)
    k = size(X)[1]
    size(X)[2] != k && error("X is of wrong size, expected ($k,$k) got $(size(X)).")
    n = div(k * (k - 1), 2)
    length(v) < n && error("The vector $(v) is too short, expected $(n) got $(length(v)).")
    m = 0
    X .= [i < j ? (m += 1; v[m]) : zero(eltype(v)) for i in 1:k, j in 1:k]
    X .= X - X'
    return X
end
function vec2skew(v, k)
    X = similar(v, k, k)
    vec2skew!(X, v)
    return X
end

"""
    ziptuples(a, b[, c[, d[, e]]])

Zips tuples `a`, `b`, and remaining in a fast, type-stable way. If they have different
lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}) where {N,M}
    ex = Expr(:tuple)
    for i in 1:min(N, M)
        push!(ex.args, :((a[$i], b[$i])))
    end
    return ex
end
@generated function ziptuples(
    a::NTuple{N,Any},
    b::NTuple{M,Any},
    c::NTuple{L,Any},
) where {N,M,L}
    ex = Expr(:tuple)
    for i in 1:min(N, M, L)
        push!(ex.args, :((a[$i], b[$i], c[$i])))
    end
    return ex
end
@generated function ziptuples(
    a::NTuple{N,Any},
    b::NTuple{M,Any},
    c::NTuple{L,Any},
    d::NTuple{K,Any},
) where {N,M,L,K}
    ex = Expr(:tuple)
    for i in 1:min(N, M, L, K)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i])))
    end
    return ex
end
@generated function ziptuples(
    a::NTuple{N,Any},
    b::NTuple{M,Any},
    c::NTuple{L,Any},
    d::NTuple{K,Any},
    e::NTuple{J,Any},
) where {N,M,L,K,J}
    ex = Expr(:tuple)
    for i in 1:min(N, M, L, K, J)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i], e[$i])))
    end
    return ex
end


text:
# Helper functions to help set material properties are defined here.

import Colors: Color, hex
export meshcolor, wireframe, edges, visible, lambert, basic, phong, normal

function meshcolor(color::Color)
    color = string("#" * hex(color))
    :color, color
end

function wireframe(wireframe::Bool=true)
    :wireframe, wireframe
end

function lambert()
    :kind, "lambert"
end

function basic()
    :kind, "basic"
end

function phong()
    :kind, "phong"
end

function normal()
    :kind, "normal"
end

function visible(visible::Bool=true)
    :hidden, !visible
end

function edges(edgeColor::Color)
    color = string("#" * hex(color))
    :edges, color
end

# To allow edges to be set to false.
function edges(edges::Bool=true)
    :edges, edges
end


text:
###img2col

function img2col(A::AbstractArray{T,3})::AbstractArray{T,2} where {T}
    S = size(A)
    vs = prod(S[1:end-1])
    m = S[end]
    return reshape(A, vs, m)
end #function img2col(A::Array{T,3})

function img2col(A::AbstractArray{T,4})::AbstractArray{T,2} where {T}
    S = size(A)
    vs = prod(S[1:end-1])
    m = S[end]
    return reshape(permutedims(A, [2, 1, 3, 4]), vs, m)
end #function img2col(A::Array{T,4})

function img2row(A::AbstractArray{T,4})::AbstractArray{T,2} where {T}
    S = size(A)
    vs = prod(S[1:end-1])
    m = S[end]
    return reshape(A, vs, m)
end #function img2col(A::Array{T,4})

function img2col(A::AbstractArray{T,5})::AbstractArray{T,2} where {T}
    S = size(A)
    vs = prod(S[1:end-1])
    m = S[end]
    return reshape(permutedims(A, [2, 1, 3, 4, 5]), vs, m)
end #function img2col(A::Array{T,5})

export img2col

### col2img

function col2img(
    Av::AbstractArray{T,2},
    outputS::Tuple{Integer,Integer,Integer},
)::AbstractArray{T,3} where {T}
    return reshape(Av, outputS)
end #function col2img(Av::Array{T,2}, c::Integer)

function col2img(
    Av::AbstractArray{T,2},
    outputS::Tuple{Integer,Integer,Integer,Integer},
)::AbstractArray{T,4} where {T}
    H, W, c, m = outputS
    return permutedims(reshape(Av, (W, H, c, m)), [2, 1, 3, 4])
end #function col2img(Av::Array{T,2}, c::Integer; H::Integer=-1, W::Integer=-1)

function col2img(
    Av::AbstractArray{T,2},
    outputS::Tuple{Integer,Integer,Integer,Integer,Integer},
)::AbstractArray{T,5} where {T}
    H, W, D, c, m = outputS
    return permutedims(reshape(Av, (W, H, D, c, m)), [2, 1, 3, 4, 5])
end #function col2img(Av::Array{T,2}, c::Integer; H::Integer=-1, W::Integer=-1, D::Integer=-1)

export col2img1D, col2img2D, col2img3D